C 00000010 C 00000020 C 00000030 C 00000040 C CONGA2 -- CONSTRAINED VERSION OF ANGEOID 00000050 C PROGRAMME TO COMPUTE GEOIDAL HEIGHTS FROM DEFLECTIONS OF THE VER00000060 C USING A LEAST-SQUARES SURFACE-FITTING TECHNIQUE 00000070 C THIS VERSION INCORPORATES HEIGHT CONSTRAINTS. COEFFICIENT C00 IS 00000080 C SOLVED FOR, NOT HELD FIXED. JUNE, 1974 00000090 C 00000100 C 00000110 C INPUT REQUIREMENTS ******* 00000120 C 00000130 C CARD ONE 00000140 C HEADER CARD -- CONTAINS PROJECT TITLE -- 20A4 00000150 C 00000160 C CARD TWO 00000170 C COORDINATES -- LAT. & LONG. , IN DEGREES, OF LOCAL ORIGIN --- 200000180 C 00000190 C CARD THREE 00000200 C LIMIT OF AREA OF INTEREST ( DX,DY,MEASURED FROM ORIGIN IN KMS.) 00000210 C AND SCALE FACTOR --- 2F12.2 00000220 C 00000230 C CARD FOUR 00000240 C DEGREE OF POLYNOMIAL TO BE USED --- I2 00000250 C 00000260 C CARD FIVE 00000270 C NP =NUMBER OF CONSTRAINTS --- I4 00000280 C 00000290 C NEXT (NP*NP)/10 CARDS: 00000300 C VARIANCE-COVARIANCE MATRIX OF HEIGHTS - IN ROW SEQUENCE - TEN ELE00000310 C PER CARD --- 10F8.3 00000320 C 00000330 C NEXT NP CARDS: 00000340 C HEIGHT CONSTRAINT INFORMATION. PER CARD: LAT., LONG.(DEGS), 00000350 C HEIGHT (METRES) --- 3F8.3 00000360 C 00000370 C NEXT CARD 00000380 C LIMITS FOR GRID VALUES - TWO LAT. VALUES, TWO LONG. VALUES, GRID00000390 C INTERVAL, ALL IN DEGREES ---- 5F8.3 00000400 C 00000410 C NEXT CARD 00000420 C COLS. 1-2 ANY INTEGER IF GRID VALUES TO BE PUNCHED. BLANK IF N00000430 C BE PUNCHED. COLS. 3-4 ANY INTEGER IF COEFFICIENTS AND 00000440 C VARIANCE-COVARIANCE MATRIX TO BE STORED ON DISC. OTHERWISE - BLAN00000450 C 00000460 C NEXT CARD 00000470 C NUMBER OF DATA SETS TO BE USED --- I2 00000480 C 00000490 C NEXT CARDS 00000500 C DATA SET REFERENCE NOS. (MAGNETIC TAPE OR DISC FILE) - ONE PER 00000510 C CARD ---- I2 00000520 C 00000530 C 00000540 IMPLICIT REAL*8 (A-H,O-Z) 00000550 DIMENSION NUM(1000),X(1000),Y(1000),WX(1000),WE(1000) 00000560 DIMENSION XSI(1000),ETA(1000) 00000570 DIMENSION NUMM(500),XX(500),YY(500),WWX(500),WWE(500) DIMENSION XXSI(500),EETA(500) DIMENSION A(100,100),FL(100),FAL(100) 00000580 DIMENSION IDSS(10) 00000590 DIMENSION NAME(20) 00000600 COMMON /BLOCK 1/ FAL,NK 00000610 COMMON /BLOCK 2/ A,NKA 00000620 COMMON /BLOCK 3/ SCALE 00000630 DATA FL/100*0.D0/ 00000640 101 FORMAT(2F12.2) 00000650 102 FORMAT(2F11.6) 00000660 103 FORMAT(5F8.2) 00000670 104 FORMAT(//,20X,'NO. OF MERIDIAN DEFLECTIONS = ',I4,/,20X,'NO. OF PR00000680 *IME VERTICAL DEFLECTIONS = ',I4,//) 00000690 105 FORMAT(3F8.3) 00000700 106 FORMAT(2I2) 00000710 107 FORMAT(I2) 00000720 108 FORMAT(20A4) 00000730 200 FORMAT(/,30X,D16.8) 00000740 201 FORMAT('1',30X,'COEFFICIENTS OF',/,20X,'BEST-FITTING POLYNOMIAL',/00000750 *//) 00000760 202 FORMAT(/,10X,I4,F15.2,2(10X,F12.2,5X,F10.2)) 00000770 203 FORMAT('1',9X,'STATION',7X,'GEOIDAL',15X,'MERIDIAN',26X,'PRIME VER00000780 *TICAL',/,10X,'NUMBER',8X,'HEIGHT',14X,'DEFLECTION',6X,'RESIDUAL',100000790 *3X,'DEFLECTION',6X,'RESIDUAL',///) 00000800 204 FORMAT(2(10X,F7.2),2F15.2,2F20.2,/) 00000810 205 FORMAT(12X,'LAT.',12X,'LONG.',9X,'HEIGHT',9X,'STD. DEV.',14X,'XSI'00000820 *,17X,'ETA',//) 00000830 206 FORMAT(//,20X,'ESTIMATED VARIANCE FACTOR=',E15.6,//) 00000840 207 FORMAT(//,20X,'DEGREE OF APPROXIMATING POLYNOMIAL=',2X,I4,//) 00000850 208 FORMAT('0',40X,'LOCAL ORIGIN',//,20X,'LATITUDE= ',F11.6,2X,'LONGIT00000860 *UDE= ',F12.6,//,20X,'LIMIT OF AREA COVERED - DX, DY =',F14.2,///) 00000870 209 FORMAT(30X,'LIMITING PRECISION :',I8,2X,'SIGNIFICANT FIGURES',//) 00000880 212 FORMAT(/,20X,'SCALE FACTOR =',F20.2,/) 00000890 213 FORMAT(/,10X,'TIME = ',I6) 00000900 214 FORMAT('1',35X,20A4,//) 00000910 215 FORMAT(///,40X,'GRID VALUES',//) 00000920 1023 FORMAT(I8,2F12.6,4F8.3) 1030 FORMAT(4F8.3) C 00000930 C CONSTANTS 00000940 C 00000950 IRC=5 00000960 IPP=6 00000970 R=6370990. 00000980 RH=206264.80625 00000990 RH2=RH**2 00001000 RHO=206264.80625/3600. 00001010 ACX=1.D0 00001020 ACE=1.D0 00001030 ACT=2.D0 MM=0 MKK=0 00001040 IMJ=0 00001050 INJ=0 00001060 DO 310 I=1,100 00001070 DO 310 J=1,100 00001080 A(I,J)=0.D0 00001110 C 00001090 310 CONTINUE 00001100 C READ HEADER CARD 00001120 READ(IRC,108) NAME WRITE(IPP,214) NAME 00001130 C 00001140 C READ COORDINATES OF LOCAL ORIGIN 00001150 C 00001160 READ(IRC,102) OLAT,OLONG 00001170 C 00001180 C READ IN LIMIT OF AREA COVERED (IN KMS.) AND SCALE FACTOR 00001190 C 00001200 READ(IRC,101) ALIMIT,SCALE 00001210 WRITE(IPP,208) OLAT,OLONG,ALIMIT 00001220 ALIMIT=ALIMIT*1000.D0 00001230 WRITE(IPP,212) SCALE 00001240 AHT=0.D0 00001250 C 00001260 C READ MAXIMUM DEGRE TO BE USED 00001270 C 00001280 READ(IRC,107) NAN 00001290 NN=NAN 00001300 NK=NN+1 00001310 WRITE(IPP,207) NN 00001320 NKA=NK*NK 00001330 NAK=NKA-1 00001340 C 00001350 C FORM NORMAL EQUATIONS FOR HEIGHT CONSTRAINTS. 00001360 C 00001370 CALL CON(AHT,NK,SCALE,OLONG,OLAT,FL,NP) 00001380 C 00001390 C READ LIMITS OF GRID EVALUATION AND GRID INTERVAL 00001400 C 00001410 READ(IRC,1030) QP1,QPP1,QL1,QLL1 READ(IRC,103) QP,QPP,QL,QLL,DEGT 00001420 IQP=IABS(IDINT((QPP-QP)/DEGT))+1 00001430 IQL=IABS(IDINT((QLL-QL)/DEGT))+1 00001440 C 00001450 C READ IN CODE FOR PUNCHING OF GRID VALUES 00001460 C 00001470 READ(IRC,106) ICODE,ICD 00001480 C 00001490 C READ IN NUMBER OF DATA SETS 00001500 C 00001510 READ(IRC,107) MJ 00001520 DO 57 MI=1,MJ 00001530 C 00001540 C READ IN DATA SET REFERENCE NUMBERS 00001550 C 00001560 READ(IRC,107) IDS 00001570 IDSS(MI)=IDS 00001580 IMM=0 INN=0 I=1 00001600 IF(IDS.EQ.IRC) GOTO 1021 C 00001610 C READ IN DATA AND TRANSFORM TO PLANE SYSTEM. 00001620 C 00001630 1 READ(IDS,END=2) NUM(I),X(I),Y(I),XSI(I),ETA(I),WX(I),WE(I) 00001640 IF(X(I).LE.QP1.OR.X(I).GT.QPP1) GOTO1 IF(Y(I).LE.QL1.OR.Y(I).GT.QLL1) GOTO1 GOTO 1022 1021 READ(IDS,1023) NUM(I),X(I),Y(I),XSI(I),ETA(I),WX(I),WE(I) NUMM(I)=NUM(I) XX(I)=X(I) YY(I)=Y(I) XXSI(I)=XSI(I) EETA(I)=ETA(I) WWX(I)=WX(I) WWE(I)=WE(I) 1022 CONTINUE IF(WX(I).LT.0.001) IMM=IMM+1 00001650 IF(WE(I).LT.0.001) INN=INN+1 00001660 IF(WX(I).LT.0.001) WX(I)=1.0D+20 00001670 IF(WE(I).LT.0.001) WE(I)=1.0D+20 00001680 WX(I)=1.D0/(WX (I)**2) 00001690 WE(I)=1.D0/(WE(I)**2) 00001700 IF(DABS(XSI(I)).GT.50.0) XSI(I)=0.D0 00001710 IF(DABS(ETA(I)).GT.50.0) ETA(I)=0.D0 00001720 Y(I)=(Y(I)-OLONG)/RHO*R*DCOS(X(I)/RHO) 00001730 X(I)=(X(I)-OLAT)/RHO*R 00001740 IF(IDS.EQ.IRC) GOTO 1025 IF(DABS(X(I)).GT.ALIMIT) GOTO1 00001750 IF(DABS(Y(I)).GT.ALIMIT) GOTO1 00001760 GOTO 1026 1025 IF(DABS(X(I)).GT.ALIMIT) GOTO 1021 IF(DABS(Y(I)).GT.ALIMIT) GOTO 1021 1026 CONTINUE IF(X(I).EQ.0.D0) X(I)=1.D0 00001770 IF(Y(I).EQ.0.D0) Y(I)=1.D0 00001780 X(I)=X(I)/SCALE 00001790 Y(I)=Y(I)/SCALE 00001800 XSI(I)=XSI(I)/RH 00001810 ETA(I)=ETA(I)/RH 00001820 I=I+1 00001830 IF(I.EQ.495) GOTO2 IF(IDS.EQ.IRC) GOTO 1021 GOTO1 00001840 2 CONTINUE 00001850 MM=I-1 00001860 IF(MM.EQ.0) GOTO57 00001870 C 00001880 C FORMATION OF MATRIX OF NORMAL EQUATIONS 00001890 C 00001900 CALL ERRSET(208,300,-1,1) 00001910 CALL CPUTIM(ITIM,IREM) 00001920 IT=ITIM/10000 00001930 WRITE(IPP,213) IT 00001940 NPP=(NAK**2+NAK)/2 00001950 DO 62 K=1,NPP 00001960 CALL INDEX(NN,K,M,N,I,J,IR,IS) 00001970 AMN=0.D0 00001980 BMN=0.D0 00001990 ISS=IS+I-2 00002000 IF(ISS.LT.0) GOTO66 00002010 IF(IS.EQ.0) GOTO66 00002020 DO 64 II=1,MM 00002030 AMN=AMN+WX(II)*X(II)**ISS*Y(II)**(IR+J) 00002040 64 CONTINUE 00002050 66 CONTINUE 00002060 ISR=IR+J-2 00002070 IF(ISR.LT.0) GOTO69 00002080 IF(IR.EQ.0) GOTO69 00002090 DO 68 II=1,MM 00002100 BMN=BMN+WE(II)*X(II)**(IS+I)*Y(II)**ISR 00002110 68 CONTINUE 00002120 69 CONTINUE 00002130 A(M+1,N+1)=A(M+1,N+1)+I*IS*AMN+J*IR*BMN 00002140 62 CONTINUE 00002150 C 00002160 C FORMATION OF OBSERVATION VECTOR 00002170 C 00002180 N=1 00002190 IS=0 00002200 DO 56 IF=1,NK 00002210 IR=0 00002220 DO 55 IE=1,NK 00002230 IF(IS.EQ.0.AND.IR.EQ.0) GOTO54 00002240 ISK=IS-1 00002250 IF(ISK.LT.0) GOTO52 00002260 FLN=0.D0 00002270 DO 50 II=1,MM 00002280 FLN=FLN+XSI(II)*WX(II)*X(II)**ISK*Y(II)**IR 00002290 50 CONTINUE 00002300 52 CONTINUE 00002310 IRK=IR-1 00002320 IF(IRK.LT.0) GOTO53 00002330 FALN=0.D0 00002340 DO 51 II=1,MM 00002350 FALN=FALN+ETA(II)*WE(II)*X(II)**IS*Y(II)**IRK 00002360 51 CONTINUE 00002370 53 CONTINUE 00002380 FL(N+1)=FL(N+1)-IS*FLN-IR*FALN 00002390 N=N+1 00002400 54 CONTINUE 00002410 IR=IR+1 00002420 55 CONTINUE 00002430 IS=IS+1 00002440 56 CONTINUE 00002450 CALL CPUTIM(ITIM,IREM) 00002460 IT=ITIM/10000 00002470 WRITE(IPP,213) IT 00002480 MKK=MKK+MM 00002490 IMJ=IMJ+MM-IMM 00002500 INJ=INJ+MM-INN 00002510 WRITE(IPP,104) IMJ,INJ 00002520 57 CONTINUE DO 471 I=1,NKA 00002540 DO 471 J=1,NKA 00002550 IF(J.GE.I) GOTO472 00002560 A(I,J)=A(J,I) 00002570 472 CONTINUE 00002580 471 CONTINUE 00002590 TEST1=A(NKA,NKA) 00002600 DO 800 M=1,NKA 00002610 800 A(M,M)=A(M,M)+4.D0*NKA*0.54D-78 00002620 C 00002630 C COMPUTE AND PRINT VECTOR OF COEFFICIENTS. 00002640 C 00002650 C 00002660 C ...CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 00002670 C 00002680 NA=NKA 00002690 DETA = A(1,1) 00002700 A(1,1) = DSQRT(A(1,1)) 00002710 DO 910 I=2,NA 00002720 910 A(I,1)=A(I,1)/A(1,1) 00002730 DO 500 J = 2,NA 00002740 SUM = 0. 00002750 J1 = J - 1 00002760 DO 920 K=1,J1 00002770 920 SUM=SUM+A(J,K)*A(J,K) 00002780 DETA = DETA * (A(J,J) - SUM) 00002790 A(J,J) = DSQRT(A(J,J) - SUM) 00002800 IF(J .EQ. NA) GOTO 500 00002810 J2 = J + 1 00002820 DO 400 I = J2,NA 00002830 SUM = 0. 00002840 DO 300 K = 1,J1 00002850 300 SUM = SUM + A(I,K) * A(J,K) 00002860 400 A(I,J) = (A(I,J) - SUM) / A(J,J) 00002870 500 CONTINUE 00002880 C 00002890 C ...INVERSION OF LOWER TRIANGULAR MATRIX 00002900 C 00002910 DO 600 I = 1,NA 00002920 600 A(I,I) = 1D0 / A(I,I) 00002930 N1 = NA - 1 00002940 DO 980 J=1,N1 00002950 J2 = J + 1 00002960 DO 980 I=J2,NA 00002970 SUM = 0. 00002980 I1 = I - 1 00002990 DO 700 K = J,I1 00003000 700 SUM = SUM + A(I,K) * A(K,J) 00003010 980 A(I,J)=-A(I,I)*SUM 00003020 C 00003030 C ...CONSTRUCTION OF INVERSE OF INPUT MATRIX 00003040 C 00003050 DO 1300 J = 1,NA 00003060 IF(J .EQ. 1) GOTO 1000 00003070 J1 = J - 1 00003080 DO 900 I = 1,J1 00003090 900 A(I,J) = A(J,I) 00003100 1000 DO 1200 I = J,NA 00003110 SUM = 0. 00003120 DO 1100 K = I,NA 00003130 1100 SUM = SUM + A(K,I) * A(K,J) 00003140 1200 A(I,J) = SUM 00003150 1300 CONTINUE 00003160 CALL CPUTIM(ITIM,IREM) 00003170 IT=ITIM/10000 00003180 WRITE(IPP,213) IT 00003190 WRITE(6,805) DETA 00003200 805 FORMAT('1',30X,'DETERMINANT= ',D15.6,//) 00003210 TEST2=A(NKA,NKA) 00003220 TEA=TEST1*TEST2 00003230 ITEA=IDINT(16.8D0-DLOG10(TEA)) 00003240 WRITE(IPP,209) ITEA 00003250 DO 510 I=1,NKA 00003260 SUM=0.D0 00003270 DO 511 K=1,NKA 00003280 511 SUM=SUM+A(I,K)*FL(K) 00003290 510 FAL(I)=SUM 00003300 WRITE(IPP,201) 00003310 DO 70 I=1,NKA 00003320 WRITE(IPP,200) FAL(I) 00003330 70 CONTINUE 00003340 IF(ICD.EQ.0) GOTO75 00003350 WRITE(20) FAL 00003360 75 CONTINUE 00003370 C 00003380 C COMPUTE AND PRINT GEOIDAL UNDULATIONS , DEFLECTIONS, RESIDUALS. 00003390 C 00003400 WRITE(IPP,203) 00003410 VARF=0. 00003420 KIA=0 00003430 REWIND 4 DO 58 MI=1,MJ 00003440 IDS=IDSS(MI) 00003450 I=1 00003470 IF(IDS.EQ.IRC) GOTO 1041 3 READ(IDS,END=4) NUM(I),X(I),Y(I),XSI(I),ETA(I),WX(I),WE(I) 00003480 IF(X(I).LE.QP1.OR.X(I).GT.QPP1) GOTO3 IF(Y(I).LE.QL1.OR.Y(I).GT.QLL1) GOTO3 GOTO 1042 1041 NUM(I)=NUMM(I) X(I)=XX(I) Y(I)=YY(I) XSI(I)=XXSI(I) ETA(I)=EETA(I) WX(I)=WWX(I) WE(I)=WWE(I) 1042 CONTINUE IF(WX(I).LT.0.001) WX(I)=1.0D+20 00003490 IF(WE(I).LT.0.001) WE(I)=1.0D+20 00003500 WX(I)=1.D0/(WX(I)**2) 00003510 WE(I)=1.D0/(WE(I)**2) 00003520 IF(DABS(XSI(I)).GT.50.0) XSI(I)=0.D0 00003530 IF(DABS(ETA(I)).GT.50.0) ETA(I)=0.D0 00003540 Y(I)=(Y(I)-OLONG)/RHO*R*DCOS(X(I)/RHO) 00003550 X(I)=(X(I)-OLAT)/RHO*R 00003560 IF(IDS.EQ.IRC) GOTO 1050 IF(DABS(X(I)).GT.ALIMIT) GOTO3 00003570 IF(DABS(Y(I)).GT.ALIMIT) GOTO3 00003580 GOTO 1051 1050 IF(DABS(X(I)).GT.ALIMIT) GOTO 1041 IF(DABS(Y(I)).GT.ALIMIT) GOTO 1041 1051 CONTINUE IF(X(I).EQ.0.D0) X(I)=1.D0 00003590 IF(Y(I).EQ.0.D0) Y(I)=1.D0 00003600 X(I)=X(I)/SCALE 00003610 Y(I)=Y(I)/SCALE 00003620 XSI(I)=XSI(I)/RH 00003630 ETA(I)=ETA(I)/RH 00003640 KIA=KIA+1 00003650 AA=KIA 00003660 LL=KIA/25 00003670 IF((AA/25-LL).LE.1E-07)WRITE(IPP,203) 00003680 XS=PXSI(X(I),Y(I)) 00003690 ET=PETA(X(I),Y(I)) 00003700 VARF=VARF+WX(I)*(XS-XSI(I))**2+WE(I)*(ET-ETA(I))**2 00003710 GEOID=SCALE*POLY(X(I),Y(I)) 00003720 XSS=(XS-XSI(I))*RH 00003730 ETS=(ET-ETA(I))*RH 00003740 XS=XS*RH 00003750 ET=ET*RH 00003760 WRITE(IPP,202) NUM(I),GEOID,XS,XSS,ET,ETS 00003770 I=I+1 00003780 IF(I.EQ.495.AND.IDS.EQ.IRC) GOTO4 IF(IDS.EQ.IRC) GOTO 1041 GOTO3 00003790 4 CONTINUE 00003800 58 CONTINUE C 00003820 C COMPUTE AND PRINT VARIANCE FACTOR AND STANDARD ERROR 00003830 C 00003840 CALL CONVAR(NP,SCALE,VRA) 00003850 VARF=VARF+VRA 00003860 VARF=VARF/(2.D0*MKK-IMM-INN+NP-(NN+1)**2) 00003870 VRA=VARF*RH2 00003880 WRITE(IPP,206) VRA 00003890 C 00003900 C VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 00003910 C 00003920 DO 73 I=1,NKA 00003930 DO 72 J=1,NKA 00003940 A(I,J)=VARF*A(I,J) 00003950 FL(J)=A(I,J) 00003960 72 CONTINUE 00003970 IF(ICD.EQ.0) GOTO76 00003980 WRITE(21) FL 00003990 76 CONTINUE 00004000 73 CONTINUE 00004010 C 00004020 C CALCULATION OF GRID VALUES 00004030 C 00004040 WRITE(IPP,215) 00004050 WRITE(IPP,205) 00004060 DO 80 I=1,IQP 00004070 QX=(QP-OLAT)/RHO*R 00004080 IF(QX.EQ.0.D0) QX=1.D0 00004090 QX=QX/SCALE 00004100 QLY=QL 00004110 DO 81 J=1,IQL 00004120 QY=(QLY-OLONG)/RHO*R*DCOS(QP/RHO) 00004130 IF(QY.EQ.0.D0) QY=1.D0 00004140 QY=QY/SCALE 00004150 XS=PXSI(QX,QY) 00004160 ET=PETA(QX,QY) 00004170 GEOID=SCALE*POLY(QX,QY) 00004180 CALL SGE(QX,QY,NK,SGI) 00004190 SGI=DSQRT(SGI) 00004200 XS=XS*RH 00004210 ET=ET*RH 00004220 C GEOCENTRIC SOLUTION C CSF=DCOS(QP/RHO) SNF=DSIN(QP/RHO) CSL=DCOS(QLY/RHO) SNL=DSIN(QLY/RHO) GEOID=GEOID+(-CSF*CSL*15.-CSF*SNL*(-165.)-SNF*(-175.)-(-71.4)+6378 *135.0*(-3.715159D-05)*SNF**2) XS=XS-(SNF*CSL*15.+SNF*SNL*(-165.)-CSF*(-175.)+2*6378135.0*SNF*CSF **(-3.715159D-05))*3.24D05/(DARSIN(1.D0)*6378135.0) ET=ET-(SNL*15.-CSL*(-165.))*3.24D05/(DARSIN(1.D0)*6378135.0) WRITE(IPP,204) QP,QLY,GEOID,SGI,XS,ET 00004230 IF (ICODE.EQ.0) GOTO82 00004240 WRITE(7,105) QP,QLY,GEOID 00004250 82 CONTINUE 00004260 QLY=QLY+DEGT 00004270 81 CONTINUE 00004280 QP=QP+DEGT 00004290 80 CONTINUE 00004300 STOP 00004310 END 00004320 SUBROUTINE INDEX(NK,K,M,N,I,J,IR,IS) 00004330 NKK=NK+1 00004340 NAK=NKK*NKK-1 00004350 AM=NAK+0.5000-SQRT((NAK+0.5000)*(NAK+0.5000)-2.000*K)-1.0E-06 00004360 M=INT(AM)+1 00004370 N=NAK-NAK*M+(M*M-M)/2+K 00004380 IQS=NKK*NKK*NKK 00004390 KS=M*NKK*NKK+N 00004400 IS=KS/IQS 00004410 IQR=NKK*NKK 00004420 KR=KS-IS*IQS 00004430 IR=KR/IQR 00004440 KI=N*NKK*NKK+M 00004450 I=KI/IQS 00004460 KJ=KI-I*IQS 00004470 J=KJ/IQR 00004480 RETURN 00004490 END 00004500 FUNCTION POLY(XX,YY) 00004510 C 00004520 C CALCULATES GEOIDAL HEIGHT. 00004530 C 00004540 IMPLICIT REAL*8 (A-H,O-Z) 00004550 DIMENSION FAL(100) 00004560 COMMON /BLOCK 1/ FAL,NK 00004570 POLY=0. 00004580 N=1 00004590 I=0 00004600 DO 10 IX=1,NK 00004610 J=0 00004620 DO 11 IY=1,NK 00004630 POLY=POLY+FAL(N)*XX**I*YY**J 00004640 N=N+1 00004650 J=J+1 00004660 11 CONTINUE 00004670 I=I+1 00004680 10 CONTINUE 00004690 RETURN 00004700 END 00004710 FUNCTION PXSI(X,Y) 00004720 C 00004730 C CALCULATES DEFLECTION IN MERIDIAN. 00004740 C 00004750 IMPLICIT REAL*8 (A-H,O-Z) 00004760 DIMENSION FAL(100) 00004770 COMMON /BLOCK 1/ FAL,NK 00004780 PXSI=0. 00004790 N=1 00004800 I=0 00004810 DO 20 IX=1,NK 00004820 J=0 00004830 CALL ERRSET(208,300,-1,1) 00004840 DO 21 IY=1,NK 00004850 IF(I.EQ.0) GOTO23 00004860 PXSI=PXSI-FAL(N)*I*X**(I-1)*Y**J 00004870 23 CONTINUE 00004880 N=N+1 00004890 J=J+1 00004900 21 CONTINUE 00004910 I=I+1 00004920 20 CONTINUE 00004930 RETURN 00004940 END 00004950 FUNCTION PETA(X,Y) 00004960 C 00004970 C CALCULATES DEFLECTION IN PRIME VERTICAL. 00004980 C 00004990 IMPLICIT REAL*8 (A-H,O-Z) 00005000 DIMENSION FAL(100) 00005010 COMMON /BLOCK 1/ FAL,NK 00005020 PETA=0. 00005030 N=1 00005040 I=0 00005050 CALL ERRSET(208,300,-1,1) 00005060 DO 20 IX=1,NK 00005070 J=0 00005080 DO 21 IY=1,NK 00005090 IF(J.EQ.0) GOTO23 00005100 PETA=PETA-FAL(N)*X**I*J*Y**(J-1) 00005110 23 CONTINUE 00005120 N=N+1 00005130 J=J+1 00005140 21 CONTINUE 00005150 I=I+1 00005160 20 CONTINUE 00005170 RETURN 00005180 END 00005190 SUBROUTINE SGE(X,Y,NK,SUM) 00005200 C 00005210 C CALCULATES VARIANCE OF GEOIDAL HEIGHT. 00005220 C 00005230 IMPLICIT REAL *8 (A-H,O-Z) 00005240 DIMENSION A(100,100),B(100),BT(100) 00005250 COMMON /BLOCK 2/ A,NAK 00005260 COMMON /BLOCK 3/ SCALE 00005270 COMMON /BLOCK 4/ B,BT 00005280 I=1 00005290 IM=0 00005300 DO 10 IX=1,NK 00005310 IN=0 00005320 DO 11 IY=1,NK 00005330 B(I)=0.D0 00005340 B(I)=X**IM*Y**IN 00005350 I=I+1 00005360 IN=IN+1 00005370 11 CONTINUE 00005380 IM=IM+1 00005390 10 CONTINUE 00005400 DO 20 I=1,NAK 00005410 SUM=0.D0 00005420 DO 21 K=1,NAK 00005430 21 SUM=SUM+A(I,K)*B(K) 00005440 20 BT(I)=SUM 00005450 SUM=0.D0 00005460 DO 22 J=1,NAK 00005470 22 SUM=SUM+B(J)*BT(J) 00005480 SUM=SUM*SCALE**2 00005490 RETURN 00005500 END 00005510 SUBROUTINE SGX(X,Y,NK,SUM) 00005520 C 00005530 C CALCULATES VARIANCE OF DEFLECTION IN MERIDIAN. 00005540 C 00005550 IMPLICIT REAL*8 (A-H,O-Z) 00005560 DIMENSION A(100,100),B(100),BT(100) 00005570 COMMON /BLOCK2/ A,NAK 00005580 COMMON /BLOCK3/ SCALE 00005590 COMMON /BLOCK 4/ B,BT 00005600 I=1 00005610 IM=0 00005620 DO 10 IX=1,NK 00005630 IN=0 00005640 DO 11 IY=1,NK 00005650 B(I)=0.D0 00005660 IF(IM.EQ.0) GOTO13 00005670 B(I)=IM*X**(IM-1)*Y**IN 00005680 13 CONTINUE 00005690 I=I+1 00005700 IN=IN+1 00005710 11 CONTINUE 00005720 IM=IM+1 00005730 10 CONTINUE 00005740 DO 20 I=1,NAK 00005750 SUM=0.D0 00005760 DO 21 K=1,NAK 00005770 21 SUM=SUM+A(I,K)*B(K) 00005780 20 BT(I)=SUM 00005790 SUM=0.D0 00005800 DO 22 J=1,NAK 00005810 22 SUM=SUM+B(J)*BT(J) 00005820 RETURN 00005830 END 00005840 SUBROUTINE SGT(X,Y ,NK,SUM) 00005850 C 00005860 C CALCULATES VARIANCE OF DEFLECTION IN PRIME VERTICAL. 00005870 C 00005880 IMPLICIT REAL*8 (A-H,O-Z) 00005890 DIMENSION A(100,100),B(100),BT(100) 00005900 COMMON /BLOCK2/ A,NAK 00005910 COMMON /BLOCK3/ SCALE 00005920 COMMON /BLOCK 4/ B,BT 00005930 I=1 00005940 IM=0 00005950 DO 10 IX=1,NK 00005960 IN=0 00005970 DO 11 IY=1,NK 00005980 B(I)=0.D0 00005990 IF(IN.EQ.0) GOTO13 00006000 B(I)=X**IM*IN*Y**(IN-1) 00006010 13 CONTINUE 00006020 I=I+1 00006030 IN=IN+1 00006040 11 CONTINUE 00006050 IM=IM+1 00006060 10 CONTINUE 00006070 DO 20 I=1,NAK 00006080 SUM=0.D0 00006090 DO 21 K=1,NAK 00006100 21 SUM=SUM+A(I,K)*B(K) 00006110 20 BT(I)=SUM 00006120 SUM=0.D0 00006130 DO 22 J=1,NAK 00006140 22 SUM=SUM+B(J)*BT(J) 00006150 RETURN 00006160 END 00006170 SUBROUTINE MATRIT(A,M,N,NZ) C THIS SUBROUTINE PRINTS DOUBLE PRECISION MATRICES C A IS THE MATRIX TO BR PRINTED C M IS ROW DIMENSION OF A C N IS COLUMN DIMENSION OF A C NZ IS THE DICLARED DIMENSION OF A IN MAIN C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(NZ,N) NC=1 NEW=0 6 NEW=NEW+10 IF(N.LT.NEW) NEW=N DO 4 I=1,M 4 READ(5,300) (A(I,K),K=NC,NEW) IF(N.EQ.NEW) RETURN NC=NC+10 GOTO 6 300 FORMAT(10F8.3) END SUBROUTINE CON(AHT,NK,SCALE,OLONG,OLAT,ATPL,NP) 00006180 IMPLICIT REAL*8(A-H,O-Z) 00006190 DIMENSION P(64,64),X(64),Y(64),AL(64),AA(64,64) DIMENSION PP(64,64),ALL(64) DIMENSION A(100,100),ATPL(100) 00006220 DIMENSION AT(100,64) COMMON /BLOCK2/ A,NAK 00006240 COMMON /BLOCK4/ P,X,Y,AL 00006250 100 FORMAT(I4) 00006260 102 FORMAT(3F8.3) 00006280 200 FORMAT(//,20X,'HEIGHT CONSTRAINTS',//,30X,'LAT.',8X,'LONG.',7X,'HE00006290 *IGHT',5X,'STD. DEV.',//) 00006300 201 FORMAT(29X,F5.2,7X,F6.2,6X,F6.2,6X,F6.2) 00006310 CALL ERRSET(208,300,-1,1) 00006320 RHO=206264.80625D0/3600.D0 00006330 RH2=206264.80625D0**2 00006340 R=6370990.D0 00006350 SCC=SCALE**2 00006360 IRC=5 00006370 IPP=6 00006380 WRITE(IPP,200) 00006390 READ(IRC,100) NP 00006400 DO 10 I=1,NP DO 10 J=1,NP AA(I,J)=0.D0 10 P(I,J)=0.D0 CALL MATRIT(P,64,64,64) DO 11 I=1,NP 00006440 READ(IRC,102) X(I),Y(I),AL(I) 00006450 STD=DSQRT(P(I,I)) 00006460 WRITE(IPP,201) X(I),Y(I),AL(I),STD 00006470 Y(I)=(Y(I)-OLONG)/RHO*R*DCOS(X(I)/RHO) 00006480 X(I)=(X(I)-OLAT)/RHO*R 00006490 IF((DABS(Y(I))).LT.1.0D-50) Y(I)=Y(I)+1.0D-40 00006500 IF((DABS(X(I))).LT.1.0D-50) X(I)=X(I)+1.0D-40 00006510 X(I)=X(I)/SCALE 00006520 Y(I)=Y(I)/SCALE 00006530 AL(I)=AL(I)/SCALE 00006540 11 CONTINUE 00006550 DO 15 I=1,NP 00006560 15 AA(I,I)=P(I,I)*RH2/SCC C FORM L FROM P 00006590 CALL CHOLD(AA,64,NP,2) DO 28 I=1,NP 00006610 ALL(I)=AL(I) 00006620 DO 28 J=1,NP 00006630 PP(I,J)=AA(I,J) 28 CONTINUE 00006650 CALL CHOLD(PP,64,NP,1) C FORM AT 00006670 DO 12 K=1,NP 00006680 N=1 00006690 I=0 00006700 DO 13 IX=1,NK 00006710 J=0 00006720 DO 14 IY=1,NK 00006730 AT(N,K)=X(K)**I*Y(K)**J 00006740 N=N+1 00006750 J=J+1 00006760 14 CONTINUE 00006770 I=I+1 00006780 13 CONTINUE 00006790 12 CONTINUE 00006800 C FORM ATL 00006810 DO 16 I=1,NAK 00006820 DO 18 K=1,NP 00006830 SUM=0.D0 00006840 DO 17 J=1,NP 00006850 IF(K.GT.J) GOTO19 00006860 SUM=SUM+AT(I,J)*PP(J,K) 00006870 19 CONTINUE 00006880 17 CONTINUE 00006890 AT(I,K)=SUM 00006900 18 CONTINUE 00006910 16 CONTINUE 00006920 C FORM ATLLTA (=ATPA) 00006930 DO 20 I=1,NAK 00006940 DO 21 K=1,NAK 00006950 SUM=0.D0 00006960 DO 22 J=1,NP 00006970 SUM=SUM+AT(I,J)*AT(K,J) 00006980 22 CONTINUE 00006990 A(I,K)=SUM 00007000 21 CONTINUE 00007010 20 CONTINUE 00007020 C FORM LTW 00007030 DO 24 K=1,NP 00007040 SUM=0.D0 00007050 DO 25 J=1,NP 00007060 IF(K.GT.J) GOTO26 00007070 SUM=SUM+PP(J,K)*ALL(J) 00007080 26 CONTINUE 00007090 25 CONTINUE 00007100 ALL(K)=SUM 00007110 24 CONTINUE 00007120 C FORM ATLLTW (=ATPW) 00007130 DO 23 K=1,NAK 00007140 SUM=0.D0 00007150 DO 27 J=1,NP 00007160 SUM=SUM+AT(K,J)*ALL(J) 00007170 27 CONTINUE 00007180 ATPL(K)=SUM 00007190 23 CONTINUE 00007200 RETURN 00007210 END 00007220 SUBROUTINE CHOLD(A,IRDA,NA,ICH) 00007230 C***********************************************************************00007240 C* *00007250 C* MATRIX INVERSION SUBROUTINE USING CHOLESKI DECOMPOSITION *00007260 C* *00007270 C* INPUTS A = ARRAY CONTAINING POSITIVE DEFINITE INPUT MATRIX *00007280 C* IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX *00007290 C* NA = SIZE OF INPUT MATRIX *00007300 C* *00007310 C* OUTPUTS A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED)*00007320 C* DETA = DETERMINANT OF INPUT MATRIX *00007330 C* *00007340 C* D. WELLS UNB NOVEMBER 1972 *00007350 C* *00007360 C***********************************************************************00007370 IMPLICIT REAL*8(A-H,O-Z) 00007380 DIMENSION A(IRDA,NA) 00007390 C 00007400 C ...CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 00007410 C 00007420 DETA = A(1,1) 00007430 A(1,1) = DSQRT(A(1,1)) 00007440 IF(NA.EQ.1) GOTO6 00007450 DO 100 I = 2,NA 00007460 100 A(I,1) = A(I,1) / A(1,1) 00007470 DO 500 J = 2,NA 00007480 SUM = 0. 00007490 J1 = J - 1 00007500 DO 200 K = 1,J1 00007510 200 SUM = SUM + A(J,K) * A(J,K) 00007520 DETA = DETA * (A(J,J) - SUM) 00007530 A(J,J) = DSQRT(A(J,J) - SUM) 00007540 IF(J .EQ. NA) GOTO 500 00007550 J2 = J + 1 00007560 DO 400 I = J2,NA 00007570 SUM = 0. 00007580 DO 300 K = 1,J1 00007590 300 SUM = SUM + A(I,K) * A(J,K) 00007600 400 A(I,J) = (A(I,J) - SUM) / A(J,J) 00007610 500 CONTINUE 00007620 6 CONTINUE 00007630 IF(ICH.EQ.1) RETURN 00007640 C 00007650 C ...INVERSION OF LOWER TRIANGULAR MATRIX 00007660 C 00007670 DO 600 I = 1,NA 00007680 600 A(I,I) = 1D0 / A(I,I) 00007690 IF(NA.EQ.1) GOTO10 00007700 N1 = NA - 1 00007710 DO 800 J = 1,N1 00007720 J2 = J + 1 00007730 DO 800 I = J2,NA 00007740 SUM = 0. 00007750 I1 = I - 1 00007760 DO 700 K = J,I1 00007770 700 SUM = SUM + A(I,K) * A(K,J) 00007780 800 A(I,J) = -A(I,I) * SUM 00007790 C 00007800 C ...CONSTRUCTION OF INVERSE OF INPUT MATRIX 00007810 C 00007820 10 DO 1300 J=1,NA 00007830 IF(J .EQ. 1) GOTO 1000 00007840 J1 = J - 1 00007850 DO 900 I = 1,J1 00007860 900 A(I,J) = A(J,I) 00007870 1000 DO 1200 I = J,NA 00007880 SUM = 0. 00007890 DO 1100 K = I,NA 00007900 1100 SUM = SUM + A(K,I) * A(K,J) 00007910 1200 A(I,J) = SUM 00007920 1300 CONTINUE 00007930 RETURN 00007940 END 00007950 SUBROUTINE CONVAR(NP,SCALE,VRA) 00007960 C 00007970 C TO COMPUTE THE CONTRIBUTION TO THE VARIANCE FACTOR FROM THE HEIG00007980 C INFORMATION (I.E. VTPV) 00007990 C 00008000 IMPLICIT REAL *8(A-H,O-Z) 00008010 DIMENSION P(64,64),X(64),Y(64),AL(64),VT(64) COMMON /BLOCK4/ P,X,Y,AL 00008030 DO 1 J=1,NP 00008040 1 AL(J)=AL(J)-POLY(X(J),Y(J)) 00008050 DO 2 I=1,NP 00008060 SUM=0.D0 00008070 DO 3 J=1,NP 00008080 3 SUM=SUM+P(I,J)*AL(J) 00008090 2 VT(I)=SUM 00008100 VRA=0.D0 00008110 DO 4 I=1,NP 00008120 4 VRA=VRA+VT(I)*AL(I) 00008130 RETURN 00008140 END 00008150