C 00000150 C PROGRAM SOLVE COMPUTES THE SOLUTION VECTOR, INVERSE OF THE NORMAL 00000160 C EQUATIONS, RESIDUALS, VARIANCE FACTOR, ABSOLUTE AND RELATIVE ERROR 00000170 C ELLIPSES FOR A PARAMETRIC LEAST SQUARES ADJUSTMENT OF A HORIZONTAL 00000180 C GEODETIC NETWORK (ON CLARKE 1866 ELLIPSOID). A-PRIORI INFORMATION 00000190 C FOR ANY NUMBER OF COORDINATES MAY BE ENTERED. THE PRESENT PROGRAM 00000200 C IS A MODIFIED VERSION OF PROGRAM 'LSAOMP' WRITTEN BY P.A. STEEVES. 00000210 C 00000220 C FOLLOWING IS A DESCRIPTION OF THE JCL, DATA INPUT AND DIMENSIONING 00000230 C INFORMATION REQUIRED. 00000240 C 00000250 C 1A INPUT DATA SETS 00000260 C 00000270 C OBSERVATION EQUATIONS (OUTPUT OF PROGRAM FORM3) 00000280 C 00000290 C //GO.FTNNF001 DD DSN=OBS.EQN, 00000300 C // VOL=SER=SEGEOM,UNIT=M2314,DISP=(OLD,KEEP) 00000310 C 00000320 C APPROXIMATE COORDINATES (OUTPUT OF PROGRAM FORM3) 00000330 C 00000340 C //GO.FTNNF001 DD DSN=APX.COORD, 00000350 C // VOL=SER=SEGEOM,UNIT=M2314,DISP=(OLD,KEEP) 00000360 C 00000370 C 1B OUTPUT DATA SETS 00000380 C 00000390 C SOLUTION VECTOR (INPUT TO PROGRAM COORDUP) 00000400 C 00000410 C //GO.FTNNF001 DD DSN=SOLN.VEC, 00000420 C // VOL=SER=SEGEOM,UNIT=M2314,DISP=(NEW,KEEP), 00000430 C // SPACE=(TRK,(1,2),RLSE),DCB=(RECFM=VBS,LRECL=8,BLKSIZE=7288) 00000440 C 00000450 C ADJUSTED COORDINATES (USED BY PROGRAM COMPARE(1)) 00000460 C 00000470 C //GO.FTNNF001 DD DSN=ADJ.COORD, 00000480 C // VOL=SER=SEGEOM,UNIT=M2314,DISP=(NEW,KEEP), 00000490 C // SPACE=(TRK,(1,2),RLSE),DCB=(RECFM=VBS,LRECL=20,BLKSIZE=7276) 00000500 C 00000510 C RESIDUALS,WEIGHTS (USED BY PROGRAM ALSAGN) 00000520 C // SPACE=(TRK,(3,24),RLSE),DCB=(RECFM=VBS,LRECL=32,BLKSIZE=7276) 00000560 C 00000530 C //GO.FTNNF001 DD DSN=RES.WTS, 00000540 C // VOL=SER=SEGEOM,UNIT=M2314,DISP=(NEW,KEEP), 00000550 C 00000570 C NORMAL EQUATIONS (INVERSE) 00000580 C 00000590 C //GO.FTNNF001 DD DSN=NORM.EQN, 00000600 C // VOL=SER=SEGEOM,UNIT=M2314,DISP=(NEW,KEEP), 00000610 C // SPACE=(TRK,(30,370),RLSE),DCB=(RECFM=VBS,LRECL=1816,BLKSIZE=728400000620 C 00000630 C SOLUTION SUMMARY (USED BY PROGRAM ALSAGN) 00000640 C 00000650 C //GO.FTNNF001 DD DSN=SUMMARY, 00000660 C // VOL=SER=SEGEOM,UNIT=M2314,DISP=(NEW,KEEP), 00000670 C // SPACE=(TRK,(1)),DCB=(RECFM=VBS,LRECL=40,BLKSIZE=7294) 00000680 C 00000690 C NOTE THE DATA SET NAME IS VARIABLE 00000700 C THE 'NN' IS REPLACED BY A VALID UNIT NUMBER (SEE DATA CARD #2 00000710 C BELOW) 00000720 C 00000730 C 00000740 C 2 CARD INPUT 00000750 C 00000760 C 2.1 JOB DESCRIPTION FORMAT 20A4 00000770 C 00000780 C 2.2 DATA SET UNIT NUMBERS FORMAT 7I2 00000790 C IN ORDER :OBS.EQN 00000800 C APX.COORD 00000810 C SOLN.VEC 00000820 C ADJ.COORD 00000830 C RES.WTS 00000840 C NORM.EQN 00000850 C SUMMARY 00000860 C 00000870 C 2.3 SEQUENCE NUMBERS OF WEIGHTED POINTS FORMAT 18I5 00000880 C 00000890 C 2.4 UPPER TRIANGULAR PORTION OF PX MATRIX FORMAT 8F10.0 00000900 C 00000910 C 2.5 NUMBER OF RELATIVE ERROR ELLIPSES DESIRED FORMAT I4 00000920 C 00000930 C 2.6 SEQUENCE NUMBERS OF STATIONS BETWEEN WHICH RELATIVE ERROR 00000940 C ELLIPSES TO BE COMPUTED FORMAT 2I5 00000950 C 00000960 C 00000970 C NOTE FOR DATA SET SEQUENCE NUMBERS DO NOT USE 5,6,7 00000980 C IF MORE THAN 18 WEIGHTED POINTS CONTINUE ON ANOTHER CARD 00000990 C WHEN ENTERING PX MATRIX THE MAIN DIAGONAL ELEMENT OF EACH ROW 00001000 C MUST START A NEW CARD -- EXAMPLE FOR 10*10 00001010 C 00001020 C 00001030 C XXXXXXXX 00001040 C XX 00001050 C XXXXXXXX 00001060 C X 00001070 C XXXXXXXX 12 CARDS IN TOTAL 00001080 C XXXXXXX 00001090 C XXXXXX 00001100 C XXXXX 00001110 C XXXX 00001120 C XXX 00001130 C XX 00001140 C X 00001150 C 00001160 C PROGRAM PX WILL PUNCH THE PX MATRIX IF IT COMES 00001170 C FROM DOPPLER SATELLITE DETERMINED COORDINATES 00001180 C 00001190 C CARD 2.5 MUST BE PRESENT EVEN IF NO ERROR ELLIPSES ARE 00001200 C REQUIRED (USE BLANK CARD) 00001210 C CARD 2.6 IS NOT REQUIRED IF THERE ARE NO ELLIPSES 00001220 C EACH ELLIPSE DESIRED REPRESENTS ONE CARD TYPE 2.6 00001230 C 00001240 C DIMENSIONING DATA 00001250 C 00001260 C THE FOLLOWING THREE DIMENSION STATEMENTS CHANGE WITH EVERY 00001270 C NETWORK 00001280 C 00001290 C 1 DIMENSION WTAZ(NK3),WTDR(NK3),WTDS(NK3),WVAZ(NK3),WVDR(NK3),WVDS(NK300001300 C ),MAZ(NK3,4),MDIR(NK3,4),MDIS(NK3,4),STOR(NK3,NK4),STORAZ(NK3,4), 00001310 C STORDR(NK3,4),STORDS(NK3,4),WEIGHT(NK3),CONST(NK3),ICOL(NK3), 00001320 C JCOL(NK4),W(NK4,NK3),AA(NK4,NK4),AL(NK4),UP(NK3,NK3),V(NK3) 00001330 C 00001340 C 2 DIMENSION B(NK2),XVEC(NK2),PHIV(NK2),BL(NK2),TA(NK5,NK2) 00001350 C 00001360 C 3 DIMENSION IVEC(NFP),INEW(NFP),XSAVE(NFPT),PX(NFPT,NFPT), 00001370 C ISTAFP(NFP),WTFP(NFPT,NFPT) 00001380 C 00001390 C WHERE NK2= NUMBER OF ELEMENTS IN SOLUTION VECTOR 00001400 C NK3= MAXIMUM NUMBER OF DIRECTION OR DISTANCE OR AZIMUTH TYPE 00001410 C OBSERVATIONS PLUS ONE (+1) AT ANY ONE STATION 00001420 C NK4= 2*NK3 00001430 C NK5= BANDWIDTH + BORDERWIDTH (NW+NB --SEE BELOW) 00001440 C NFP= NUMBER OF POINTS WITH A-PRIORI DATA 00001450 C NFPT= 2*NFP 00001460 C 00001470 C IN ADDITION THE FOLLOWING CARDS ARE REQUIRED TO BE PUNCHED 00001480 C EACH RUN 00001490 C 00001500 C NK1=XXXX 00001510 C NK2=XXXX 00001520 C NK3=XXXX 00001530 C NK4=XXXX 00001540 C NK5=XXXX 00001550 C NW=XXXX 00001560 C NB=XXXX 00001570 C NFP=XXXX 00001580 C NFPT=XXXX 00001590 C N=XXXX 00001600 C NBB=XXXX 00001610 C NSPACE=XXXX 00001620 C 00001630 C WHERE NK1= NUMBER OF NETWORK STATIONS 00001640 C NW= BANDWIDTH 00001650 C NB= BORDERWIDTH 00001660 C N= NK2-NB 00001670 C NBB= (NB**2+NB)/2 00001680 C NSPACE= N-NW+1 00001690 C 00001700 C ALL OF THE ABOVE CARDS CAN BE PUNCHED BY PROGRAM FORM3 BY SPECIFYING 00001710 C 00001720 C IPUN=1 ==> DIMENSION 1 ,NK1,NK2,NK3,NK4 00001730 C =2 ==> DIMENSION 2 ,NK5,NW,NB,N,NBB,NSPACE 00001740 C =3 ==> DIMENSION 3, NFP,NFPT 00001750 C =4 ==> BOTH 1 & 2 ABOVE 00001760 C =5 ==> BOTH 2 & 3 ABOVE 00001770 C =6 ==> BOTH 1 & 3 ABOVE 00001780 C =7 ==> 1 & 2 & 3 ABOVE 00001790 C 00001800 C 00001810 C ***** IMPORTANT NOTE ***** 00001820 C 00001830 C ALL DIMENSIONS ARE FIXED -- NO OVERDIMENSIONING IS ALLOWED 00001840 C 00001850 C 00001860 C 00001870 C THE FOLLOWING CONDITION CODES APPLY 00001880 C 00001890 C 1. STOP SOLUTION NOT CONVERGED - ONLY SOLUTION VECTOR COMPUTED 00001900 C 00001910 C 2. STOP 100 SOLUTION CONVERGED - SOLUTION VECTOR,INVERSE,RESIDUALS, 00001920 C VARIANCE FACTOR,ERROR ELLIPSES 00001930 C COMPUTED 00001940 C 00001950 C 3. STOP 200 SOLUTION DIVERGED - ONLY SOLUTION COMPUTED 00001960 C 00001970 C THE FOLLOWING CRITERIA IS USED IN CHECKING THE SOLUTION VECTOR 00001980 C 00001990 C ALL X(I) LESS THAN 0.001 ARC-SECONDS = STOP 100 00002000 C ANY X(I) GREATER THAN 0.001 BUT ALL X(I) LESS THAN 5.0 ARC-SECONDS 00002010 C STOP 00002020 C ANY X(I) GREATER THAN 5.0 ARC-SECONDS = STOP 200 00002030 C 00002040 C NOTE X(I) IN ABSOLUTE VALUE 00002050 C 00002060 C 00002070 C FOR FULL DETAILS SEE 00002080 C 00002090 C 1. PROGRAM PACKAGE FOR THE RIGOROUS COMPUTATION OF HORIZONTAL 00002100 C GEODETIC NETWORKS -- C.A.CHAMBERLAIN 00002110 C 00002120 C 2. LEAST SQUARES ADJUSTMENT OF HORIZONTAL CONTROL ON A MAPPING 00002130 C PLANE -- P.A.STEEVES 00002140 C 00002150 C 00002160 C 00002170 C 00002180 IMPLICIT REAL*8(A-H,O-Z) 00002190 DIMENSION DJOB(20),ITIME(17),IOUNIT(7) 00002200 DIMENSION WTAZ( 8),WTDR( 8),WTDS( 8),WVAZ( 8),WVDR( 8),WVDS( 8), 00002210 1MAZ( 8,4),MDIR( 8,4),MDIS( 8,4),STOR( 8,16),STORAZ( 8,4), 00002220 2 STORDR( 8,4),STORDS( 8,4), 00002230 3WEIGHT( 8),CONST( 8),ICOL( 8),JCOL(16),W(16, 8),AA(16,16),AL(16), 00002240 4 UP( 8, 8),V( 8) 00002250 DIMENSION B( 106),XVEC( 106),PHIV( 106),BL( 106),TA( 84, 106) 00002260 DIMENSION IVEC( 1),INEW( 1),XSAVE( 2),PX( 2, 2), 00002270 1 ISTAFP( 1),WTFP( 2, 2) 00002280 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00002290 CALL CPUTIM(ITIME( 1)) 00002300 CALL WTO('GOOD EVENING \') 00002310 NK1= 53 00002320 NK2= 106 00002330 NK3= 8 00002340 NK4= 16 00002350 NK5= 84 00002360 NW= 18 00002370 NB= 66 00002380 N= 40 00002390 NBB= 2211 00002400 NSPACE= 23 00002410 NFP= 1 00002420 NFPT= 2 00002430 C 00002440 C READ DATA CARDS 1 AND 2 00002450 C 00002460 READ(5,1000)DJOB 00002470 PRINT 1010,DJOB 00002480 READ 1020,IOUNIT 00002490 C 00002500 C READ IN THE SEQUENCE NUMBERS AND THEN THE UPPER TRIANGULAR PORTION OF00002510 C THE WEIGHTS TO BE ASSIGNED TO THE WEIGHTED STATIONS (ROW BY ROW) 00002520 C 00002530 READ(5,1030)(ISTAFP(I),I=1,NFP) 00002540 DO 10 I=1,NFPT 00002550 10 READ(5,1040)(WTFP(I,J),J=I,NFPT) 00002560 C 00002570 C FILL IN LOWER TRIANGULAR PORTION OF WEIGHT MATRIX 00002580 C 00002590 DO 20 I=1,NFPT 00002600 K=I+1 00002610 IF(K.GT.NFPT)GO TO 30 00002620 DO 20 J=K,NFPT 00002630 20 WTFP(J,I)=WTFP(I,J) 00002640 30 PRINT 1050,(ISTAFP(I),I=1,NFP) 00002650 PRINT 1060 00002660 DO 40 I=1,NFPT 00002670 40 PRINT 1070,(WTFP(I,J),J=I,NFPT) 00002680 C 00002690 C CALL U N B S A D 00002700 C 00002710 CALL UNBSAD(NK1,NK2,NK3,NK4,NK5,NFP,NFPT,WTAZ,WTDR,WTDS,WVAZ, 00002720 1WVDR,WVDS,MAZ,MDIR,MDIS,STOR,STORAZ,STORDR,STORDS,WEIGHT,CONST, 00002730 2ICOL,JCOL,W,AA,AL,UP,V,IVEC,INEW,XSAVE,PX,ISTAFP,WTFP,B,XVEC, 00002740 3PHIV,BL,TA) 00002750 PRINT 1080 00002760 DO 50 I=1,17 00002770 Z=DFLOAT(ITIME(I))/10000.0D0 00002780 50 PRINT 1090,I,Z 00002790 PRINT 1100 00002800 STOP 100 00002810 1000 FORMAT(20A4) 00002820 1010 FORMAT('1',10X,20A4///) 00002830 1020 FORMAT(7I2) 00002840 1030 FORMAT(6I5) 00002850 1040 FORMAT(8F10.0) 00002860 1050 FORMAT(30X,'SEQUENCE NUMBERS OF FIXED STATIONS'//10X,10I5//10X,'PX00002870 1 MATRIX TO BE ADDED') 00002880 1060 FORMAT('0',15X,'PX MATRIX TO BE ADDED'//) 00002890 1070 FORMAT(1X,10F10.1) 00002900 1080 FORMAT('-',10X,'CPU TIME IN SECONDS BETWEEN WTO''S') 00002910 1090 FORMAT('0',18X,I4,F15.4) 00002920 1100 FORMAT('0',15X,'COORDINATES HAVE CONVERGED') 00002930 END 00002940 SUBROUTINE UNBSAD(NK1,NK2,NK3,NK4,NK5,NFP,NFPT,WTAZ,WTDR,WTDS,WVAZ00002950 1,WVDR,WVDS,MAZ,MDIR,MDIS,STOR,STORAZ,STORDR,STORDS,WEIGHT,CONST, 00002960 2ICOL,JCOL,W,AA,AL,UP,V,IVEC,INEW,XSAVE,PX,ISTAFP,WTFP,B,XVEC, 00002970 3PHIV,BL,TA) 00002980 IMPLICIT REAL*8(A-H,O-Z) 00002990 REAL*4 WV,EV,ASEC,BSEC 00003000 DIMENSION ITIME(17),IOUNIT(7) 00003010 DIMENSION WTAZ(NK3),WTDR(NK3),WTDS(NK3),WVAZ(NK3), 00003020 1WVDR(NK3),WVDS(NK3),MAZ(NK3,4),MDIR(NK3,4),MDIS(NK3,4), 00003030 3STOR(NK3,NK4),STORAZ(NK3,4),STORDR(NK3,4),STORDS(NK3,4), 00003040 4WEIGHT(NK3),CONST(NK3),ICOL(NK3),JCOL(NK4),W(NK4,NK3),AA(NK4,NK4),00003050 5AL(NK4),UP(NK3,NK3),V(NK3),IVEC(NFP),INEW(NFP),XSAVE(NFPT), 00003060 6PX(NFPT,NFPT),ISTAFP(NFP),WTFP(NFPT,NFPT),B(NK2),XVEC(NK2), 00003070 7PHIV(NK1),BL(NK2),TA(NK5,NK2) 00003080 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00003090 C 00003100 C CONSTANTS 00003110 C 00003120 RAD=206264.8062470963D0 00003130 CFACT=2.449 00003140 C 00003150 C ZERO NORMAL EQUATION MATRIX (TA) 00003160 C 00003170 DO 10 I=1,NK2 00003180 DO 10 J=1,NK5 00003190 10 TA(J,I) = 0.0D0 00003200 DO 20 I=1,NK2 00003210 20 BL(I) = 0.0D0 00003220 C 00003230 C READ IN OBSERVATION EQUATIONS FROM SEQUENTIAL DATA FILE 00003240 C PERTAINING TO ONE OBSERVATION STATION AT A TIME 00003250 C 00003260 C THE OBSERVATION EQUATIONS ARE MIXED OF AZIMUTHS, DIRECTIONS AND 00003270 C DISTANCES SO FIRST MUST BE SORTED INTO THREE GROUPS, STORAZ, 00003280 C STORDR, STORDS 00003290 C 00003300 C THE NUMBERS I1,I2,I3,I4,I5 ARE THE COLUMN NUMBERS OF THE OBSERVATI00003310 C EQUATION COEFICIENTS T1,T2,T3,T4,T5 00003320 C EV = THE WEIGHT OF OBSERVATION EQUATION 00003330 C WV = THE MISCLOSURE OF THE OBSERVATION EQUATION 00003340 C NCODE IS AN INDICATOR OF THE TYPE OF OBSERVATION EQUATION 00003350 C 00003360 C NCODE: DIR = 30 00003370 C DIST = 40 00003380 C AZMTH = 50 00003390 C 00003400 C NEQ = THE NUMBER OF OBSERVATION EQUATIONS FROM THIS PARTICULAR 00003410 C STATION 00003420 C 00003430 IUNIT=IOUNIT(1) 00003440 REWIND IUNIT 00003450 READ(IUNIT) NOSETS 00003460 DO 110 IJK = 1,NOSETS 00003470 I = 0 00003480 J = 0 00003490 K = 0 00003500 MEQ = 1 00003510 READ(IUNIT)NEQ,T1,T2,T3,T4,T5,WV,EV,I1,I2,I3,I4,I5,NCODE 00003520 GO TO 40 00003530 30 MEQ = MEQ + 1 00003540 READ(IUNIT) T1,T2,T3,T4,T5,WV,EV,I1,I2,I3,I4,I5,NCODE 00003550 40 IF(NCODE.EQ.30) GO TO50 00003560 IF(NCODE.EQ.40) GO TO 60 00003570 I = I + 1 00003580 C 00003590 C AZIMUTHS 00003600 C 00003610 CALL READAZ(STORAZ,T1,T2,T3,T4,I1,I2,I3,I4,MAZ,WTAZ,WVAZ,EV,WV,I,N00003620 1K3) 00003630 GO TO 70 00003640 50 J = J + 1 00003650 C 00003660 C DIRECTIONS 00003670 C 00003680 CALL READDR(STORDR,T1,T2,T3,T4,I1,I2,I3,I4,MDIR,WTDR,WVDR,EV,WV,J,00003690 1NK3) 00003700 GO TO 70 00003710 60 K = K + 1 00003720 C 00003730 C DISTANCES 00003740 C 00003750 CALL READDS(STORDS,T1,T2,T3,T4,I1,I2,I3,I4,MDIS,WTDS,WVDS,EV,WV,K,00003760 1NK3) 00003770 70 CONTINUE 00003780 IF(MEQ.LT.NEQ) GO TO 30 00003790 C 00003800 C EACH SET OF SORTED OBSERVATION EQUATIONS IS NOW EXPANDED INTO A 00003810 C FOURTH BLOCK AND THE CONTRIBUTION TO THE NORMAL EQUATIONS FOUND 00003820 C AND ADDED TO THE NORMAL EQUATION MATRIX 00003830 C 00003840 IF(I.EQ.0) GO TO 80 00003850 C 00003860 C AZIMUTHS 00003870 C 00003880 II = 2*(I+1) 00003890 JN = 2 00003900 KN = 1 00003910 CALL SORT(I,WEIGHT,WTAZ,CONST,WVAZ,STOR,STORAZ,JN,KN,ICOL,MAZ,NK3,00003920 1NK4) 00003930 CALL DIRNOR(STOR,WEIGHT,CONST,AL,AA,I,II,UP,W,NK3,NK4,1) 00003940 CALL COLUMN(II,ICOL,JCOL,NK3,NK4) 00003950 CALL ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5) 00003960 CALL ADDVEC(BL,AL,II,JCOL,NK2,NK4) 00003970 80 CONTINUE 00003980 IF(J.EQ.0) GO TO 90 00003990 C 00004000 C DIRECTIONS 00004010 C 00004020 II = 2*(J+1) 00004030 JN = 2 00004040 KN = 1 00004050 CALL SORT(J,WEIGHT,WTDR,CONST,WVDR,STOR,STORDR,JN,KN,ICOL,MDIR,NK300004060 1,NK4) 00004070 CALL DIRNOR(STOR,WEIGHT,CONST,AL,AA,J,II,UP,W,NK3,NK4,0) 00004080 CALL COLUMN(II,ICOL,JCOL,NK3,NK4) 00004090 CALL ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5) 00004100 CALL ADDVEC(BL,AL,II,JCOL,NK2,NK4) 00004110 90 CONTINUE 00004120 IF(K.EQ.0) GO TO 100 00004130 C 00004140 C DISTANCES 00004150 C 00004160 II = 2*(K+1) 00004170 JN = 2 00004180 KN=1 00004190 CALL SORT(K,WEIGHT,WTDS,CONST,WVDS,STOR,STORDS,JN,KN,ICOL,MDIS,NK300004200 1,NK4) 00004210 CALL DISNOR(STOR,WEIGHT,CONST,AL,AA,K,II,NK3,NK4) 00004220 CALL COLUMN(II,ICOL,JCOL,NK3,NK4) 00004230 CALL ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5) 00004240 CALL ADDVEC(BL,AL,II,JCOL,NK2,NK4) 00004250 100 CONTINUE 00004260 110 CONTINUE 00004270 CALL CPUTIM(ITIME( 2)) 00004280 CALL WTO('THE NORMALS HAVE BEEN FORMED \') 00004290 C 00004300 C ADD PX TO NORMAL EQUATIONS 00004310 C 00004320 CALL PXADD(NK1,NK2,NK5,NFP,NFPT,INEW,IVEC,ISTAFP,WTFP,PX,TA) 00004330 C 00004340 C NORMAL EQUATION MATRIX IS FORMED 00004350 C CALL ARROW FOR SOLUTION OF THE NORMAL EQUATIONS 00004360 C 00004370 CALL CPUTIM(ITIME( 3)) 00004380 CALL WTO('ARROW HAS BEEN CALLED \') 00004390 CALL ARROW(NK2,NK5,NFP,NFPT,ISTOP,TA,B,XVEC,BL,XMAX,ISTAFP,XSAVE) 00004400 CALL CPUTIM(ITIME(13)) 00004410 CALL WTO('I HAVE RETURNED FROM ARROW \') 00004420 C 00004430 C CREATE DATA SET FOR SOLUTION VECTOR 00004440 C 00004450 IUNIT=IOUNIT(3) 00004460 REWIND IUNIT 00004470 DO 120 LM = 1,NK2 00004480 XVEC(LM)=-XVEC(LM) 00004490 120 WRITE(IUNIT)XVEC(LM) 00004500 PRINT 1000 00004510 IUNIT=IOUNIT(2) 00004520 JUNIT=IOUNIT(4) 00004530 REWIND IUNIT 00004540 K1 = 0 00004550 C 00004560 C COMPUTE ADJUSTED COORDINATES 00004570 C 00004580 DO 140 KIJ = 1,NK1 00004590 IFLAG=0 00004600 READ(IUNIT) IGEOD,PHI,ALONG 00004610 K1 = K1 + 1 00004620 PHI=PHI+XVEC(K1)/RAD 00004630 PHIV(KIJ)=PHI 00004640 CALL RADARC(PHI,I,J,ASEC) 00004650 K1=K1+1 00004660 ALONG=ALONG+XVEC(K1)/RAD 00004670 CALL RADARC(ALONG,K,L,BSEC) 00004680 K2 = K1/2 00004690 IF((K2/25*25).EQ.K2)PRINT 1000 00004700 DO 130 ICK=1,NFP 00004710 IF(ISTAFP(ICK).EQ.K2)IFLAG=1 00004720 130 CONTINUE 00004730 PRINT 1010,K2,IGEOD,I,J,ASEC,K,L,BSEC,XVEC(K1 - 1),XVEC(K1) 00004740 IF(IFLAG.EQ.1)PRINT 1020 00004750 C 00004760 C IF SOLUTION CONVERGED STORE ADJUSTED COORDINATES 00004770 C 00004780 IF(ISTOP.EQ.0)WRITE(JUNIT)IGEOD,PHI,ALONG 00004790 140 CONTINUE 00004800 CALL CPUTIM(ITIME(14)) 00004810 CALL WTO('THE ADJUSTED COORDINATES HAVE BEEN COMPUTED \') 00004820 C 00004830 C ERROR STOP FOR DIVERGENCE 00004840 C 00004850 IF(XMAX.GT.5.0)STOP 200 00004860 C 00004870 C IF ISTOP=1 SOLUTION NOT CONVERGED YET 00004880 C 00004890 IF(ISTOP.EQ.1)STOP 00004900 C 00004910 C READ IN OBSERVATION EQUATIONS FROM SEQUENTIAL DATA FILE 00004920 C PERTAINING TO ONE OBSERVATION STATION AT A TIME, THE FILE 00004930 C IS GONE THROUGH TWICE; THE FIRST TIME THROUGH THE AZIMUTH AND 00004940 C DIRECTION OBSERVATION ARE COLLECTED SO THAT THE RESIDUALS 00004950 C OF THESE MAY BE COMPUTED, THE SECOND TIME THROUGH IS FOR THE 00004960 C DISTANCE RESIDUAL COMPUTATION. A DATA SET IS CREATED 00004970 C CONTAINING THE RESIDUALS 00004980 C 00004990 N30=0 00005000 N40=0 00005010 N50=0 00005020 NEQT=0 00005030 NORU=0 00005040 PRINT 1030 00005050 IR = 0 00005060 VARFAC = 0.0D0 00005070 IUNIT=IOUNIT(1) 00005080 REWIND IUNIT 00005090 JUNIT=IOUNIT(5) 00005100 READ(IUNIT) NOSETS 00005110 DO 280 IJK = 1,NOSETS 00005120 MEQ = 1 00005130 I = 0 00005140 J = 0 00005150 READ(IUNIT)NEQ,T1,T2,T3,T4,T5,WV,EV,I1,I2,I3,I4,I5,NCODE 00005160 NEQT=NEQT+NEQ 00005170 IF(DABS(T5).GT.0.5)NORU=NORU+1 00005180 GO TO 160 00005190 150 MEQ = MEQ + 1 00005200 READ(IUNIT) T1,T2,T3,T4,T5,WV,EV,I1,I2,I3,I4,I5,NCODE 00005210 160 IF(NCODE.EQ.30) GO TO 170 00005220 IF(NCODE.EQ.40) GO TO 180 00005230 N50=N50+1 00005240 I = I + 1 00005250 CALL READAZ(STORAZ,T1,T2,T3,T4,I1,I2,I3,I4,MAZ,WTAZ,WVAZ,EV,WV,I,N00005260 1K3) 00005270 GO TO 180 00005280 170 J = J + 1 00005290 N30=N30+1 00005300 CALL READDR(STORDR,T1,T2,T3,T4,I1,I2,I3,I4,MDIR,WTDR,WVDR,EV,WV,J,00005310 1NK3) 00005320 180 CONTINUE 00005330 IF(MEQ.LT.NEQ) GO TO 150 00005340 C 00005350 C FORM RESIDUALS OF AZIMUTHS AND CONTRIBUTION TO VARIANCE FACTOR 00005360 C 00005370 IF(I.EQ.0) GO TO 210 00005380 DO 200 IN = 1,I 00005390 IR = IR + 1 00005400 V(IN) = 0.0D0 00005410 DO 190 INN = 1,4 00005420 190 V(IN) = STORAZ(IN,INN)*XVEC(MAZ(IN,INN)) + V(IN) 00005430 V(IN) = V(IN) + WVAZ(IN) 00005440 ISTA = MAZ(IN,2)/2 00005450 JSTA = MAZ(IN,4)/2 00005460 IF((IR/25*25).EQ.IR)PRINT 1030 00005470 PRINT 1040,IN,ISTA,JSTA,WVAZ(IN),V(IN) 00005480 NCODE=50 00005490 WRITE(JUNIT)IN,ISTA,JSTA,V(IN),WTAZ(IN),NCODE 00005500 200 VARFAC = VARFAC + V(IN)*WTAZ(IN)*V(IN) 00005510 C 00005520 C FORM RESIDUALS OF DIRECTIONS AND CONTRIBUTION TO VARIANCE FACTOR 00005530 C 00005540 210 IF(J.EQ.0) GO TO 280 00005550 SUM = 0.0D0 00005560 DO 220 IN = 1,J 00005570 220 SUM = SUM + WTDR(IN) 00005580 SUM = -1.0D0/SUM 00005590 DO 240 IN = 1,J 00005600 V(IN) = 0.0D0 00005610 DO 230 INN = 1,4 00005620 230 V(IN) = STORDR(IN,INN)*XVEC(MDIR(IN,INN)) + V(IN) 00005630 240 V(IN) = V(IN) + WVDR(IN) 00005640 ASUM = 0.0D0 00005650 DO 250 IN = 1,J 00005660 250 ASUM = ASUM + WTDR(IN)*V(IN) 00005670 DO 260 IN = 1,J 00005680 260 V(IN) = V(IN) + SUM*ASUM 00005960 DO 270 IN = 1,J 00005700 IR = IR + 1 00005710 ISTA = MDIR(IN,2)/2 00005720 JSTA = MDIR(IN,4)/2 00005730 IF((IR/25*25).EQ.IR)PRINT 1030 00005740 PRINT 1050,IN,ISTA,JSTA,WVDR(IN),V(IN) 00005750 NCODE=30 00005760 WRITE(JUNIT)IN,ISTA,JSTA,V(IN),WTDR(IN),NCODE 00005770 270 VARFAC = VARFAC + V(IN)*WTDR(IN)*V(IN) 00005780 280 CONTINUE 00005790 C 00005800 C RE-READ OBSERVATIONS LOOKING FOR DISTANCES 00005810 C 00005820 REWIND IUNIT 00005830 READ(IUNIT) NOSETS 00005840 DO 340 IJK = 1,NOSETS 00005850 MEQ = 1 00005860 K = 0 00005870 READ(IUNIT)NEQ,T1,T2,T3,T4,T5,WV,EV,I1,I2,I3,I4,I5,NCODE 00005880 GO TO 300 00005890 290 MEQ = MEQ + 1 00005900 READ(IUNIT) T1,T2,T3,T4,T5,WV,EV,I1,I2,I3,I4,I5,NCODE 00005910 300 IF(NCODE.NE.40)GO TO 310 00005920 N40=N40+1 00005930 K = K + 1 00005940 CALL READDS(STORDS,T1,T2,T3,T4,I1,I2,I3,I4,MDIS,WTDS,WVDS,EV,WV,K,00005950 1NK3) 00005960 310 CONTINUE 00005970 IF(MEQ.LT.NEQ) GO TO 290 00005980 C 00005990 C FORM RESIDUALS OF DISTANCES AND CONTRIBUTION TO VARIANCE FACTOR 00006000 C 00006010 IF(K.EQ.0) GO TO 340 00006020 DO 330 IN = 1,K 00006030 IR = IR + 1 00006040 V(IN) = 0.0D0 00006050 DO 320 INN = 1,4 00006060 320 V(IN) = V(IN) + STORDS(IN,INN)*XVEC(MDIS(IN,INN)) 00006070 V(IN) = V(IN) + WVDS(IN) 00006080 ISTA = MDIS(IN,2)/2 00006090 JSTA = MDIS(IN,4)/2 00006100 IF((IR/25*25).EQ.IR)PRINT 1030 00006110 PRINT 1060,IN,ISTA,JSTA,WVDS(IN),V(IN) 00006120 NCODE=40 00006130 WRITE(JUNIT)IN,ISTA,JSTA,V(IN),WTDS(IN),NCODE 00006140 330 VARFAC = VARFAC + V(IN)*WTDS(IN)*V(IN) 00006150 340 CONTINUE 00006160 CALL CPUTIM(ITIME(15)) 00006170 CALL WTO('THE RESIDUALS HAVE BEEN COMPUTED \') 00006180 NUN=NK2+NORU 00006190 NDF=NEQT-NUN+2*NFP 00006200 C 00006210 C COMPUTE XT*PX*X FOR POINTS WITH APRIORI INFORMATION 00006220 C 00006230 J=1 00006240 DO 360 K=1,NFPT 00006250 SUM=0.0D0 00006260 DO 350 I=1,NFPT 00006270 350 SUM=SUM+XSAVE(I)*WTFP(I,K) 00006280 BL(J)=SUM 00006290 360 J=J+1 00006300 XFAC=0.0D0 00006310 DO 370 I=1,NFPT 00006320 370 XFAC=XFAC+BL(I)*XSAVE(I) 00006330 PRINT 1070,VARFAC,NUN,XFAC,NFPT 00006340 C 00006350 C VARFAC=VT*P*V + XT*PX*X 00006360 C 00006370 VARFAC=VARFAC+XFAC 00006380 VARFAC=VARFAC/DFLOAT(NDF) 00006390 PRINT 1080,VARFAC 00006400 C 00006410 C CREATE 'SUMMARY' AND 'NORMALS' DATA SETS 00006420 C 00006430 IUNIT=IOUNIT(6) 00006440 JUNIT=IOUNIT(7) 00006450 WRITE(IUNIT)NK2,NK5,NW,NB 00006460 DO 380 I=1,NK2 00006470 380 WRITE(IUNIT)(TA(J,I),J=1,NK5) 00006480 WRITE(JUNIT)N30,N40,N50,NEQT,NK1,NUN,NDF,VARFAC 00006490 PRINT 1090 00006500 C 00006510 C STANDARD DEVIATIONS OF ADJUSTED COORDINATES 00006520 C 00006530 I=1 00006540 390 SP=DSQRT(TA(1,I)) 00006550 I=I+1 00006560 SL=DSQRT(TA(1,I)) 00006570 IF((J/25*25).EQ.J)PRINT 1090 00006580 PRINT 1100,J,SP,SL 00006590 I=I+1 00006600 J=J+1 00006610 IF(I.LT.NK2)GO TO 390 00006620 C 00006630 C ABSOLUTE AND RELATIVE ERROR ELLIPSES 00006640 C 00006650 CALL CPUTIM(ITIME(16)) 00006660 CALL WTO('ABSOLU HAS BEEN CALLED \') 00006670 CALL ABSOLU(TA,NK1,NK2,NK5,CFACT,PHIV) 00006680 CALL RELATI(TA,NK1,NK2,NK5,CFACT,PHIV) 00006690 PRINT 1110 00006700 WRITE(6,1120)N30,N40,N50,NEQT,NUN,NDF,VARFAC 00006710 CALL CPUTIM(ITIME(17)) 00006720 CALL WTO ('GOOD NIGHT \') 00006730 RETURN 00006740 1000 FORMAT('1', ///,35X,'A D J U S T E D C O O R D I N A T E S', 00006750 1 ' ',//,11X,'SEQ. #',8X,'GEOD. #',8X,'LATITUDE',14X,'LONGITU00006760 2DE',13X,'D. LAT.',3X,'D. LONG.', 00006770 3 ' ',/,36X,' DEG MIN SEC',7X,' DEG MIN SEC',10X,'SEC00006780 4',8X,'SEC',//) 00006790 1010 FORMAT(' ',/,I15,I15,I9,I5,F9.4,I9,I5,F9.4,F15.4,F11.4) 00006800 1020 FORMAT('+',1X,'WGHT. PT.') 00006810 1030 FORMAT('1',/////,46X,'R E S I D U A L S', 00006820 1 ' ',//,30X,'INST', 4X,'TARGET', 4X,'MISSCLOSURE', 7X,'RESID00006830 2UAL', 8X,'CODE') 00006840 1040 FORMAT(' ',/,I28,I6,I10,2F15.1,' AZIMUTH') 00006850 1050 FORMAT(' ',/,I28,I6,I10,2F15.1,' DIRECTION') 00006860 1060 FORMAT(' ',/,I28,I6,I10,2F15.3,' DISTANCE') 00006870 1070 FORMAT('1',30X,'VT*P*V =',F15.5,5X,'N-U =',I7//30X,'XT*PX*X =', 00006880 1F15.5,5X,'UX =',I7) 00006890 1080 FORMAT('-',/////47X,'VARIANCE FACTOR =',F15.5) 00006900 1090 FORMAT('1'///29X,'STANDARD DEVIATIONS OF ADJUSTED COORDINATES'/// 00006910 120X,'STATION NO.',10X,'SIGMA LATITUDE',10X,'SIGMA LONGTITUDE'/ 00006920 241X,'(ARC-SECONDS)',12X,'(ARC-SECONDS)') 00006930 1100 FORMAT('0',22X,I5,16X,F7.3,18X,F7.3) 00006940 1110 FORMAT('1',25X,'**** SUMMARY ****'///) 00006950 1120 FORMAT(15X,'NO. DIRECTION EQN.=',I9// 16X,'NO. DISTANCE EQN.=', 00006960 1I9//17X,'NO. AZIMUTH EQN.=',I9//19X,'TOTAL NO. EQN.=',I9//15X, 00006970 2'TOTAL NO. UNKNOWNS=',I9//15X,'DEGREES OF FREEDOM=',I9//18X, 00006980 3'VARIANCE FACTOR=',F9.3) 00006990 END 00007000 C 00007010 C ******************************************************************00007020 C * *00007030 C * S U B R O U T I N E R E A D D R *00007040 C * *00007050 C * THIS SUBROUTINE TAKES THE DESIGN MATRIX COEFFICIENTS THAT *00007060 C * CAME FROM THE STORAGE DEVICE AND STORES THEM IN A MATRIX *00007070 C * FOR DIRECTION TYPE EQUATIONS. *00007080 C * *00007090 C ******************************************************************00007100 C 00007110 SUBROUTINE READDR(S,T1,T2,T3,T4,I1,I2,I3,I4,MAZ,WT,WV,E,V,I,NK3) 00007120 IMPLICIT REAL*8(A-H,O-Z) 00007130 DIMENSION S(NK3,4),MAZ(NK3,4),WT(NK3),WV(NK3) 00007140 REAL*4 V,E 00007150 S(I,1) = T1 00007160 S(I,2) = T2 00007170 S(I,3) = T3 00007180 S(I,4) = T4 00007190 MAZ(I,1) = I1 00007200 MAZ(I,2) = I2 00007210 MAZ(I,3) = I3 00007220 MAZ(I,4) = I4 00007230 WV(I) = DBLE(V) 00007240 WT(I) = DBLE(E) 00007250 RETURN 00007260 END 00007270 C 00007280 C ******************************************************************00007290 C * *00007300 C * S U B R O U T I N E R E A D D S *00007310 C * *00007320 C * THIS SUBROUTINE TAKES THE DESIGN MATRIX COEFFICIENTS THAT *00007330 C * CAME FROM THE STORAGE DEVICE AND STORES THEM IN A MATRIX *00007340 C * FOR DISTANCE TYPE EQUATIONS. *00007350 C * *00007360 C ******************************************************************00007370 C 00007380 SUBROUTINE READDS(S,T1,T2,T3,T4,I1,I2,I3,I4,MAZ,WT,WV,E,V,I,NK3) 00007390 IMPLICIT REAL*8(A-H,O-Z) 00007400 DIMENSION S(NK3,4),MAZ(NK3,4),WT(NK3),WV(NK3) 00007410 REAL*4 V,E 00007420 S(I,1) = T1 00007430 S(I,2) = T2 00007440 S(I,3) = T3 00007450 S(I,4) = T4 00007460 MAZ(I,1) = I1 00007470 MAZ(I,2) = I2 00007480 MAZ(I,3) = I3 00007490 MAZ(I,4) = I4 00007500 WV(I) = DBLE(V) 00007510 WT(I) = DBLE(E) 00007520 RETURN 00007530 END 00007540 C 00007550 C ******************************************************************00007560 C * *00007570 C * S U B R O U T I N E R E A D A Z *00007580 C * *00007590 C * THIS SUBROUTINE TAKES THE DESIGN MATRIX COEFFICIENTS THAT *00007600 C * CAME FROM THE STORAGE DEVICE AND STORES THEM IN A MATRIX *00007610 C * FOR AZIMUTH TYPE EQUATIONS. *00007620 C * *00007630 C ******************************************************************00007640 C 00007650 SUBROUTINE READAZ(S,T1,T2,T3,T4,I1,I2,I3,I4,MAZ,WT,WV,E,V,I,NK3) 00007660 IMPLICIT REAL*8(A-H,O-Z) 00007670 DIMENSION S(NK3,4),MAZ(NK3,4),WT(NK3),WV(NK3) 00007680 REAL*4 V,E 00007690 S(I,1) = T1 00007700 S(I,2) = T2 00007710 S(I,3) = T3 00007720 S(I,4) = T4 00007730 MAZ(I,1) = I1 00007740 MAZ(I,2) = I2 00007750 MAZ(I,3) = I3 00007760 MAZ(I,4) = I4 00007770 WV(I) = DBLE(V) 00007780 WT(I) = DBLE(E) 00007790 RETURN 00007800 END 00007810 C 00007820 C ******************************************************************00007830 C * *00007840 C * S U B R O U T I N E S O R T *00007850 C * *00007860 C * THIS SUBROUTINE SORTS THE PARTIAL NORMAL EQUATION MATRIX *00007870 C * FOR THE ADDITION TO THE NORMAL EQUATION MATRIX. *00007880 C * *00007890 C ******************************************************************00007900 C 00007910 SUBROUTINE SORT(I,WEIGHT,WT,CONST,WV,SS,S,JN,KN,ICOL,MAZ,NK3,NK4) 00007920 IMPLICIT REAL*8(A-H,O-Z) 00007930 DIMENSION WEIGHT(NK3),CONST(NK3),WT(NK3),WV(NK3),SS(NK3,NK4) 00007940 DIMENSION S(NK3,4),MAZ(NK3,4),ICOL(NK3) 00007950 DO 10 IL = 1,NK4 00007960 DO 10 IN = 1,NK3 00007970 10 SS(IN,IL) = 0.0D0 00007980 ICOL(1) = MAZ(1,1) 00007990 DO 20 IN = 1,I 00008000 WEIGHT(IN) = WT(IN) 00008010 CONST(IN) = WV(IN) 00008020 SS(IN,1) = S(IN,1) 00008030 SS(IN,2) = S(IN,2) 00008040 JN = JN + 1 00008050 SS(IN,JN) = S(IN,3) 00008060 JN = JN + 1 00008070 SS(IN,JN) = S(IN,4) 00008080 KN = KN + 1 00008090 20 ICOL(KN) = MAZ(IN,3) 00008100 RETURN 00008110 END 00008120 C 00008130 C ******************************************************************00008140 C * *00008150 C * S U B R O U T I N E D I R N O R *00008160 C * *00008170 C * THIS SUBROUTINE ACCEPTS THE COMPACTED DESIGN MATRIX FROM *00008180 C * THE PARTICULAR OBSERVATION STATION RESULTING FROM EITHER *00008190 C * DIRECTIONS OR AZIMUTHS. IF THE OBSERVATION EQUATIONS ARE *00008200 C * FROM DIRECTIONS, THE ORIENTATION UNKNOWNS ARE ELIMINATED. *00008210 C * THE PARTIAL NORMAL EQUATION MATRIX AND CONSTANT VECTOR *00008220 C * ARE FORMED. *00008230 C * *00008240 C ******************************************************************00008250 C 00008260 SUBROUTINE DIRNOR(STOR,WEIGHT,CONST,AL,AA,NODIR,II,UP,W,NK3,NK4,IA00008270 1Z) 00008280 IMPLICIT REAL*8(A-H,O-Z) 00008290 DIMENSION STOR(NK3,NK4),WEIGHT(NK3),CONST(NK3),AL(NK4) 00008300 DIMENSION AA(NK4,NK4),UP(NK3,NK3),W(NK4,NK3) 00008310 SUM=0.0D0 00008320 C 00008330 C MULTIPLY ATUP BY CONSTANT VECTOR 'W' 00008340 C 00008350 DO 10 I=1, NODIR 00008360 DO 10 J=1, NODIR 00008370 10 UP(J,I) = 0.0D0 00008380 IF(IAZ.EQ.1) GO TO 40 00008390 DO 20 I=1,NODIR 00008400 20 SUM=SUM+WEIGHT(I) 00008410 DO 30 I=1,NODIR 00008420 DO 30 J=1,NODIR 00008430 UP(I,J)=-WEIGHT(I)*WEIGHT(J)/SUM 00008440 30 UP(J,I)=UP(I,J) 00008450 40 DO 50 I=1,NODIR 00008460 50 UP(I,I)=UP(I,I)+WEIGHT(I) 00008470 C 00008480 C MULTIPLY A TRANSPOSED BY MATRIX UP 00008490 C 00008500 SUM=0 00008510 DO 70 I=1,II 00008520 DO 70 J=1,NODIR 00008530 DO 60 K=1,NODIR 00008540 60 SUM=SUM+STOR(K,I)*UP(K,J) 00008550 W(I,J)=SUM 00008560 70 SUM=0.0D0 00008570 C 00008580 C MULTIPLY ATUP BY A MATRIX (PARTIAL) 00008590 SUM=0.0D0 00008600 DO 90 I=1,II 00008610 DO 90 J=1,II 00008620 DO 80 K=1,NODIR 00008630 80 SUM=SUM+W(I,K)*STOR(K,J) 00008640 AA(I,J)=SUM 00008650 90 SUM=0.0D0 00008660 C 00008670 C MULTIPLY ATUP BY CONSTANT VECTOR 'W' 00008680 C 00008690 SUM=0.0D0 00008700 DO 110 I=1,II 00008710 DO 100 K=1,NODIR 00008720 100 SUM=SUM+W(I,K)*CONST(K) 00008730 AL(I)=SUM 00008740 110 SUM=0.0D0 00008750 RETURN 00008760 END 00008770 C 00008780 C ******************************************************************00008790 C * *00008800 C * S U B R O U T I N E D I S N O R *00008810 C * *00008820 C * THIS SUBROUTINE ACCEPTS THE COMPACTED DESIGN MATRIX FROM *00008830 C * THE PARTICULAR OBSERVATION STATION, RESULTING FROM *00008840 C * DISTANCES. THE PARTIAL NORMAL EQUATION MATRIX AND CONSTANT *00008850 C * VECTOR ARE FORMED. *00008860 C * *00008870 C ******************************************************************00008880 C 00008890 SUBROUTINE DISNOR(STOR,WEIGHT,CONST,AL,AA,NODIST,II,NK3,NK4) 00008900 IMPLICIT REAL*8(A-H,O-Z) 00008910 DIMENSION STOR(NK3,NK4),WEIGHT(NK3),CONST(NK3),AL(NK4),AA(NK4,NK4)00008920 C 00008930 C MULTIPLY A TRANSPOSED BY P BY A MATRIX PARTIAL 00008940 C 00008950 SUM=0.0D0 00008960 DO 20 I=1,II 00008970 DO 20 J=1,II 00008980 DO 10 K=1,NODIST 00008990 10 SUM=SUM+STOR(K,I)*WEIGHT(K)*STOR(K,J) 00009000 AA(I,J)=SUM 00009010 20 SUM=0.0D0 00009020 C 00009030 C MULTIPLY A TRANSPOSED BY P BY CONSTANT VECTOR 'W' 00009040 C 00009050 SUM = 0.0D0 00009060 DO 40 I = 1,II 00009070 DO 30 K = 1,NODIST 00009080 30 SUM = SUM + STOR(K,I)*WEIGHT(K)*CONST(K) 00009090 AL(I) = SUM 00009100 40 SUM = 0.0D0 00009110 RETURN 00009120 END 00009130 C 00009140 C ******************************************************************00009150 C * *00009160 C * S U B R O U T I N E A D D N O R *00009170 C * *00009180 C * THIS SUBROUTINE ADDS THE PARTIAL NORMAL EQUATION MATRIX *00009190 C * TO THE NORMAL EQUATION MATRIX. THE NORMAL EQUATION MATRIX *00009200 C * IS COMPACTED INTO THE ARROWHEAD PATTERN. THIS IS *00009210 C * ACCOMPLISHED WITH THE USE OF NSPACE, NSPACE EQUALS THE *00009220 C * NUMBER OF ZEROS BETWEEN THE EDGE OF THE BAND AND THE EDGE *00009230 C * OF THE BORDER IN THE FIRST ROW OF THE MATRIX. NSPACE IS *00009240 C * SUBTRACTED FROM THE COLUMN INDICIE FOR THE FIRST ROW, FOR *00009250 C * THE SECOND ROW WE MUST ADD ONE, THIRD ROW ADD TWO, ETC. *00009260 C * *00009270 C ******************************************************************00009280 C 00009290 SUBROUTINE ADDNOR(TA,AA,II,JCOL,NK2,NK4,NK5) 00009300 IMPLICIT REAL*8(A-H,O-Z) 00009310 DIMENSION TA(NK5,NK2),AA(NK4,NK4),JCOL(NK4),ITIME(17),IOUNIT(7) 00009320 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00009330 C 00009340 C ADD THE SUB NORMALS TO THE NORMALS 00009350 C 00009360 ICON=NW-1 00009370 DO 20 I=1,II 00009380 IJC=JCOL(I) 00009390 NWIJC=ICON+IJC 00009400 DO 20 K=1,II 00009410 KJC=JCOL(K) 00009420 IF(KJC.LT.IJC) GO TO 20 00009430 IF(KJC.LE.NWIJC)GO TO 10 00009440 NSKIP=NSPACE-IJC 00009450 IF(NSKIP.LT.0) NSKIP=0 00009460 KJC=KJC-NSKIP 00009470 10 IKJC=KJC-IJC+1 00009480 TA(IKJC,IJC)=TA(IKJC,IJC)+AA(I,K) 00009490 20 CONTINUE 00009500 RETURN 00009510 END 00009520 C 00009530 C ******************************************************************00009540 C * *00009550 C * S U B R O U T I N E C O L U M N *00009560 C * *00009570 C * THIS SUBROUTINE COMPUTES THE COLUMN INDICIES FOR THE *00009580 C * ADDITION OF THE PARTIAL NORMAL EQUATION MATRIX TO THE *00009590 C * NORMAL EQUATION MATRIX. *00009600 C * *00009610 C ******************************************************************00009620 C 00009630 SUBROUTINE COLUMN(II,ICOL,JCOL,NK3,NK4) 00009640 DIMENSION ICOL(NK3),JCOL(NK4) 00009650 C 00009660 C COMPUTE COLUMN INDICES FOR THE ADDITION OF THE SUB NORMALS 00009670 C 00009680 IT=II/2 00009690 K=0 00009700 IF(IT.EQ.0) GO TO 20 00009710 DO 10 I=1,IT 00009720 K=K+1 00009730 JCOL(K)=ICOL(I) 00009740 K=K+1 00009750 10 JCOL(K)=ICOL(I)+1 00009760 RETURN 00009770 20 JCOL(1)=ICOL(1) 00009780 JCOL(2)=ICOL(1)+1 00009790 RETURN 00009800 END 00009810 C 00009820 C ******************************************************************00009830 C * *00009840 C * S U B R O U T I N E A D D V E C *00009850 C * *00009860 C * THIS SUBROUTINE ADDS THE PARTIAL CONSTANT VECTOR ATPW TO *00009870 C * THE CONSTANT VECTOR ATPW. *00009880 C * *00009890 C ******************************************************************00009900 C 00009910 SUBROUTINE ADDVEC(BL,AL,II,JCOL,NK2,NK4) 00009920 IMPLICIT REAL*8(A-H,O-Z) 00009930 DIMENSION BL(NK2),AL(NK4),JCOL(NK4) 00009940 C 00009950 C ADD THE SUB VECTOR BL TO THE VECTOR 'L' 00009960 C 00009970 DO 10 I=1,II 00009980 IJC=JCOL(I) 00009990 10 BL(IJC)=BL(IJC)+AL(I) 00010000 RETURN 00010010 END 00010020 C 00010030 C SUBROUTINE 'RADARC' CONVERTS RADIANS TO DEGREES MINUTES AND 00010040 C SECONDS. FOR NEGATIVE ANGLES ONLY THE LEFTMOST NONZERO VALUE IS 00010050 C NEGATIVE (EGS. -50,15,30.5 ; 0,-35,30.0 ; 0,0,-50.5) 00010060 C 00010070 C NOTE: THE 0.0005 VALUE IS TO GUARD AGAINST ROUNDOFF 00010080 C 00010090 C INPUT: A = RADIAN VALUE OF ANGLE (REAL*8) 00010100 C 00010110 C OUTPUT: I = DEGREES (INTEGER) 00010120 C J = MINUTES (INTEGER) 00010130 C S = SECONDS (REAL*4) 00010140 C 00010150 SUBROUTINE RADARC(A,I,J,S) 00010160 DOUBLE PRECISION A,SEC,AD,AJ,RHO, SIGN 00010170 DATA RHO/206264.8062470963D0/ 00010180 C 00010190 C CHECK SIGN OF 'A' -- SET SIGN=-1 IF NEGATIVE AND CONVERT 'A' TO 00010200 C POSITIVE VALUE 00010210 C 00010220 SIGN=1.0D0 00010230 IF(A.LT.0.0)SIGN=-1.0D0 00010240 IF(SIGN.LT.0.0)A=-A 00010250 C 00010260 C CONVERT 'A' TO ARCSECONDS 00010270 C 00010280 SEC=A*RHO+0.0005D0 00010290 C 00010300 C FIND INTEGER DEGREES 00010310 C 00010320 I=SEC/3600.0D0 00010330 AD=I 00010340 C 00010350 C FIND INTEGER MINUTES 00010360 C 00010370 J=SEC/60.0D0-AD*60.0D0 00010380 AJ=J 00010390 C 00010400 C FIND REAL*4 SECONDS 00010410 C 00010420 S=SEC-AD*3600.0D0-AJ*60.0D0-0.0005D0 00010430 C 00010440 C SET LEFTMOST VALUE NEGATIVE IF SIGN=-1 00010450 C 00010460 IF(I.NE.0)GO TO 20 00010470 IF(J.EQ.0)GO TO 10 00010480 J=J*SIGN 00010490 GO TO 30 00010500 10 S=S*SIGN 00010510 GO TO 30 00010520 20 I=I*SIGN 00010530 C 00010540 C CONVERT 'A' BACK TO NEGATIVE IF SIGN=-1 00010550 C 00010560 30 IF(SIGN.LT.0.0)A=-A 00010570 RETURN 00010580 END 00010590 SUBROUTINE ARROW(NNB,NWB,NFP,NFPT,ISTOP,A,B,XVEC,W,XMAX,ISTAFP, 00010600 1XSAVE) 00010610 C 00010620 C ******************************************************************00010630 C * *00010640 C * S U B R O U T I N E A R R O W *00010650 C * *00010660 C * THIS SUBROUTINE SOLVES A BANDED AND BORDERED SYMMETRIC *00010670 C * POSITIVE DEFINITE SET OF NORMAL EQUATIONS. THE BANDED AND *00010680 C * BORDERED PORTION ONLY OF THE VARIANCE CO-VARIANCE MATRIX *00010690 C * IS SOLVED FOR,TOGETHER WITH THE X-VECTOR. THE METHOD USED *00010700 C * IS THAT OF CHOLESKI. *00010710 C * *00010720 C ******************************************************************00010730 C 00010740 REAL*8 A(NWB,NNB),B(NNB),XVEC(NNB),W(NNB),CONV,DABS,XMAX, 00010750 1XSAVE(NFPT) 00010760 DIMENSION ISTAFP(NFP),ITIME(17),IOUNIT(7) 00010770 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00010780 C 00010790 C THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES 00010800 C 1) TRIANG 00010810 C 2) FORWAR AND BACK 00010820 C 3) ARRINV 00010830 C 00010840 C FIND THE CHOLESKI SQUARE ROOT 00010850 C 00010860 CALL TRIANG(A,NNB,NWB) 00010870 CALL CPUTIM(ITIME( 8)) 00010880 CALL WTO('I HAVE RETURNED FROM TRIANG \') 00010890 C 00010900 C CLEAR X AND B VECTORS 00010910 C 00010920 DO 10 I = 1,NNB 00010930 XVEC(I) = 0.0D0 00010940 10 B(I) = 0.0D0 00010950 C 00010960 C COMPUTE FORWARD SOLUTION OF X-VECTOR 00010970 C 00010980 CALL FORWAR(A,W,B,NNB,NWB) 00010990 CALL CPUTIM(ITIME( 9)) 00011000 CALL WTO('I HAVE RETURNED FROM FORWAR \') 00011010 C 00011020 C COMPUTE BACK SOLUTION FOR X-VECTOR 00011030 C 00011040 CALL BACK(A,B,XVEC,NNB,NWB) 00011050 CALL CPUTIM(ITIME(10)) 00011060 CALL WTO('I HAVE RETURNED FROM BACK \') 00011070 C 00011080 C CHECK SOLUTION VECTOR FOR CONVERGENCE. CHECK ONLY AT POINTS 00011090 C WITHOUT A-PRIORI KNOWLEDGE (CORRESPONDING ELEMENTS OF SOLUTION 00011100 C VECTOR SAVED IN XSAVE) 00011110 C CRITERIA FOR CONVERGENCE = 0.001 ARC-SECONDS 00011120 C CRITERIA FOR DIVERGENCE = 5.0 ARC-SECONDS 00011130 C FIND INVERSE ONLY IF ISTOP=0 (IE. CONVERGENCE) 00011140 C 00011150 ISTOP=0 00011160 XMAX=0.0D0 00011170 CONV=0.001D0 00011180 DO 20 LM=1,NFP 00011190 LM2=2*LM 00011200 IR1=2*ISTAFP(LM) 00011210 IR2=IR1-1 00011220 XSAVE(LM2)=XVEC(IR1) 00011230 XSAVE(LM2-1)=XVEC(IR2) 00011240 XVEC(IR2)=0.0D0 00011250 20 XVEC(IR1)=0.0D0 00011260 DO 30 I=1,NNB 00011270 IF(DABS(XVEC(I)).GE.CONV)GO TO 40 00011280 30 CONTINUE 00011290 GO TO 60 00011300 40 ISTOP=1 00011310 CALL WTO('SOLUTION HAS NOT CONVERGED--ANOTHER ITERATION IS REQUIRE00011320 1D \') 00011330 DO 50 I=1,NNB 00011340 IF(DABS(XVEC(I)).GT.XMAX)XMAX=DABS(XVEC(I)) 00011350 50 CONTINUE 00011360 RETURN 00011370 C 00011380 C ELEMENT (N,N) OF VARIANCE CO-VARIANCE IS (RESTORE XVEC BEFORE 00011390 C COMPUTATION) 00011400 C 00011410 60 DO 70 LM=1,NFP 00011420 LM2=2*LM 00011430 IR1=2*ISTAFP(LM) 00011440 IR2=IR1-1 00011450 XVEC(IR2)=XSAVE(LM2-1) 00011460 70 XVEC(IR1)=XSAVE(LM2) 00011470 A(1,NNB)=1.0D0/A(1,NNB)/A(1,NNB) 00011480 MW = NWB - 1 00011490 I = NNB 00011500 J1 = 1 00011510 I2 = NWB - 1 00011520 C 00011530 C COMPUTE LOWER BLOCK OF VARIANCE CO-VARIANCE MATRIX 00011540 C 00011550 CALL CPUTIM(ITIME(11)) 00011560 CALL WTO('ARRINV HAS BEEN CALLED THE FIRST TIME \') 00011570 CALL ARRINV(A,B,NNB,NWB,MW,I,J1,I2) 00011580 J1 = NWB 00011590 MW = NNB - 1 00011600 C 00011610 C COMPUTE REMAINING PORTION OF THE VARIANCE CO-VARIANCE PERTAINING 00011620 C TO THE BANDED AND BORDERED PORTION OF THE NORMAL EQUATION MATRIX 00011630 C 00011640 CALL CPUTIM(ITIME(12)) 00011650 CALL WTO('ARRINV HAS BEEN CALLED THE SECOND TIME \') 00011660 CALL ARRINV(A,B,NNB,NWB,MW,I,J1,I2) 00011670 RETURN 00011680 END 00011690 SUBROUTINE TRIANG(A,NNB,NWB) 00011700 C 00011710 C ******************************************************************00011720 C * *00011730 C * S U B R O U T I N E T R I A N G *00011740 C * *00011750 C * THIS SUBROUTINE TAKES THE BANDED AND BORDERED NORMAL EQUATION*00011760 C * MATRIX AND PERFORMS THE CHOLESKI SQUARE ROOT (THAT IS *00011770 C * DECOMPOSES THE MATRIX INTO TWO TRIANGULAR MATRICES L AND LT; *00011780 C * WHERE LT IS THE TRANSPOSE OF L AND LT*L = N) *00011790 C * *00011800 C ******************************************************************00011810 C 00011820 C 00011830 C FIND SQRT OF FIRST DIAGONAL COEFFICIENT 00011840 C 00011850 REAL*8 A(NWB,NNB), DSQRT 00011860 DIMENSION ITIME(17),IOUNIT(7) 00011870 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00011880 A(1,1) = DSQRT(A(1,1)) 00011890 KW = NW + NB 00011900 II = 0 00011910 C 00011920 C DIVIDE REMAINING COEFFICIENTS IN FIRST ROW BY A(1,1) 00011930 C 00011940 DO 10 J = 2,NWB 00011950 10 A(J,1) = A(J,1)/A(1,1) 00011960 MW = N - NW + 1 00011970 C 00011980 C FIND THE CHOLESKI SQUARE ROOT FOR THE REMAINDER OF THE BAND AND 00011990 C BORDER DOWN TO WHERE THE BAND BEGINS TO DECREASE IN SIZE. 00012000 C 00012010 CALL CPUTIM(ITIME( 4)) 00012020 CALL WTO('THE FIRST CALL TO TRI1 IS BEING MADE \') 00012030 DO 20 I = 2,MW 00012040 I1 = I - 1 00012050 IF(I.GT.NW) I1 = NW - 1 00012060 CALL TRI1(A,I,NWB,I1,II,NNB) 00012070 20 CONTINUE 00012080 CALL CPUTIM(ITIME( 5)) 00012090 CALL WTO('THE FIRST LOOP THROUGH TRI1 IS FINISHED \') 00012100 JW = NW + 1 00012110 L1 = NW - 1 00012120 C 00012130 C FIND THE CHOLESKI SQUARE ROOT FOR THE REMAINDER OF THE BAND AND 00012140 C BORDER DOWN TO WHERE THE BAND VANISHES. 00012150 C 00012160 DO 30 II = 1,L1 00012170 JW = JW - 1 00012180 I = II + MW 00012190 KW = KW - 1 00012200 CALL TRI1(A,I,NWB,L1,II,NNB) 00012210 30 CONTINUE 00012220 CALL CPUTIM(ITIME( 6)) 00012230 CALL WTO('THE SECOND LOOP TO TRI1 IS FINISHED \') 00012240 MW = N + 1 00012250 ND = NB + 1 00012260 C 00012270 C UPDATE N22 TO EQUAL N22 - L12'*L12 00012280 C 00012290 KNW=0 00012300 DO 40 K=1,N 00012310 IF(K.GT.N-NW+1)KNW=KNW+1 00012320 DO 40 I=1,NB 00012330 IN=I+N 00012340 INW=I+NW 00012350 INWKNW=INW-KNW 00012360 DO 40 J = I,NB 00012370 JNW = J + NW 00012380 JI = J - I + 1 00012390 40 A(JI,IN) = A(JI,IN) - A(INWKNW ,K)*A(JNW - KNW,K) 00012400 CALL CPUTIM(ITIME( 7)) 00012410 CALL WTO('N22 IS UPDATED \') 00012420 C 00012430 C FIND CHOLESKI SQUARE ROOT FOR N22 00012440 C 00012450 DO 70 I = MW,NNB 00012460 ND = ND - 1 00012470 DO 50 J = 1,ND 00012480 IF(J.EQ.1) A(J,I) = DSQRT(A(J,I)) 00012490 IF(J.NE.1) A(J,I) = A(J,I)/A(1,I) 00012500 50 CONTINUE 00012510 IF(I.EQ.NNB) GO TO 70 00012520 DO 60 K = 2,ND 00012530 DO 60 J = K,ND 00012540 60 A(J-K+1,I+K-1) = A(J-K+1,I+K-1) - A(K,I)*A(J,I) 00012550 70 CONTINUE 00012560 RETURN 00012570 END 00012580 SUBROUTINE TRI1(A,I,NWB,I1,II,NNB) 00012590 C 00012600 C ******************************************************************00012610 C * *00012620 C * S U B R O U T I N E T R I 1 *00012630 C * *00012640 C * THIS SUBROUTINE IS CALLED BY TRIANG AND IS USED TO FIND *00012650 C * THE CHOLESKI SQUARE ROOT. THE SUBROUTINE DOES A DIFFERENT *00012660 C * JOB EACH TIME IT IS CALLED SO CHECK THE COMMENTS IN TRIANG *00012670 C * FOR THE BASIC IDEA BEHIND IT. *00012680 C * *00012690 C ******************************************************************00012700 C 00012710 DIMENSION ITIME(17),IOUNIT(7) 00012720 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00012730 REAL*8 A(NWB,NNB),DSQRT 00012740 IF(II.LE.1) J1 = NW - 1 00012750 IF(II.GT.1) J1 = J1 - 1 00012760 DO 10 K = 1,I1 00012770 IK = I - K 00012780 K2 = K + 1 00012790 10 A(1,I) = A(1,I) - A(K2,IK)*A(K2,IK) 00012800 A(1,I) = DSQRT(A(1,I)) 00012810 IF(II.EQ.I1) GO TO 40 00012820 IF(II.EQ.0) A(NW,I) = A(NW,I)/A(1,I) 00012830 IF(II.EQ.0) K1 = 1 00012840 IF(II.NE.0) K1 = II 00012850 I2 = I1 00012860 IF(I.GE.NW) I2 = I2 - 1 00012870 I3 = 0 00012880 IF(I.LT.(NW-1)) I3 = NW - I - 1 00012890 DO 20 L = 1,I2 00012900 I3 = I3 + 1 00012910 IF((II.GT.1).AND.(I3.GT.(NW-II-1))) I3 = I3 - 1 00012920 DO 20 J = 1,I3 00012930 J2 = J + 1 00012940 I4 = I2 + 2 - L 00012950 IL=I+L-I2-1 00012960 20 A(J2,I) = A(J2,I) - A(I4,IL)*A(I4+J,IL) 00012970 DO 30 J = 1,I3 00012980 J2 = J + 1 00012990 30 A(J2,I) = A(J2,I)/A(1,I) 00013000 40 CONTINUE 00013010 L3 = 0 00013020 DO 50 L = 1,I1 00013030 IF(L.LE.II) L3 = L 00013040 DO 50 J = 1,NB 00013050 NWJ = NW - II + J 00013060 IL = I - L 00013070 50 A(NWJ,I) = A(NWJ,I) - A(NWJ+L3,IL)*A(L+1,IL) 00013080 DO 60 J = 1,NB 00013090 NWJ = NW - II + J 00013100 60 A(NWJ,I) = A(NWJ,I)/A(1,I) 00013110 RETURN 00013120 END 00013130 SUBROUTINE FORWAR(A,W,B,NNB,NWB) 00013140 C 00013150 C ******************************************************************00013160 C * *00013170 C * S U B R O U T I N E F O R W A R *00013180 C * *00013190 C * THIS SUBROUTINE COMPUTES THE FORWARD SOLUTION OF THE *00013200 C * X-VECTOR FOR THE BANDED AND BORDERED NORMAL EQUATIONS, *00013210 C * METHOD OF CHOLESKI. *00013220 C * *00013230 C ******************************************************************00013240 C 00013250 REAL*8 A(NWB,NNB),W(NNB),B(NNB) 00013260 DIMENSION ITIME(17),IOUNIT(7) 00013270 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00013280 B(1) = W(1)/A(1,1) 00013290 DO 10 K = 1,NB 00013300 IK = K + N 00013310 10 B(IK) = B(IK) + B(1)*A(NW+K,1) 00013320 KW = NW 00013330 MW = N - NW + 1 00013340 DO 40 I = 2,N 00013350 IF(I.GT.MW) KW = KW - 1 00013360 I1 = I - 1 00013370 IF(I.GT.NW) I1 = NW - 1 00013380 DO 20 K = 1,I1 00013390 IK = I - K 00013400 20 B(I) = B(I) + B(IK)*A(K+1,IK) 00013410 B(I) = (W(I) - B(I))/A(1,I) 00013420 DO 30 K = 1,NB 00013430 IK = K + N 00013440 30 B(IK) = B(IK) + B(I)*A(KW+K,I) 00013450 40 CONTINUE 00013460 NBB = NB 00013470 MW = N + 1 00013480 DO 60 I = MW,NNB 00013490 NBB =NBB - 1 00013500 B(I) = (W(I) - B(I))/A(1,I) 00013510 IF(NBB.EQ.0) GO TO 60 00013520 DO 50 K = 1,NBB 00013530 IK = I + K 00013540 50 B(IK) = B(IK) + B(I)*A(K+1,I) 00013550 60 CONTINUE 00013560 RETURN 00013570 END 00013580 SUBROUTINE BACK(A,B,XVEC,NNB,NWB) 00013590 C 00013600 C ******************************************************************00013610 C * *00013620 C * S U B R O U T I N E B A C K *00013630 C * *00013640 C * THIS SUBROUTINE COMPUTES THE BACK SOLUTION OF THE X-VECTOR *00013650 C * FOR THE BANDED AND BORDERED NORMAL EQUATIONS, METHOD OF *00013660 C * CHOLESKI *00013670 C * *00013680 C ******************************************************************00013690 C 00013700 REAL*8 A(NWB,NNB),B(NNB),XVEC(NNB) 00013710 DIMENSION ITIME(17),IOUNIT(7) 00013720 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00013730 XVEC(NNB) = B(NNB)/A(1,NNB) 00013740 DO 20 I = 2,NB 00013750 I1 = I -1 00013760 I2 = NNB - I + 1 00013770 DO 10 K = 1,I1 00013780 10 XVEC(I2) = XVEC(I2) + XVEC(I2+K)*A(K+1,I2) 00013790 20 XVEC(I2) = (B(I2)-XVEC(I2))/A(1,I2) 00013800 DO 40 I = 1,NW 00013810 I2 = N - I + 1 00013820 I1 = I1 + 1 00013830 DO 30 K = 1,I1 00013840 30 XVEC(I2) = XVEC(I2) + XVEC(I2+K)*A(K+1,I2) 00013850 40 XVEC(I2) = (B(I2) - XVEC(I2))/A(1,I2) 00013860 MW = N - NW 00013870 I1 = NW - 1 00013880 DO 70 I = 1,MW 00013890 I2 = I2 - 1 00013900 DO 50 K = 1,NB 00013910 50 XVEC(I2) = XVEC(I2) + XVEC(NNB+1-K)*A(NWB+1-K,I2) 00013920 DO 60 K = 1,I1 00013930 60 XVEC(I2) = XVEC(I2) + XVEC(I2+K)*A(K+1,I2) 00013940 70 XVEC(I2) = (B(I2) - XVEC(I2))/A(1,I2) 00013950 RETURN 00013960 END 00013970 SUBROUTINE ARRINV(A,B,NNB,NWB,MW,I,J1,I2) 00013980 C 00013990 C ******************************************************************00014000 C * *00014010 C * S U B R O U T I N E A R R I N V *00014020 C * *00014030 C * THIS SUBROUTINE TAKES THE CHOLESKI DECOMPOSED MATRIX AND *00014040 C * DOES A BACK SOLUTION FOR THE BANDED AND BORDERED PORTION *00014050 C * OF THE VARIANCE CO-VARIANCE MATRIX. *00014060 C * *00014070 C ******************************************************************00014080 C 00014090 REAL*8 A(NWB,NNB),B(NNB),P 00014100 DIMENSION ITIME(17),IOUNIT(7) 00014110 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00014120 DO 140 II = J1,MW 00014130 I = I - 1 00014140 IK = 0 00014150 I1 = II 00014160 IF(II.GE.NWB) I1 = NW - 1 00014170 DO 10 J=1,I1 00014180 10 B(J)=0.0D0 00014190 DO 20 K=1,I1 00014200 DO 20 J=K,I1 00014210 20 B(J)=B(J)+A(K+1,I)*A(J-K+1,I+K) 00014220 DO 60 J=1,I1 00014230 M = I + J 00014240 IF(J.EQ.I1) GO TO 40 00014250 L = J + 1 00014260 DO 30 K = L,I1 00014270 KI = K + 1 00014280 30 B(J) = B(J) + A(KI,I)*A(KI-J,M) 00014290 40 IF(II.LT.NWB) GO TO 60 00014300 L = I1 + 1 00014310 IF((II-J).LT.I2) IK = IK - 1 00014320 DO 50 K = L,I2 00014330 KI = K + 1 00014340 50 B(J) = B(J) + A(KI,I)*A(KI+IK,M) 00014350 60 B(J) = -B(J)/A(1,I) 00014360 IF(II.LT.NWB) GO TO 110 00014370 I3 = I1 + 1 00014380 I4 = NWB - 1 00014390 DO 70 J=I3,I4 00014400 70 B(J)=0.0D0 00014410 IK=0 00014420 DO 80 K=1,I1 00014430 IF((II-K).LT.I2)IK=IK-1 00014440 DO 80 J=I3,I4 00014450 80 B(J)=B(J)+A(K+1,I)*A(J+1+IK,I+K) 00014460 DO 100 J = I3,I4 00014470 L = 1 00014480 JJ = J - I3 + 1 00014490 DO 90 K = I3,I4 00014500 B(J) = B(J) + A(K+1,I)*A(JJ,N+L) 00014510 L = L + 1 00014520 JJ = JJ - 1 00014530 IF(L.GT.(J-I3+1)) JJ = JJ + 2 00014540 IF(L.GT.(J-I3+1)) L = L - 1 00014550 90 CONTINUE 00014560 100 B(J) = -B(J)/A(1,I) 00014570 110 P = 0.0D0 00014580 K = I2 00014590 IF(II.LT.NWB) K = I1 00014600 DO 120 J = 1,K 00014610 120 P = P + B(J)*A(J+1,I) 00014620 A(1,I) = (1.0D0/A(1,I)-P)/A(1,I) 00014630 DO 130 J = 1,K 00014640 130 A(J+1,I) = B(J) 00014650 140 CONTINUE 00014660 RETURN 00014670 END 00014680 SUBROUTINE ABSOLU(TA,NK1,NK2,NK5,CFACT,PHIV) 00014690 C 00014700 C SUBROUTINE 'ABSOLU' COMPUTES THE ABSOLUTE ERROR ELLIPSES 00014710 C (IN METRES) AT EACH NETWORK STATION. THE ANGLE IS W.R.T. GEODETIC 00014720 C NORTH AT THE COMPUTATION POINT. 00014730 C 00014740 C INPUT: 00014750 C TA = INVERSE OF NORMAL EQUATIONS 00014760 C NK1 = NUMBER OF NETWORK STATIONS 00014770 C NK2 = NUMBER OF ELEMENTS IN SOLUTION VECTOR 00014780 C NK5 = BANDWIDTH + BORDERWIDTH 00014790 C CFACT = SCALING FACTOR (F-STAT) 00014800 C PHIV = VECTOR NK1 LONG WITH GEODETIC LATITUDES 00014810 C 00014820 C OUTPUT: LINE PRINTER ONLY 00014830 C 00014840 IMPLICIT REAL*8(A-H,O-Z) 00014850 DIMENSION TA(NK5,NK2),PHIV(NK1) 00014860 C 00014870 C PRIME VERTICAL AND MERIDIAN RADII OF CURVATURE 00014880 C 00014890 RNORM(X)=AR/DSQRT(1.0D0-E2*DSIN(X)**2) 00014900 RMER(X)=AR*(1.0D0-E2)/DSQRT(1.0D0-E2*DSIN(X)**2)**3 00014910 RHO=206264.8062470963D0**2 00014920 C 00014930 C CLARKE 1866 ELLIPSOID 00014940 C 00014950 AR=6378206.4D0 00014960 BR=6356583.8D0 00014970 E2=(AR*AR-BR*BR)/AR/AR 00014980 J=1 00014990 PRINT 1000 00015000 PRINT 1010 00015010 I=1 00015020 K = NK1 - NK2/2 00015030 C 00015040 C COMPUTE RADII OF CURVATURE AT STATION J 00015050 C 00015060 10 RM=RMER(PHIV(J)) 00015070 RN=RNORM(PHIV(J)) 00015080 CP=DCOS(PHIV(J)) 00015090 C 00015100 C QYY = SIGMA**PHI IN METRES**2 00015110 C QXY = SIGMA(PHI,LAMDA) IN METRES**2 00015120 C QXX = SIGMA**2 LAMDA IN METRES**2 00015130 C 00015140 QYY=TA(1,I)/RHO*RM*RM 00015150 QXY=TA(2,I)/RHO*RM*RN*CP 00015160 I=I+1 00015170 QXX=TA(1,I)*RN*CP/RHO*RN*CP 00015180 C 00015190 C A = SEMI-MAJOR AXIS IN METRES (SCALED) 00015200 C G = SEMI-MINOR AXIS IN METRES (SCALED) 00015210 C C = ANGLE FROM NORTH IN DEGREES 00015220 C 00015230 A=DSQRT((QXX+QYY)/2.0D0+DSQRT((QYY-QXX)**2/4.0D0+QXY**2)) 00015240 G=DSQRT((QXX+QYY)/2.0D0-DSQRT((QYY-QXX)**2/4.0D0+QXY**2)) 00015250 C= 2.0*QXY 00015260 D= (QYY-QXX) 00015270 C=-DATAN2(C,D)*57.2957795D0/2.0D0 00015280 A = A*CFACT 00015290 G = G*CFACT 00015300 K=K+1 00015310 PRINT 1020,K,A,G,C 00015320 I=I+1 00015330 J=J+1 00015340 IF((J/26*26).EQ.J) PRINT 1000 00015350 IF((J/26*26).EQ.J) PRINT 1010 00015360 IF(I.LT.NK2) GO TO 10 00015370 RETURN 00015380 1000 FORMAT('1',/////,39X,'ABSOLUTE 95 % ERROR ELLIPSES',//) 00015390 1010 FORMAT(15X,'STATION NO.',10X,'SEMI MAJOR AXIS',10X,'SEMI MINOR AXI00015400 .S',15X,'ANGLE'/40X,'(METRES)',17X,'(METRES)'///) 00015410 1020 FORMAT(I22,F24.4,F25.4,F25.4,/) 00015420 END 00015430 SUBROUTINE RELATI(TA,NK1,NK2,NK5,CFACT,PHIV) 00015440 C 00015450 C SUBROUTINE 'RELATI' COMPUTES THE RELATIVE ERROR ELLIPSE (IN 00015460 C METRES) BETWEEN TWO POINTS. THE ANGLE IS W.R.T. GEODETIC NORTH 00015470 C OF THE MID-POINT OF THE LINE 00015480 C 00015490 C INPUT: 00015500 C TA = INVERSE OF NORMAL EQUATIONS (NK5 * NK2) 00015510 C NK1 = NUMBER OF NETWORK POINTS 00015520 C NK2 = NUMBER OF ELEMENTS IN SOLUTION VECTOR 00015530 C NK5 = BANDWIDTH PLUS BORDERWIDTH 00015540 C CFACT = SCALING FACTOR FOR 95 % (F DIST) 00015550 C PHIV = VECTOR NK1 LONG WITH LATITUDES OF NETWORK STATIONS 00015560 C 00015570 C THE NUMBER OF ELLIPSES REQUIRED AND BETWEEN WHICH STATIONS IS READ 00015580 C FROM CARDS (SEE COMMENTS IN MAIN PROGRAM FOR FORMAT) 00015590 C 00015600 C OUTPUT: LINE PRINTER ONLY 00015610 C 00015620 IMPLICIT REAL*8(A-H,O-Z) 00015630 DIMENSION TA(NK5,NK2),PHIV(NK1),ITIME(17),IOUNIT(7) 00015640 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00015650 C 00015660 C PRIME VERTICAL AND MERIDIAN RADII OF CURVATURE OF REFERENCE ELLIPSOID00015670 C AT LATITUDE X 00015680 C 00015690 RNORM(X)=AR/DSQRT(1.0D0-E2*DSIN(X)**2) 00015700 RMER(X)=AR*(1.0D0-E2)/DSQRT(1.0D0-E2*DSIN(X)**2)**3 00015710 C 00015720 C CLARKE 1866 ELLIPSOID 00015730 C 00015740 AR=6378206.4D0 00015750 BR=6356583.8D0 00015760 E2=(AR*AR-BR*BR)/AR/AR 00015770 RHO=206264.8062470963D0**2 00015780 PRINT 1000 00015790 PRINT 1010 00015800 MMN = 2*NK1 - NK2 00015810 C 00015820 C NELIPS = NUMBER OF ELLIPSES REQUIRED (READ FROM DATA CARD) 00015830 C 00015840 READ 1020,NELIPS 00015850 IF(NELIPS.EQ.0)GO TO 30 00015860 K = 0 00015870 C 00015880 C READ FROM DATA CARDS STATION NUMBERS OF POINTS FOR WHICH ELLIPSE 00015890 C REQUIRED (NOTE SEQUENCE NUMBERS) 00015900 C 00015910 10 READ 1030,NSTA1,NSTA2 00015920 C 00015930 C NSTA1 MUST BE LESS THAN NSTA2 00015940 C 00015950 IF(NSTA1.LT.NSTA2)GO TO 20 00015960 NN=NSTA1 00015970 NSTA1=NSTA2 00015980 NSTA2=NN 00015990 C 00016000 C COMPUTE ALL RADII OF CURVATURE 00016010 C 00016020 20 RMI=RMER(PHIV(NSTA1)) 00016030 RNI=RNORM(PHIV(NSTA1)) 00016040 RMJ=RMER(PHIV(NSTA2)) 00016050 RNJ=RNORM(PHIV(NSTA2)) 00016060 RMEAN=(RMI*RNI*DCOS(PHIV(NSTA1))+RMJ*RNJ*DCOS(PHIV(NSTA2)))/2.0D0/00016070 .RHO 00016080 C 00016090 C A,BB,C = 2 * 2 SUBMATRIX FOR NSTA1 - SCALED TO METRES 00016100 C 00016110 I = 2*NSTA1 - 1 - MMN 00016120 A=TA(1,I)/RHO*RMI*RMI 00016130 BB=TA(2,I)/RHO*RMI*RNI*DCOS(PHIV(NSTA1)) 00016140 I = I + 1 00016150 C=TA(1,I)/RHO*(RNI*DCOS(PHIV(NSTA1)))**2 00016160 I = I - 1 00016170 J = (NSTA2 - NSTA1)*2 + 1 00016180 IF(J.GT.NW.AND.I.LT.(N-NW+1))J=J-NSPACE+I 00016190 C 00016200 C D,E,F,G = 2 * 2 COVARIANCE MATRIX NSTA1,NSTA2 - SCALED TO METRES 00016210 C 00016220 D=TA(J,I)*RMEAN 00016230 J = J + 1 00016240 E=TA(J,I)*RMEAN 00016250 I = I + 1 00016260 IF(J.LE.NW) J = J - 1 00016270 F=TA(J,I)*RMEAN 00016280 J = J - 1 00016290 G=TA(J,I)*RMEAN 00016300 C 00016310 C H,U,Q = 2 * 2 SUBMATRIX FOR NSTA2 - SCALED TO METRES 00016320 C 00016330 I = 2*NSTA2 - 1 - MMN 00016340 H=TA(1,I)/RHO*RMJ*RMJ 00016350 U=TA(2,I)/RHO*RMJ*RNJ*DCOS(PHIV(NSTA2)) 00016360 I = I + 1 00016370 Q=TA(1,I)/RHO*(RNJ*DCOS(PHIV(NSTA2)))**2 00016380 QDXDX = H - 2.0D0*D + A 00016390 QDYDY = Q - 2.0D0*F + C 00016400 QDXDY = U + BB - E - G 00016410 C 00016420 C A = SEMI-MAJOR AXIS 00016430 C BB = SEMI-MINOR AXIS 00016440 C C = ANGLE FROM NORTH IN DEGREES 00016450 C 00016460 A = DSQRT((QDXDX + QDYDY)/2.0D0 + DSQRT((QDYDY - QDXDX)**2/4.0D0 00016470 1+ QDXDY**2)) 00016480 BB = DSQRT((QDXDX + QDYDY)/2.0D0 - DSQRT((QDYDY - QDXDX)**2/4.0D0 00016490 1+ QDXDY**2)) 00016500 C = 2.0D0*QDXDY 00016510 D= (QDYDY-QDXDX) 00016520 C = -DATAN2(C,D)*57.2957795D0/2.0D0 00016530 C 00016540 C SCALE AXIS OF ELLIPSE TO 95 % 00016550 C 00016560 A = A*CFACT 00016570 BB = BB*CFACT 00016580 PRINT 1040,NSTA1,NSTA2,A,BB,C 00016590 K = K + 1 00016600 IF(K.LT.NELIPS)GO TO 10 00016610 RETURN 00016620 C 00016630 C NO ERROR ELLIPSES REQUIRED 00016640 C 00016650 30 PRINT 1050 00016660 RETURN 00016670 1000 FORMAT('1',/////,39X,'RELATIVE 95 % ERROR ELLIPSES',//) 00016680 1010 FORMAT(10X,'STATION TO STATION',10X,'SEMI MAJOR AXIS',10X,'SEMI MI00016690 .NOR AXIS',15X,'ANGLE'/42X,'(METRES)',17X,'(METRES)'///) 00016700 1020 FORMAT(I4) 00016710 1030 FORMAT(2I5) 00016720 1040 FORMAT(I17,I13,F20.4,F25.4,F25.4) 00016730 1050 FORMAT('-',10X,'NO RELATIVE ERROR ELLIPSES ARE REQUIRED') 00016740 END 00016750 SUBROUTINE PXADD(NK1,NK2,NK5,NFP,NFPT,INEW,IVEC,ISTAFP,WTFP,PX,TA)00016760 C 00016770 C SUBROUTINE 'PXADD' ADDS THE MATRIX 'WTFP' ONTO THE NORMAL 00016780 C EQUATIONS (TA) 00016790 C 00016800 C INPUT: 00016810 C NK1 = NUMBER OF NETWORK STATIONS 00016820 C NK2 = NUMBER OF ELEMENTS IN SOLUTION VECTOR 00016830 C NK5 = BANDWIDTH + BORDERWIDTH 00016840 C NFP = NUMBER OF WEIGHTED POINTS 00016850 C NFPT = 2*NFP 00016860 C INEW = SCRATCH VECTOR NFP LONG 00016870 C IVEC = SCRATCH VECTOR NFP LONG 00016880 C ISTAFP = VECTOR NFP LONG WITH SEQUENCE NUMBERS OF WEIGHTED 00016890 C POINTS 00016900 C WTFP = WEIGHT (PX) MATRIX (NFPT * NFPT) 00016910 C PX = SCRATCH MATRIX (NFPT * NFPT) 00016920 C TA = NORMAL EQUATIONS (NK5 * NK2) 00016930 C 00016940 C OUTPUT: 00016950 C TA = NORMAL EQUATIONS WITH PX ADDED 00016960 C 00016970 IMPLICIT REAL*8 (A-H,O-Z) 00016980 DIMENSION TA(NK5,NK2),ISTAFP(NFP),WTFP(NFPT,NFPT),INEW(NFP), 00016990 1IVEC(NFP),PX(NFPT,NFPT),ITIME(17),IOUNIT(7) 00017000 COMMON NW,NB,N,NBB,NSPACE,IOUNIT,ITIME 00017010 C 00017020 C RE-ORDER ISTAFP AND WTFP SO THAT SEQUENCE NUMBERS ARE IN 00017030 C ASCENDING ORDER 00017040 C 00017050 CALL PXORD(NFP,NFPT,ISTAFP,WTFP,PX,INEW,IVEC) 00017060 MMN=2*NK1-NK2 00017070 KK=0 00017080 DO 40 L=1,NFP 00017090 K=L+1 00017100 LL=L+1 00017110 IF(LL.GT.NFP)LL=NFP 00017120 DO 30 M=LL,NFP 00017130 IF(M.NE.LL)GO TO 10 00017140 C 00017150 C ADD ON MAIN DIAGONAL 2 * 2 SUBMATRIX (UPPER HALF) 00017160 C 00017170 I=2*ISTAFP(L)-1-MMN 00017180 TA(1,I)=TA(1,I)+WTFP(L+KK,L+KK) 00017190 TA(2,I)=TA(2,I)+WTFP(L+KK,L+KK+1) 00017200 I=I+1 00017210 TA(1,I)=TA(1,I)+WTFP(L+KK+1,L+KK+1) 00017220 I=I-1 00017230 IF(K.GT.NFP)RETURN 00017240 C 00017250 C ADD ON 2 * 2 SUBMATRICES WORKING ACROSS THE ROWS 00017260 C 00017270 IF(M.EQ.LL)GO TO 20 00017280 10 I=2*ISTAFP(L)-1-MMN 00017290 20 J=(ISTAFP( M)-ISTAFP(L))*2+1 00017300 IF(J.GT.NW.AND.I.LT.(N-NW+1))J=J-NSPACE+I 00017310 TA(J,I)=TA(J,I)+WTFP(L+KK,2*M-1) 00017320 J=J+1 00017330 TA(J,I)=TA(J,I)+WTFP(L+KK,2*M) 00017340 I=I+1 00017350 IF(J.LE.NW)J=J-1 00017360 TA(J,I)=TA(J,I)+WTFP(L+KK+1,2*M) 00017370 J=J-1 00017380 30 TA(J,I)=TA(J,I)+WTFP(L+KK+1,2*M-1) 00017390 40 KK=KK+1 00017400 STOP 300 00017410 C 00017420 C ERROR EXIT - SUBSCRIPTING ERROR 00017430 C 00017440 END 00017450 SUBROUTINE PXORD(N,N2,ISTA,PXMAT,PX,INEW,IVEC) 00017460 C 00017470 C SUBROUTINE PXORD REORDERS STATIONS IN INCREASING SEQUENCE NUMBERS AND00017480 C THE WEIGHT MATRIX FOR SUBROUTINE PXADD 00017490 C 00017500 C INPUT: 00017510 C N = NUMBER OF WEIGHTED STATIONS 00017520 C N2 = 2*N 00017530 C ISTA = VECTOR N LONG WITH UNORDERED STATION NUMBERS 00017540 C PXMAT = UNORDERED PX MATRIX (N2 * N2) 00017550 C PX = SCRATCH MATRIX (N2 * N2) 00017560 C INEW = SCRATCH VECTOR N LONG 00017570 C IVEC = SCRATCH VECTOR N LONG 00017580 C 00017590 C OUTPUT: 00017600 C ISTA = ORDERED STATION NUMBERS (ASCENDING) 00017610 C PXMAT = ORDER%D PX MATRIX 00017620 C 00017630 IMPLICIT REAL*8 (A-H,O-Z) 00017640 DIMENSION PXMAT(N2,N2),PX(N2,N2),ISTA(N),INEW(N),IVEC(N) 00017650 DATA IDUM/999999/ 00017660 C 00017670 C ONLY 1 STATION - NO RE-ORDERING REQUIRED 00017680 C 00017690 IF(N.EQ.1)RETURN 00017700 C 00017710 C CHECK IF RE-ORDERING REQUIRED - IF NO RETURN 00017720 C 00017730 DO 10 I=2,N 00017740 K=I-1 00017750 DO 10 J=1,K 00017760 IF(ISTA(I).LT.ISTA(J))GO TO 20 00017770 10 CONTINUE 00017780 RETURN 00017790 C 00017800 C REORDER SEQUENCE NUMBERS 00017810 C 00017820 20 K=1 00017830 30 IMIN=ISTA(1) 00017840 JMIN=1 00017850 DO 40 I=2,N 00017860 IF(ISTA(I).GE.IMIN)GO TO 40 00017870 JMIN=I 00017880 IMIN=ISTA(I) 00017890 40 CONTINUE 00017900 INEW(K)=IMIN 00017910 IVEC(K)=JMIN 00017920 ISTA(JMIN)=IDUM 00017930 K=K+1 00017940 IF(K.LE.N)GO TO 30 00017950 C 00017960 C PUT REORDERED STATIONS BACK INTO INPUT MATRIX 00017970 C 00017980 DO 50 I=1,N 00017990 ISTA(I)=INEW(I) 00018000 50 INEW(I)=IVEC(I) 00018010 C 00018020 C INTERCHANGE COLUMNS OF WEIGHT MATRIX (PX) AND STORE IN TEMPORARY 00018030 C ARRAY 00018040 C 00018050 DO 60 I=1,N 00018060 K=2*I 00018070 L=K-1 00018080 KK=2*IVEC(I) 00018090 LL=KK-1 00018100 DO 60 J=1,N2 00018110 PX(J,L)=PXMAT(J,LL) 00018120 PX(J,K)=PXMAT(J,KK) 00018130 60 CONTINUE 00018140 C 00018150 C INTERCHANGE ROWS OF TEMPORARY ARRAY AND PUT BACK INTO WEIGHT (PX) 00018160 C MATRIX 00018170 C 00018180 DO 70 I=1,N 00018190 K=2*I 00018200 L=K-1 00018210 KK=2*INEW(I) 00018220 LL=KK-1 00018230 DO 70 J=1,N2 00018240 PXMAT(L,J)=PX(LL,J) 00018250 PXMAT(K,J)=PX(KK,J) 00018260 70 CONTINUE 00018270 RETURN 00018280 END 00018290