C 00008210 C 00008220 C 00008230 C ANGEOID 00008240 C PROGRAMME TO COMPUTE GEOIDAL HEIGHTS FROM DEFLECTIONS OF THE VER00008250 C USING A LEAST-SQUARES SURFACE-FITTING TECHNIQUE 00008260 C 00008270 C C. L. MERRY JUNE 1972 00008280 C 00008290 C . . . MODIFIED MARCH 1973 00008300 C 00008310 C ..... VERSION TO USE TWO OR MORE DATA SETS. ---- MAY, 1973 00008320 C 00008330 C 00008340 C INPUT REQUIREMENTS ******* 00008350 C 00008360 C CARD ONE 00008370 C HEADER CARD -- CONTAINS PROJECT TITLE -- 20A4 00008380 C 00008390 C CARD TWO 00008400 C COORDINATES ( LAT. AND LONG. IN DEGREES, GEOIDAL HEIGHT IN METRES00008410 C OF DATUM ORIGIN --- 2F11.6,F6.2 00008420 C 00008430 C CARD THREE 00008440 C LIMIT OF AREA OF INTEREST ( DX,DY,MEASURED FROM ORIGIN IN KMS.) 00008450 C AND SCALE FACTOR --- 2F12.2 00008460 C 00008470 C CARD FOUR 00008480 C DEGREE OF POLYNOMIAL TO BE USED --- I2 00008490 C 00008500 C CARD FIVE 00008510 C LIMITS FOR GRID VALUES - TWO LAT. VALUES, TWO LONG. VALUES, GRID00008520 C INTERVAL, ALL IN DEGREES ---- 5F8.2 00008530 C 00008540 C CARD SIX 00008550 C COLS. 1-2 ANY INTEGER IF GRID VALUES TO BE PUNCHED. BLANK IF N00008560 C BE PUNCHED. COLS. 3-4 ANY INTEGER IF COEFFICIENTS AND 00008570 C VARIANCE-COVARIANCE MATRIX TO BE STORED ON DISC. OTHERWISE - BLAN00008580 C 00008590 C CARD SEVEN 00008600 C NUMBER OF DATA SETS TO BE USED --- I2 00008610 C 00008620 C CARD EIGHT ON 00008630 C DATA SET REFERENCE NOS. (MAGNETIC TAPE OR DISC FILE) - ONE PER 00008640 C CARD ---- I2 00008650 C 00008660 C 00008670 IMPLICIT REAL*8 (A-H,O-Z) 00008680 DIMENSION NUM(300),X(300),Y(300),WX(300),WE(300),XSI(300),ETA(300)00008690 DIMENSION A(168,168),FL(168),FAL(168) 00008700 DIMENSION IDSS(10) 00008710 DIMENSION NAME(20) 00008720 COMMON /BLOCK 1/ FAL,NK 00008730 COMMON /BLOCK 2/ A,NAK 00008740 COMMON /BLOCK 3/ SCALE 00008750 DATA FL/168*0.D0/ 00008760 101 FORMAT(2F12.2) 00008770 102 FORMAT(2F11.6,F6.2) 00008780 103 FORMAT(5F8.2) 00008790 104 FORMAT(//,20X,'NO. OF MERIDIAN DEFLECTIONS = ',I4,/,20X,'NO. OF PR00008800 *IME VERTICAL DEFLECTIONS = ',I4,//) 00008810 105 FORMAT(3F8.3) 00008820 106 FORMAT(2I2) 00008830 107 FORMAT(I2) 00008840 108 FORMAT(20A4) 00008850 200 FORMAT(/,30X,D16.8) 00008860 201 FORMAT('1',30X,'COEFFICIENTS OF',/,20X,'BEST-FITTING POLYNOMIAL',/00008870 *//) 00008880 202 FORMAT(/,10X,I8,F15.2,2(10X,F12.2,5X,F10.2)) 00008890 203 FORMAT('1',9X,'STATION',7X,'GEOIDAL',15X,'MERIDIAN',26X,'PRIME VER00008900 *TICAL',/,10X,'NUMBER',8X,'HEIGHT',14X,'DEFLECTION',6X,'RESIDUAL',100008910 *3X,'DEFLECTION',6X,'RESIDUAL',///) 00008920 204 FORMAT(2(10X,F7.2),2F15.2,2F20.2,/) 00008930 205 FORMAT(12X,'LAT.',12X,'LONG.',9X,'HEIGHT',9X,'STD. DEV.',14X,'XSI'00008940 *,17X,'ETA',//) 00008950 206 FORMAT(//,20X,'ESTIMATED VARIANCE FACTOR=',E15.6,//) 00008960 207 FORMAT(//,20X,'DEGREE OF APPROXIMATING POLYNOMIAL=',2X,I4,//) 00008970 208 FORMAT('0',30X,10X,'DATUM ORIGIN',//,10X,'LATITUDE= ',F11.6,2X,'L00008980 *ONGITUDE=',F13.6,3X,'GEOIDAL HEIGHT=',F8.2,//,20X,'LIMIT OF AREA C00008990 *OVERED - DX, DY =',F14.2,///) 00009000 209 FORMAT(30X,'LIMITING PRECISION :',I8,2X,'SIGNIFICANT FIGURES',//) 00009010 212 FORMAT(/,20X,'SCALE FACTOR =',F20.2,/) 00009020 213 FORMAT(/,10X,'TIME = ',I6) 00009030 214 FORMAT('1',35X,20A4,//) 00009040 215 FORMAT(///,40X,'GRID VALUES',//) 00009050 C 00009060 C CONSTANTS 00009070 C 00009080 IRC=5 00009090 IPP=6 00009100 R=6370990. 00009110 RH=206264.80625 00009120 RHO=206264.80625/3600. 00009130 MM=0 00009140 MKK=0 00009150 IMJ=0 00009160 INJ=0 00009170 DO 310 I=1,168 00009180 DO 310 J=1,168 00009190 A(I,J)=0.D0 00009200 310 CONTINUE 00009210 C 00009220 C READ HEADER CARD 00009230 C 00009240 READ(IRC,108) NAME 00009250 WRITE(IPP,214) NAME 00009260 C 00009270 C READ IN COORDINATES AND HEIGHT OF ORIGIN. 00009280 C 00009290 READ(IRC,102) OLAT,OLONG,AHT 00009300 C 00009310 C READ IN LIMIT OF AREA COVERED (IN KMS.) AND SCALE FACTOR 00009320 C 00009330 READ(IRC,101) ALIMIT,SCALE 00009340 WRITE(IPP,208) OLAT,OLONG,AHT,ALIMIT 00009350 ALIMIT=ALIMIT*1000.D0 00009360 WRITE(IPP,212) SCALE 00009370 C 00009380 C READ MAXIMUM DEGRE TO BE USED 00009390 C 00009400 READ(IRC,107) NAN 00009410 NN=NAN 00009420 NK=NN+1 00009430 WRITE(IPP,207) NN 00009440 NAK=NK*NK-1 00009450 C 00009460 C READ LIMITS OF GRID EVALUATION AND GRID INTERVAL 00009470 C 00009480 READ(IRC,103) QP,QPP,QL,QLL,DEGT 00009490 IQP=IABS(IDINT((QPP-QP)/DEGT))+1 00009500 IQL=IABS(IDINT((QLL-QL)/DEGT))+1 00009510 C 00009520 C READ IN CODE FOR PUNCHING OF GRID VALUES 00009530 C 00009540 READ(IRC,106) ICODE,ICD 00009550 C 00009560 C READ IN NUMBER OF DATA SETS 00009570 C 00009580 READ(IRC,107) MJ 00009590 DO 57 MI=1,MJ 00009600 C 00009610 C READ IN DATA SET REFERENCE NUMBERS 00009620 C 00009630 READ(IRC,107) IDS 00009640 IDSS(MI)=IDS 00009650 REWIND IDS 00009660 IMM=0 00009670 INN=0 00009680 I=1 00009690 C 00009700 C READ IN DATA AND TRANSFORM TO PLANE SYSTEM. 00009710 C 00009720 1 READ(IDS,END=2) NUM(I),X(I),Y(I),XSI(I),ETA(I),WX(I),WE(I) 00009730 IF(WX(I).LT.0.001) IMM=IMM+1 00009740 IF(WE(I).LT.0.001) INN=INN+1 00009750 IF(WX(I).LT.0.001) WX(I)=1.0D+20 00009760 IF(WE(I).LT.0.001) WE(I)=1.0D+20 00009770 WX(I)=1.D0/(WX(I)**2) 00009780 WE(I)=1.D0/(WE(I)**2) 00009790 IF(DABS(XSI(I)).GT.50.0) XSI(I)=0.D0 00009800 IF(DABS(ETA(I)).GT.50.0) ETA(I)=0.D0 00009810 Y(I)=(Y(I)-OLONG)/RHO*R*DCOS(X(I)/RHO) 00009820 X(I)=(X(I)-OLAT)/RHO*R 00009830 IF(DABS(X(I)).GT.ALIMIT) GOTO1 00009840 IF(DABS(Y(I)).GT.ALIMIT) GOTO1 00009850 IF(X(I).EQ.0.D0) X(I)=1.D0 00009860 IF(Y(I).EQ.0.D0) Y(I)=1.D0 00009870 X(I)=X(I)/SCALE 00009880 Y(I)=Y(I)/SCALE 00009890 XSI(I)=XSI(I)/RH 00009900 ETA(I)=ETA(I)/RH 00009910 I=I+1 00009920 GOTO1 00009930 2 CONTINUE 00009940 MM=I-1 00009950 IF(MM.EQ.0) GOTO57 00009960 C 00009970 C FORMATION OF MATRIX OF NORMAL EQUATIONS 00009980 C 00009990 CALL ERRSET(208,300,-1,1) 00010000 CALL CPUTIM(ITIM,IREM) 00010010 IT=ITIM/10000 00010020 WRITE(IPP,213) IT 00010030 NPP=(NAK**2+NAK)/2 00010040 DO 62 K=1,NPP 00010050 CALL INDEX(NN,K,M,N,I,J,IR,IS) 00010060 AMN=0.D0 00010070 BMN=0.D0 00010080 ISS=IS+I-2 00010090 IF(ISS.LT.0) GOTO66 00010100 IF(IS.EQ.0) GOTO66 00010110 DO 64 II=1,MM 00010120 AMN=AMN+WX(II)*X(II)**ISS*Y(II)**(IR+J) 00010130 64 CONTINUE 00010140 66 CONTINUE 00010150 ISR=IR+J-2 00010160 IF(ISR.LT.0) GOTO69 00010170 IF(IR.EQ.0) GOTO69 00010180 DO 68 II=1,MM 00010190 BMN=BMN+WE(II)*X(II)**(IS+I)*Y(II)**ISR 00010200 68 CONTINUE 00010210 69 CONTINUE 00010220 A(M,N)=A(M,N)+I*IS*AMN+J*IR*BMN 00010230 62 CONTINUE 00010240 C 00010250 C FORMATION OF OBSERVATION VECTOR 00010260 C 00010270 N=1 00010280 IS=0 00010290 DO 56 IF=1,NK 00010300 IR=0 00010310 DO 55 IE=1,NK 00010320 IF(IS.EQ.0.AND.IR.EQ.0) GOTO54 00010330 ISK=IS-1 00010340 IF(ISK.LT.0) GOTO52 00010350 FLN=0.D0 00010360 DO 50 II=1,MM 00010370 FLN=FLN+XSI(II)*WX(II)*X(II)**ISK*Y(II)**IR 00010380 50 CONTINUE 00010390 52 CONTINUE 00010400 IRK=IR-1 00010410 IF(IRK.LT.0) GOTO53 00010420 FALN=0.D0 00010430 DO 51 II=1,MM 00010440 FALN=FALN+ETA(II)*WE(II)*X(II)**IS*Y(II)**IRK 00010450 51 CONTINUE 00010460 53 CONTINUE 00010470 FL(N)=FL(N)-IS*FLN-IR*FALN 00010480 N=N+1 00010490 54 CONTINUE 00010500 IR=IR+1 00010510 55 CONTINUE 00010520 IS=IS+1 00010530 56 CONTINUE 00010540 CALL CPUTIM(ITIM,IREM) 00010550 IT=ITIM/10000 00010560 WRITE(IPP,213) IT 00010570 MKK=MKK+MM 00010580 IMJ=IMJ+MM-IMM 00010590 INJ=INJ+MM-INN 00010600 WRITE(IPP,104) IMJ,INJ 00010610 57 CONTINUE 00010620 DO 471 I=1,NAK 00010630 DO 471 J=1,NAK 00010640 IF(J.GE.I) GOTO472 00010650 A(I,J)=A(J,I) 00010660 472 CONTINUE 00010670 471 CONTINUE 00010680 TEST1=A(NAK,NAK) 00010690 DO 800 M=1,NAK 00010700 800 A(M,M)=A(M,M)+4.D0*NAK*0.54D-78 00010710 C 00010720 C COMPUTE AND PRINT VECTOR OF COEFFICIENTS. 00010730 C 00010740 C 00010750 C ...CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX 00010760 C 00010770 NA=NAK 00010780 DETA = A(1,1) 00010790 A(1,1) = DSQRT(A(1,1)) 00010800 DO 910 I=2,NA 00010810 910 A(I,1)=A(I,1)/A(1,1) 00010820 DO 500 J = 2,NA 00010830 SUM = 0. 00010840 J1 = J - 1 00010850 DO 920 K=1,J1 00010860 920 SUM=SUM+A(J,K)*A(J,K) 00010870 DETA = DETA * (A(J,J) - SUM) 00010880 A(J,J) = DSQRT(A(J,J) - SUM) 00010890 IF(J .EQ. NA) GOTO 500 00010900 J2 = J + 1 00010910 DO 400 I = J2,NA 00010920 SUM = 0. 00010930 DO 300 K = 1,J1 00010940 300 SUM = SUM + A(I,K) * A(J,K) 00010950 400 A(I,J) = (A(I,J) - SUM) / A(J,J) 00010960 500 CONTINUE 00010970 C 00010980 C ...INVERSION OF LOWER TRIANGULAR MATRIX 00010990 C 00011000 DO 600 I = 1,NA 00011010 600 A(I,I) = 1D0 / A(I,I) 00011020 N1 = NA - 1 00011030 DO 980 J=1,N1 00011040 J2 = J + 1 00011050 DO 980 I=J2,NA 00011060 SUM = 0. 00011070 I1 = I - 1 00011080 DO 700 K = J,I1 00011090 700 SUM = SUM + A(I,K) * A(K,J) 00011100 980 A(I,J)=-A(I,I)*SUM 00011110 C 00011120 C ...CONSTRUCTION OF INVERSE OF INPUT MATRIX 00011130 C 00011140 DO 1300 J = 1,NA 00011150 IF(J .EQ. 1) GOTO 1000 00011160 J1 = J - 1 00011170 DO 900 I = 1,J1 00011180 900 A(I,J) = A(J,I) 00011190 1000 DO 1200 I = J,NA 00011200 SUM = 0. 00011210 DO 1100 K = I,NA 00011220 1100 SUM = SUM + A(K,I) * A(K,J) 00011230 1200 A(I,J) = SUM 00011240 1300 CONTINUE 00011250 CALL CPUTIM(ITIM,IREM) 00011260 IT=ITIM/10000 00011270 WRITE(IPP,213) IT 00011280 WRITE(6,805) DETA 00011290 805 FORMAT('1',30X,'DETERMINANT= ',D15.6,//) 00011300 TEST2=A(NAK,NAK) 00011310 TEA=TEST1*TEST2 00011320 ITEA=IDINT(16.8D0-DLOG10(TEA)) 00011330 WRITE(IPP,209) ITEA 00011340 DO 510 I=1,NAK 00011350 SUM=0.D0 00011360 DO 511 K=1,NAK 00011370 511 SUM=SUM+A(I,K)*FL(K) 00011380 510 FAL(I)=SUM 00011390 WRITE(IPP,201) 00011400 DO 70 I=1,NAK 00011410 WRITE(IPP,200) FAL(I) 00011420 70 CONTINUE 00011430 IF(ICD.EQ.0) GOTO75 00011440 WRITE(20) FAL 00011450 75 CONTINUE 00011460 C 00011470 C COMPUTE AND PRINT GEOIDAL UNDULATIONS , DEFLECTIONS, RESIDUALS. 00011480 C 00011490 WRITE(IPP,203) 00011500 VARF=0. 00011510 KIA=0 00011520 DO 58 MI=1,MJ 00011530 IDS=IDSS(MI) 00011540 REWIND IDS 00011550 I=1 00011560 3 READ(IDS,END=4) NUM(I),X(I),Y(I),XSI(I),ETA(I),WX(I),WE(I) 00011570 IF(WX(I).LT.0.001) WX(I)=1.0D+20 00011580 IF(WE(I).LT.0.001) WE(I)=1.0D+20 00011590 WX(I)=1.D0/(WX(I)**2) 00011600 WE(I)=1.D0/(WE(I)**2) 00011610 IF(DABS(XSI(I)).GT.50.0) XSI(I)=0.D0 00011620 IF(DABS(ETA(I)).GT.50.0) ETA(I)=0.D0 00011630 Y(I)=(Y(I)-OLONG)/RHO*R*DCOS(X(I)/RHO) 00011640 X(I)=(X(I)-OLAT)/RHO*R 00011650 IF(DABS(X(I)).GT.ALIMIT) GOTO3 00011660 IF(DABS(Y(I)).GT.ALIMIT) GOTO3 00011670 IF(X(I).EQ.0.D0) X(I)=1.D0 00011680 IF(Y(I).EQ.0.D0) Y(I)=1.D0 00011690 X(I)=X(I)/SCALE 00011700 Y(I)=Y(I)/SCALE 00011710 XSI(I)=XSI(I)/RH 00011720 ETA(I)=ETA(I)/RH 00011730 KIA=KIA+1 00011740 AA=KIA 00011750 LL=KIA/25 00011760 IF((AA/25-LL).LE.1E-07)WRITE(IPP,203) 00011770 XS=PXSI(X(I),Y(I)) 00011780 ET=PETA(X(I),Y(I)) 00011790 VARF=VARF+WX(I)*(XS-XSI(I))**2+WE(I)*(ET-ETA(I))**2 00011800 GEOID=AHT+SCALE*POLY(X(I),Y(I)) 00011810 XSS=(XS-XSI(I))*RH 00011820 ETS=(ET-ETA(I))*RH 00011830 XS=XS*RH 00011840 ET=ET*RH 00011850 WRITE(IPP,202) NUM(I),GEOID,XS,XSS,ET,ETS 00011860 I=I+1 00011870 GOTO3 00011880 4 CONTINUE 00011890 58 CONTINUE 00011900 C 00011910 C COMPUTE AND PRINT VARIANCE FACTOR AND STANDARD ERROR 00011920 C 00011930 VARF=VARF/(2.D0*MKK-IMM-INN-(NN+1)**2+1.D0) 00011940 VRA=VARF*RH**2 00011950 WRITE(IPP,206) VRA 00011960 C 00011970 C VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 00011980 C 00011990 DO 73 I=1,NAK 00012000 DO 72 J=1,NAK 00012010 A(I,J)=VARF*A(I,J) 00012020 FL(J)=A(I,J) 00012030 72 CONTINUE 00012040 IF(ICD.EQ.0) GOTO76 00012050 WRITE(21) FL 00012060 76 CONTINUE 00012070 73 CONTINUE 00012080 C 00012090 C CALCULATION OF GRID VALUES 00012100 C 00012110 WRITE(IPP,215) 00012120 WRITE(IPP,205) 00012130 DO 80 I=1,IQP 00012140 QX=(QP-OLAT)/RHO*R 00012150 IF(QX.EQ.0) QX=1.D0 00012160 QX=QX/SCALE 00012170 QLY=QL 00012180 DO 81 J=1,IQL 00012190 QY=(QLY-OLONG)/RHO*R*DCOS(QP/RHO) 00012200 IF(QY.EQ.0) QY=1.D0 00012210 QY=QY/SCALE 00012220 XS=PXSI(QX,QY) 00012230 ET=PETA(QX,QY) 00012240 GEOID=AHT+SCALE*POLY(QX,QY) 00012250 CALL SGE(QX,QY,NK,SGI) 00012260 SGI=DSQRT(SGI) 00012270 XS=XS*RH 00012280 ET=ET*RH 00012290 WRITE(IPP,204) QP,QLY,GEOID,SGI,XS,ET 00012300 IF (ICODE.EQ.0) GOTO82 00012310 WRITE(7,105) QP,QLY,GEOID 00012320 82 CONTINUE 00012330 QLY=QLY+DEGT 00012340 81 CONTINUE 00012350 QP=QP+DEGT 00012360 80 CONTINUE 00012370 STOP 00012380 END 00012390 SUBROUTINE INDEX(NK,K,M,N,I,J,IR,IS) 00012400 NKK=NK+1 00012410 NAK=NKK*NKK-1 00012420 AM=NAK+0.5000-SQRT((NAK+0.5000)*(NAK+0.5000)-2.000*K)-1.0E-06 00012430 M=INT(AM)+1 00012440 N=NAK-NAK*M+(M*M-M)/2+K 00012450 IQS=NKK*NKK*NKK 00012460 KS=M*NKK*NKK+N 00012470 IS=KS/IQS 00012480 IQR=NKK*NKK 00012490 KR=KS-IS*IQS 00012500 IR=KR/IQR 00012510 KI=N*NKK*NKK+M 00012520 I=KI/IQS 00012530 KJ=KI-I*IQS 00012540 J=KJ/IQR 00012550 RETURN 00012560 END 00012570 FUNCTION POLY(XX,YY) 00012580 C 00012590 C CALCULATES GEOIDAL HEIGHT. 00012600 C 00012610 IMPLICIT REAL*8 (A-H,O-Z) 00012620 DIMENSION FAL(168) 00012630 COMMON /BLOCK 1/ FAL,NK 00012640 POLY=0. 00012650 N=1 00012660 I=0 00012670 DO 10 IX=1,NK 00012680 J=0 00012690 DO 11 IY=1,NK 00012700 IF(I.EQ.0.AND.J.EQ.0) GOTO12 00012710 POLY=POLY+FAL(N)*XX**I*YY**J 00012720 N=N+1 00012730 12 CONTINUE 00012740 J=J+1 00012750 11 CONTINUE 00012760 I=I+1 00012770 10 CONTINUE 00012780 RETURN 00012790 END 00012800 FUNCTION PXSI(X,Y) 00012810 C 00012820 C CALCULATES DEFLECTION IN MERIDIAN. 00012830 C 00012840 IMPLICIT REAL*8 (A-H,O-Z) 00012850 DIMENSION FAL(168) 00012860 N=1 0 12 8Z PXSI=0. 00012880 COMMON /BLOCK 1/ FAL,NK 00012870 I=0 00012900 DO 20 IX=1,NK 00012910 J=0 00012920 CALL ERRSET(208,300,-1,1) 00012930 DO 21 IY=1,NK 00012940 IF(I.EQ.0.AND.J.EQ.0) GOTO22 00012950 IF(I.EQ.0) GOTO23 00012960 PXSI=PXSI-FAL(N)*I*X**(I-1)*Y**J 00012970 23 CONTINUE 00012980 N=N+1 00012990 22 CONTINUE 00013000 J=J+1 00013010 21 CONTINUE 00013020 I=I+1 00013030 20 CONTINUE 00013040 RETURN 00013050 END 00013060 FUNCTION PETA(X,Y) 00013070 C 00013080 C CALCULATES DEFLECTION IN PRIME VERTICAL. 00013090 C 00013100 IMPLICIT REAL*8 (A-H,O-Z) 00013110 DIMENSION FAL(168) 00013120 COMMON /BLOCK 1/ FAL,NK 00013130 PETA=0. 00013140 N=1 00013150 I=0 00013160 CALL ERRSET(208,300,-1,1) 00013170 DO 20 IX=1,NK 00013180 J=0 00013190 DO 21 IY=1,NK 00013200 IF(I.EQ.0.AND.J.EQ.0) GOTO22 00013210 IF(J.EQ.0) GOTO23 00013220 PETA=PETA-FAL(N)*X**I*J*Y**(J-1) 00013230 23 CONTINUE 00013240 N=N+1 00013250 22 CONTINUE 00013260 21 CONTINUE 00013280 I=I+1 00013290 20 CONTINUE 00013300 RETURN 00013310 END 00013320 SUBROUTINE SGE(X,Y,NK,SUM) 00013330 C 00013340 C CALCULATES VARIANCE OF GEOIDAL HEIGHT. 00013350 C 00013360 IMPLICIT REAL *8 (A-H,O-Z) 00013370 DIMENSION A(168,168),B(168),BT(168) 00013380 COMMON /BLOCK 2/ A,NAK 00013390 COMMON /BLOCK 3/ SCALE 00013400 COMMON /BLOCK 4/ B,BT 00013410 I=1 00013420 IM=0 00013430 DO 10 IX=1,NK 00013440 IN=0 00013450 DO 11 IY=1,NK 00013460 B(I)=0.D0 00013470 IF(IM.EQ.0.AND.IN.EQ.0) GOTO12 00013480 B(I)=X**IM*Y**IN 00013490 I=I+1 00013500 12 CONTINUE 00013510 IN=IN+1 00013520 11 CONTINUE 00013530 IM=IM+1 00013540 10 CONTINUE 00013550 DO 20 I=1,NAK 00013560 SUM=0.D0 00013570 DO 21 K=1,NAK 00013580 21 SUM=SUM+A(I,K)*B(K) 00013590 20 BT(I)=SUM 00013600 SUM=0.D0 00013610 DO 22 J=1,NAK 00013620 22 SUM=SUM+B(J)*BT(J) 00013630 SUM=SUM*SCALE**2 00013640 RETURN 00013650 END 00013660 SUBROUTINE SGX(X,Y,NK,SUM) 00013670 C 00013680 C CALCULATES VARIANCE OF DEFLECTION IN MERIDIAN. 00013690 C 00013700 IMPLICIT REAL*8 (A-H,O-Z) 00013710 DIMENSION A(168,168),B(168),BT(168) 00013720 COMMON /BLOCK2/ A,NAK 00013730 COMMON /BLOCK3/ SCALE 00013740 COMMON /BLOCK 4/ B,BT 00013750 I=1 00013760 IM=0 00013770 DO 10 IX=1,NK 00013780 IN=0 00013790 DO 11 IY=1,NK 00013800 B(I)=0.D0 00013810 IF(IM.EQ.0.AND.IN.EQ.0) GOTO12 00013820 IF(IM.EQ.0) GOTO13 00013830 B(I)=IM*X**(IM-1)*Y**IN 00013840 13 CONTINUE 00013850 I=I+1 00013860 12 CONTINUE 00013870 IN=IN+1 00013880 11 CONTINUE 00013890 IM=IM+1 00013900 10 CONTINUE 00013910 DO 20 I=1,NAK 00013920 SUM=0.D0 00013930 DO 21 K=1,NAK 00013940 21 SUM=SUM+A(I,K)*B(K) 00013950 20 BT(I)=SUM 00013960 SUM=0.D0 00013970 DO 22 J=1,NAK 00013980 22 SUM=SUM+B(J)*BT(J) 00013990 RETURN 00014000 END 00014010 SUBROUTINE SGT(X,Y ,NK,SUM) 00014020 C 00014030 C CALCULATES VARIANCE OF DEFLECTION IN PRIME VERTICAL. 00014040 C 00014050 IMPLICIT REAL*8 (A-H,O-Z) 00014060 DIMENSION A(168,168),B(168),BT(168) 00014070 COMMON /BLOCK2/ A,NAK 00014080 COMMON /BLOCK3/ SCALE 00014090 COMMON /BLOCK 4/ B,BT 00014100 I=1 00014110 IM=0 00014120 DO 10 IX=1,NK 00014130 IN=0 00014140 DO 11 IY=1,NK 00014150 B(I)=0.D0 00014160 IF(IM.EQ.0.AND.IN.EQ.0) GOTO12 00014170 IF(IN.EQ.0) GOTO13 00014180 B(I)=X**IM*IN*Y**(IN-1) 00014190 13 CONTINUE 00014200 I=I+1 00014210 12 CONTINUE 00014220 IN=IN+1 00014230 11 CONTINUE 00014240 IM=IM+1 00014250 10 CONTINUE 00014260 DO 20 I=1,NAK 00014270 SUM=0.D0 00014280 DO 21 K=1,NAK 00014290 21 SUM=SUM+A(I,K)*B(K) 00014300 20 BT(I)=SUM 00014310 SUM=0.D0 00014320 DO 22 J=1,NAK 00014330 22 SUM=SUM+B(J)*BT(J) 00014340 RETURN 00014350 END 00014360