C***********************************************************************00000050 C* *00000060 C* DOVEGUN2 *00000070 C* *00000080 C* COMPUTATION OF DEFLECTIONS OF THE VERTICAL (XSI & ETA) AND GEOID *00000090 C* UNDULATIONS (N) AND THEIR RESPECTIVE VARIANCE-COVARIANCE MATRICIES *00000100 C* USING THE LEAST-SQUARES SURFACE FITTING POLYNOMIAL OF THE GEOID *00000110 C* COMPUTED FROM OBSERVED DEFLECTIONS OF THE VERTICAL. *00000120 C* THE POLYNOMIAL COEFFICIENTS AND THEIR VARIANCES AND COVARIANCES *00000130 C* ARE OBTAINED FROM A PROGRAMME BY C.L.MERRY - DETERMINATION OF THE *00000140 C* GEOID FROM DEFLECTIONS OF THE VERTICAL USING A LEAST-SQUARES SURFACE*00000150 C* FITTING TECHNIQUE. *00000160 C* *00000170 C***********************************************************************00000180 C* INPUT DATA *00000190 C* *00000200 C* CARD 1 COLS 1-80 JOB DESCRIPTION (20A4) *00000210 C* *00000220 C* CARD 2 DATUM ORIGIN INFORMATION *00000230 C* COLS 1-6 STATION NUMBER (I6) *00000240 C* COLS 7-11 LATITUDE DEGS (I5) 00000250 C* COLS 12-14 MINS (I3) 00000260 C* COLS 15-22 SECS (F8.4) 00000270 C* COLS 23-27 LONGITUDE (+VE E) DEGS (I5) 00000280 C* COLS 28-30 MINS (I3) 00000290 C* COLS 31-38 SECS (F8.4) 00000300 C* COLS 39-46 GEOID UNDULATION (METRES) (F8.3) 00000310 C* COLS 47-53 XSI (F6.2) 00000320 C* COLS 54-60 ETA (F6.2) 00000330 C* *00000340 C* CARD 3 COLS 1-2 DEGREE OF THE POLYNOMIAL USED (I2) *00000350 C* COLS 4-15 SCALE FACTOR (F12.2) *00000360 C* NOTE : SEE PROGRAMME DOCUEMENTATION OF 'ANGEOID' FOR DETAILS RE THE *00000370 C* DEGREE OF THE POLYNOMIAL BEING USED HEREIN AND THE SCALE FACTOR *00000380 C* BEING USED. *00000390 C* CARD 3A COLS 1-2 PUNCH CODE (I2) *00000400 C* IF NIT BLANK, GEOIDAL HEIGHTS AND THEIR ERROR COVARIANCE *00000410 C* MATRIX ARE PUNCHED. *00000420 C* *00000430 C* CARD 4 STATION DATA ( 1 SET OF STATION INFORMATION PER CARD *00000440 C* WHERE N,XSI, AND ETA ARE TO BE EVALUATED. *00000450 C* COLS 1-6 STATION NUMBER (I6) *00000460 C* COLS 7-11 LATITUDE DEGS (I5) 00000470 C* COLS 12-14 MINS (I3) 00000480 C* COLS 15-22 SECS (F8.4) 00000490 C* COLS 23-27 LONGITUDE (+VE E) DEGS (I5) 00000500 C* COLS 28-30 MINS (I3) 00000510 C* COLS 31-38 SECS (F8.4) 00000520 C* NOTE: THE LAST CARD OF THIS TYPE MUST HAVE ONLY THE STATION *00000530 C* NUMBER 000000 PUNCHED IN COLUMNS 1-6 *00000540 C* *00000550 C* READ THE VARIANCE-COVARIANCE MATRIX FROM DISC DATA SET (RESULTS *00000560 C* FROM 'ANGEOID') *00000570 C* *00000580 C* READ THE POLYNOMIAL COEFFICIENTS FROM DISC DATA SET ( RESULTS FROM*00000590 C* 'ANGEOID') *00000600 C* *00000610 C* NOTE: SEE BEGINNING OF THE PROGRAMME FOR THE REQUIRED JCL *00000620 C***********************************************************************00000630 C* *00000640 IMPLICIT REAL*8 (A-H,O-Z) 00000650 REAL*8 X(50),Y(50),B(50,100),SIGN(50,50),SIGC(100,100),C(100) 00000660 REAL*8 SIG(50,100) 00000670 REAL*4 DJOB(20),SLATO,SLONO,SLAT,SLON,XSI(50),ETA(50),GN(50) 00000680 REAL*4 STDN(50),STDX(50),STDE(50) 00000690 DIMENSION IST(50) 00000700 C* *00000710 C* REQUIRED CONSTANTS *00000720 C* NOTE : (1) MAX MUST BE CHANGED TO AGREE WITH THE MAXIMUM ALLOWABLE *00000730 C* POINTS ( X(I),Y(I) ) *00000740 C* (2) R (MEAN EARTH RADIUS ) MUST AGREE WITH THAT USED IN THE *00000750 C* PROGRAMME USED TO GENERATE THE GEOID POLYNOMIAL *00000760 C* *00000770 MAX=50 00000780 R=6370990.D0 00000790 RHO=206264.8062471D0 00000800 PI=3.1415926536D0 00000810 RH=RHO/3600.D0 00000820 C* *00000830 C* READ DATUM INFORMATION *00000840 C* *00000850 READ(5,1) DJOB 00000860 1 FORMAT (20(A4)) 00000870 WRITE (6,40) DJOB 00000880 40 FORMAT('1',////,10X,20A4,////) 00000890 READ(5,2) ISTO,LATDO,LATMO,SLATO,LONDO,LONMO,SLONO,GNO,XSIO,ETAO 00000900 2 FORMAT(I6,2(I5,I3,F8.4),F8.3,2F7.2) 00000910 CALL ARCRAD (LATDO,LATMO,SLATO,RXO) 00000920 CALL ARCRAD (LONDO,LONMO,SLONO,RYO) 00000930 IF(RYO)110,111,111 00000940 110 WRITE(6,112) 00000950 112 FORMAT(' ',9X,'****NOTE**** THE LONGITUDE OF THE ORIGIN HAS BEEN M00000960 1ADE POSITIVE (EAST) FOR COMPUTATION PURPOSES',////) 00000970 RYO=2.D0*PI+RYO 00000980 C* *00000990 C* READ POLYNOMIAL INFORMATION *00001000 C* *00001010 111 READ(5,3) IDP,SF 00001020 3 FORMAT(I2,F13.2) 00001030 ID=(1+IDP)**2-1 00001040 I=1 00001050 C* 00001060 C* READ PUNCH CODE 00001070 C* 00001080 READ(5,200) ICODE 00001090 200 FORMAT(I2) 00001100 C* *00001110 C* WRITE JOB HEADING AND DATUM ORIGIN INFORMATION *00001120 C* *00001130 WRITE(6,41) ISTO,LATDO,LATMO,SLATO, LONDO,LONMO,SLONO,GNO,XSIO, 00001140 1ETAO 00001150 41 FORMAT(' ',9X,'DATUM ORIGIN DATA',//,10X,'STATION : ',I6,3X,'LAT.:00001160 1 ',I4,I3,F7.3,3X,'LONG.(+VE E.) : ',I4,I3,F7.3,/,10X,'GEOID HT.(M)00001170 2 : ',F8.3,3X,'DEFL.MERID.(XSI) : ',F6.3,' (SECS)',3X,'DEFL.PR.VERT00001180 3.(ETA) : ',F6.3,' (SECS)',////) 00001190 C* *00001200 C* WRITE POLYNOMIAL INFORMATION *00001210 C* *00001220 WRITE(6,42) IDP,ID,SF 00001230 42 FORMAT(' ',9X,'POLYNOMIAL DATA',//,10X,'DEGREE OF THE POLYNOMIAL :00001240 1 ',I2,/,10X,'NUMBER OF COEFFICIENTS : ',I3,/,10X,'SCALE FACTOR APP00001250 2LIED : ',F13.2,//////) 00001260 C* *00001270 C* READ STATION NUMBERS AND COORDS WHERE GEOID HTS AND DEFLECTIONS *00001280 C* ARE TO BE COMPUTED *00001290 C* *00001300 5 READ(5,4) IST(I),LATD,LATM,SLAT,LOND,LONM,SLON 00001310 4 FORMAT(I6,2(I5,I3,F8.4)) 00001320 IF(IST(I).EQ.0) GO TO 28 00001330 CALL ARCRAD (LATD,LATM,SLAT,RX) 00001340 CALL ARCRAD(LOND,LONM,SLON,RY) 00001350 IF(RY)80,81,81 00001360 80 WRITE(6,82) IST(I) 00001370 82 FORMAT(' ',9X,'****IMPORTANT MESSAGE**** THE LONGITUDE OF STATION'00001380 1,I7,' HAS BEEN CHANGED FROM A NEGATIVE VALUE TO ITS',/,36X,'EQUIVA00001390 2LENT POSITIVE (EAST) VALUE',//) 00001400 RY=2.D0*PI+RY 00001410 81 IF(DABS(RX-RXO).LE.1.D0/RHO) GO TO 60 00001420 GO TO 62 00001430 60 WRITE(6,61) IST(I) 00001440 61 FORMAT(' ',9X,'****IMPORTANT MESSAGE**** THE LATITUDE OF STATION',00001450 1I7,' HAS BEEN SHIFTED 1 ARCSEC NORTH FOR COMPUTATIONS',/,36X,'IT W00001460 2AS COINCIDENT WITH THE LATITUDE OF THE ORIGIN.',//) 00001470 RX=RX+1.D0/RHO 00001480 62 X(I)=R*(RX-RXO)/SF 00001490 IF(DABS(RY-RYO).LE.1.D0/RHO) GO TO 63 00001500 GO TO 64 00001510 63 WRITE(6,65) IST(I) 00001520 65 FORMAT(' ',9X,'****IMPORTANT MESSAGE**** THE LONGITUDE OF STATION'00001530 1,I7,' HAS BEEN SHIFTED 1 ARCSEC EAST FOR COMPUTATIONS',/,36X,'IT W00001540 2AS COINCIDENT WITH THE LONGITUDE OF THE ORIGIN',//) 00001550 RY=RY+1.D0/RHO 00001560 64 Y(I)=(R*DCOS(RX)*(RY-RYO))/SF 00001570 I=I+1 00001580 IF(I.LE.MAX) GO TO 5 00001590 READ(5,25) IDUM 00001600 25 FORMAT(I6) 00001610 IF(IDUM.EQ.0) GO TO 28 00001620 GO TO 101 00001630 28 IQ=I-1 00001640 C* *00001650 C* READ VAR-COVAR MATRIX OF THE POLYNOMIAL COEFFICIENTS FROM DISC *00001660 C* *00001670 DO 7 I=1,ID 00001680 READ(12) C 00001690 DO 8 J=1,ID 00001700 8 SIGC(I,J)=C(J) 00001710 7 CONTINUE 00001720 C* *00001730 C* READ POLYNOMIAL COEFFICIENTS FROM DISC *00001740 C* *00001750 READ(11) C 00001760 WRITE(6,27) IQ 00001770 27 FORMAT(' ',////,9X,'NUMBER OF STATIONS AT WHICH N,XSI,&ETA ARE TO 00001780 1BE COMPUTED : ',I3) 00001790 WRITE(6,45) 00001800 45 FORMAT('1',5X,'POLYNOMIAL COEFFICIENTS',//) 00001810 CALL MATRIT(C,ID,1,143) 00001820 C* *00001830 C* CALL SUBROUTINE GUN TO COMPUTE N'S AND FORM THE TRANSFORMATION *00001840 C* MATRIX REQUIRED TO COMPUTE THE VARIANCE-COVARIANCE MATRIX OF THE N'S00001850 C* *00001860 CALL GUN (IQ,SF,GNO,X,Y,C,IDP,GN,B,ID,50) 00001870 IF(ICODE.EQ.0) GOTO205 00001880 DO 206 I=1,IQ 00001890 RX=(X(I)*SF)/R+RXO 00001900 RY=(Y(I)*SF)/(R*DCOS(RX))+RYO 00001910 RX=RX*RH 00001920 RY=RY*RH 00001930 WRITE(7,202) RX,RY,GN(I) 00001940 202 FORMAT(3F8.3) 00001950 206 CONTINUE 00001960 205 CONTINUE 00001970 C* *00001980 C* COMPUTE THE VAR-COVAR MATRIX OF THE COMPUTED GEOID UNDULATIONS *00001990 C* COMPUTE THE STANDARD DEVIATIONS OF THE RESPECTIVE GEOID HEIGHTS *00002000 C* *00002010 DO 9 I=1,IQ 00002020 DO 9 J=1,ID 00002030 9 SIG(I,J)=0.D0 00002040 DO 89 I=1,IQ 00002050 DO 89 J=1,IQ 00002060 89 SIGN(I,J)=0.D0 00002070 DO 10 I=1,IQ 00002080 DO 10 J=1,ID 00002090 DO 10 K=1,ID 00002100 10 SIG(I,J)=B(I,K)*SIGC(K,J)+SIG(I,J) 00002110 DO 11 I=1,IQ 00002120 DO 11 J=1,IQ 00002130 DO 11 K=1,ID 00002140 11 SIGN(I,J)=SIG(I,K)*B(J,K)+SIGN(I,J) 00002150 DO 12 I=1,IQ 00002160 DO 12 J=1,IQ 00002170 12 SIGN(I,J)=SIGN(I,J)*SF**2 00002180 DO 79 I=1,IQ 00002190 IF(ICODE.EQ.0) GOTO208 00002200 201 FORMAT(10F8.3) 00002210 WRITE(7,201) (SIGN(I,J),J=1,10) 00002220 WRITE(7,201) (SIGN(I,J),J=11,IQ) 00002230 208 CONTINUE 00002240 79 STDN(I)=DSQRT(SIGN(I,I)) 00002250 C* *00002260 C* WRITE THE VAR-COVAR MATRIX OF THE COMPUTED GEOID UNDULATIONS *00002270 C* *00002280 WRITE(6,50) 00002290 50 FORMAT('1',10X,'VARIANCE-COVARIANCE MATRIX OF THE COMPUTED GEOID H00002300 1EIGHTS',///) 00002310 CALL MATRIT(SIGN,IQ,IQ,50) 00002320 C* *00002330 C* CALL DOVE1 TO COMPUTE THE DEFLECTION OF THE VERTICAL IN THE *00002340 C* MERIDAN PLANE, AND THE TRANSFORMATION MATRIX REQ'D TO COMPUTE *00002350 C* THE VAR-COVAR MATRIX OF THESE DEFLECTIONS. *00002360 C* *00002370 CALL DOVE1(IQ,X,Y,C,IDP,XSI,B,ID,50) 00002380 DO 16 I=1,IQ 00002390 DO 16 J= 1,ID 00002400 16 SIG( I,J)=0.D0 00002410 DO 17 I=1,IQ 00002420 DO 17 J=1,IQ 00002430 17 SIGN(I,J)=0.D0 00002440 DO 13 I=1,IQ 00002450 DO 13 J=1,ID 00002460 DO 13 K=1,ID 00002470 13 SIG(I,J)=B(I,K)*SIGC(K,J)+SIG(I,J) 00002480 DO 14 I=1,IQ 00002490 DO 14 J=1,IQ 00002500 DO 14 K=1,ID 00002510 14 SIGN(I,J)=SIG(I, K)*B(J,K)+SIGN(I,J) 00002520 DO 75 I=1,IQ 00002530 DO 75 J=1,IQ 00002540 75 SIGN(I,J)=SIGN(I,J)*RHO**2 00002550 DO 77 I=1,IQ 00002560 77 STDX(I)=DSQRT(SIGN(I,I)) 00002570 C* *00002580 C* WRITE THE VAR-COVAR MATRIX OF THE MERID COMPONENT OF THE *00002590 C* COMPUTED DEFLECTION *00002600 C* *00002610 WRITE(6,52) 00002620 52 FORMAT('1',10X,'VARIANCE-COVARIANCE MATRIX OF THE COMPUTED DEFLECT00002630 1ION IN THE MERIDIAN(XSI)',///) 00002640 CALL MATRIT(SIGN,IQ,IQ,50) 00002650 C* *00002660 C* CALL DOVE2 TO COMPUTE THE DEFLECTION OF THE VERTICAL IN THE PRIME *00002670 C* VERTICAL, AND THE TRANSFORMATION MATRIX REQ'D TO COMPUTE THE *00002680 C* ASSOCIATED VAR-COVAR MATRIX *00002690 C* *00002700 CALL DOVE2(IQ,X,Y,C,IDP,ETA,B,ID,50) 00002710 DO 18 I=1,IQ 00002720 DO 18 J=1,ID 00002730 18 SIG(I,J)=0.D0 00002740 DO 19 I=1,IQ 00002750 DO 19 J=1,IQ 00002760 19 SIGN(I,J)=0.D0 00002770 DO 20 I=1,IQ 00002780 DO 20 J= 1,ID 00002790 DO 20 K=1,ID 00002800 20 SIG(I,J)=B(I,K)*SIGC(K,J)+SIG(I,J) 00002810 DO 21 I=1,IQ 00002820 DO 21 J=1,IQ 00002830 DO 21 K= 1,ID 00002840 21 SIGN(I,J)=SIG(I,K)*B(J,K)+SIGN(I,J) 00002850 DO 76 I=1,IQ 00002860 DO 76 J=1,IQ 00002870 76 SIGN(I,J)=SIGN(I,J)*RHO**2 00002880 DO 78 I=1,IQ 00002890 78 STDE(I)=DSQRT(SIGN(I,I)) 00002900 C* *00002910 C* WRITE THE VAR-COVAR MATRIX OF THE DEFLECTIONS IN THE PRIME *00002920 C* VERTICAL (ETA) *00002930 C* *00002940 WRITE(6,53) 00002950 53 FORMAT('1',10X,'VARIANCE-COVARIANCE MATRIX OF THE DEFLECTIONS IN T00002960 1HE PRIME VERTICAL (ETA)',///) 00002970 CALL MATRIT(SIGN,IQ,IQ,50) 00002980 C* *00002990 C* WRITE THE GEOID HTS AND THE DEFLECTIONS AT THE SELECTED STATIONS *00003000 C* *00003010 WRITE(6,54) 00003020 54 FORMAT('1',//,20X,'GEOID HEIGHTS AND DEFLECTIONS OF THE VERTICAL A00003030 1T SELECTED POSITIONS',////,1X,'STATION',5X,'LATITUDE',8X,'LONGITUD00003040 2E',6X,'GEOID HT.',3X,'STD DEV',5X,'XSI',4X,'STD DEV',5X,'ETA',4X,'00003050 3STD DEV',/,30X,'(+VE E)',9X,'(M)',9X,'(M)',6X,'(SECS)',3X,'(SECS)'00003060 4,4X,'(SECS)',3X,'(SECS)',//) 00003070 DO 23 I=1,IQ 00003080 RX=((X(I)*SF)/R)+RXO 00003090 RY=((Y(I)*SF)/(R*DCOS(RX)))+RYO 00003100 CALL RADARC(RX,LATD,LATM,SLAT) 00003110 CALL RADARC(RY,LOND,LONM,SLON) 00003120 23 WRITE(6,55) IST(I),LATD,LATM,SLAT,LOND,LONM,SLON,GN(I),STDN(I), 00003130 1XSI(I),STDX(I),ETA(I),STDE(I) 00003140 55 FORMAT(2X,I6,I5,I3,F8.4,I6,I3,F8.4,F10.2,F11.2,2(F11.2,F8.2)) 00003150 GO TO 999 00003160 101 WRITE(6,26) 00003170 26 FORMAT(' ',////,5X,'*** YOU ARE TRYING TO READ IN MORE STATIONS TH00003180 1AN ARE ALLOWED - CHECK MAX DIMENSIONS OF X&Y ***') 00003190 999 STOP 00003200 END 00003210 SUBROUTINE ARCRAD (I,M,S,RADS) 00003220 C *$****************************************************************00003230 C * SUBROUTINE 'ARCRAD' CONVERTS SINGLE PRECISION DEGREES,MINUTES *00003240 C * AND SECONDS TO DOUBLE PRECISION RADIANS *00003250 C *$**************************************************************@@00003260 REAL*8 RADS,D,Q,SEC,DABS,RHO/57.29577951308233D0/ 00003270 D=I 00003280 Q=M 00003290 SEC=S 00003300 IF(I)1,2,3 00003310 1 SIGN=-1. 00003320 GO TO 7 00003330 2 IF(M)4,5,3 00003340 4 SIGN=-1. 00003350 GO TO 7 00003360 5 IF(S)6,99,3 00003370 6 SIGN=-1. 00003380 GO TO 7 00003390 3 SIGN=1. 00003400 7 RADS=((DABS(D)+DABS(Q/60.)+DABS(SEC/3600.))/RHO)*SIGN 00003410 GO TO 999 00003420 99 RADS=0.D0 00003430 999 CONTINUE 00003440 RETURN 00003450 END 00003460 SUBROUTINE RADARC (A,I,J,S) 00003470 C ******************************************************************00003480 C * SUBROUTINE 'RADARC' CONVERTS DOUBLE PRECISION RADIANS TO SINGLE*00003490 C * PRECISION DEGREES,MINUTES AND SECONDS *00003500 C ******************************************************************00003510 REAL*8 A,SEC,AD,AJ,RHO/206264.8062470964D0/ 00003520 IF(A)2,1,1 00003530 1 SIGN=1. 00003540 GO TO 3 00003550 2 SIGN=-1. 00003560 A=-A 00003570 3 SEC=A*RHO+0.0005 00003580 I=SEC/3600. 00003590 AD=I 00003600 J=SEC/60.-AD*60. 00003610 AJ=J 00003620 S=SEC-AD*3600.-AJ*60.-0.0005 00003630 IF(I)5,4,5 00003640 4 IF(J)6,7,6 00003650 6 J=J*SIGN 00003660 GO TO 10 00003670 7 S=S*SIGN 00003680 GO TO 10 00003690 5 I=I*SIGN 00003700 10 CONTINUE 00003710 RETURN 00003720 END 00003730 SUBROUTINE GUN(IQ,SF,GO,X,Y,C,IDP,G,B,ID,NR) 00003740 C***********************************************************************00003750 C* SUBROUTINE 'GUN' COMPUTES THE N'S AND FORMS THE TRANSFORMATION *00003760 C* MATRIX REQUIRED TO COMPUTE THE VARIANCE-COVARIANCE MATRIX OF THE N'S 00003770 C***********************************************************************00003780 IMPLICIT REAL*8 (A-H,O-Z) 00003790 REAL*8 X(IQ),Y(IQ),C(ID),B(NR,ID) 00003800 REAL*4 G(IQ) 00003810 CALL ERRSET(208,300,-1,1) 00003820 IN=IDP+1 00003830 DO 1 I=1,IQ 00003840 SUM=0.D0 00003850 M=0 00003860 DO 2 J=1,IN 00003870 DO 3 K=1,IN 00003880 IF(J.EQ.1.AND.K.EQ.1) GO TO 3 00003890 M=M+1 00003900 B(I,M)=X(I)**(J-1)*Y(I)**(K-1) 00003910 SUM=SUM+B(I,M)*C(M) 00003920 3 CONTINUE 00003930 2 CONTINUE 00003940 1 G(I)=GO+SUM*SF 00003950 RETURN 00003960 END 00003970 SUBROUTINE DOVE1(IQ,X,Y,C,IDP,XSI,B,ID,NR) 00003980 C***********************************************************************00003990 C* SUBROUTINE 'DOVE1' COMPUTES THE DEFLECTION OF THE VERTICAL IN THE *00004000 C* PLANE, AND THE TRANSFORMATION MATRIX REQUIRED TO COMPUTE *00004010 C* THE VARIANCE-COVARIANCE MATRIX OF THESE REFLECTIONS 00004020 C***********************************************************************00004030 IMPLICIT REAL*8(A-H,O-Z) 00004040 REAL*8 X(IQ),Y(IQ),C(ID),B(NR,ID),RHO/206264.8062471D0/ 00004050 REAL*4 XSI(IQ) 00004060 CALL ERRSET(208,300,-1,1) 00004070 IN=IDP+1 00004080 DO 1 I=1,IQ 00004090 XSI(I)=0.D0 00004100 M=0 00004110 DO 2 J=1,IN 00004120 DO 3 K=1,IN 00004130 IF(J.EQ.1.AND.K.EQ.1) GO TO 3 00004140 M=M+1 00004150 B(I,M)=(J-1)*X(I)**(J-2)*Y(I)**(K-1) 00004160 XSI(I)=XSI(I)-C(M)*B(I,M) 00004170 3 CONTINUE 00004180 2 CONTINUE 00004190 1 XSI(I)=XSI(I)*RHO 00004200 RETURN 00004210 END 00004220 SUBROUTINE DOVE2(IQ,X,Y,C,IDP,ETA,B,ID,NR) 00004230 C***********************************************************************00004240 C* SUBROUTINE 'DOVE2' COMPUTES THE DEFLECTION OF THE VERICAL IN THE *00004250 C* PRIME VERTICAL AND THE TRANSFORMATION MATRIX REQUIRED TO COMPUTE THE*00004260 C* ASSOCIATED VARIANCE-COVARIANCE MATRIS *00004270 C***********************************************************************00004280 IMPLICIT REAL*8(A-H,O-Z) 00004290 REAL*8 X(IQ),Y(IQ),C(ID),B(NR,ID),RHO/206264.8062471D0/ 00004300 REAL*4 ETA(IQ) 00004310 CALL ERRSET(208,300,-1,1) 00004320 IN=IDP+1 00004330 DO 1 I=1,IQ 00004340 ETA(I)=0.D0 00004350 M=0 00004360 DO 2 J=1,IN 00004370 DO 3 K=1,IN 00004380 IF(J.EQ.1.AND.K.EQ.1) GO TO 3 00004390 M=M+1 00004400 B(I,M)=(K-1)*X(I)**(J-1)*Y(I)**(K-2) 00004410 ETA(I)=ETA(I)-C(M)*B(I,M) 00004420 3 CONTINUE 00004430 2 CONTINUE 00004440 1 ETA(I)=ETA(I)*RHO 00004450 RETURN 00004460 END 00004470 SUBROUTINE MATRIT(A,M,N,NZ) 00004480 C THIS SUBROUTINE PRINTS DOUBLE PRECISION MATRICES 5 COLUMNS AND 40 ROWS00004490 C PAGE 00004500 C A IS THE MATRIX TO BE PRINTED 00004510 C M IS THE ROW DIMENSION OF A 00004520 C N IS THE COLUMN DIMENSION OF A 00004530 C NZ IS THE MAXIMUM ROW DIMENSION OF A IN MAIN 00004540 C IN CASE OF MORE THAN 40 ROWS, THE MATRIX COLUMN CONTINUES IN THE NEXT00004550 C OF THE SAME PAGE 00004560 C A HEADING SHOULD BE PRINTED IN THE MAIN PROGRAM . THE MATRIX WILL STA00004570 C PRINTING ON THE SAME PAGE 00004580 IMPLICIT REAL*8(A-H,O-Z) 00004590 DIMENSION A(NZ,N),IR(2),NC1(5) 00004600 C INITIALIZE 'INT' TO 0 TO INDICATE THAT OUTPUT IS NOT TO START ON A NEW00004610 INT=0 00004620 NC=1 00004630 IF(M.GT.40) GO TO 1 00004640 C PROCEDURE TO WRITE MATRIX IF NO. OF ROWS < OR = 40 00004650 NEW=0 00004660 C PRINT COLUMN NUMBERS 00004670 6 NEW=NEW+5 00004680 IF(N.LT.NEW) NEW=N 00004690 IF (INT.EQ.0) WRITE(6,100)(I,I=NC,NEW) 00004700 IF (INT.NE.0) WRITE(6,200)(I,I=NC,NEW) 00004710 C PRINT M ROWS OF MATRIX 00004720 DO 4 I=1,M 00004730 4 WRITE(6,300) I,(A(I,K),K=NC,NEW) 00004740 IF (N.EQ.NEW) RETURN 00004750 INT=1 00004760 NC=NC+5 00004770 GO TO 6 00004780 C PROCEDURE TO PRINT MATRIX IF NO. OF ROWS IS > 40 00004790 C DETERMINE THE NUMBER OF PAGE COLUMNS EACH MATRIX COLUMN IS TO OCCUPY (00004800 1 NCL=M/40 00004810 IF(MOD(M,40).NE.0) NCL=NCL+1 00004820 MN=1 00004830 NPER=0 00004840 LEFT=0 00004850 25 IF(NCL-5) 83,84,85 00004860 83 NPER=5-NCL 00004870 84 MCL=NCL 00004880 GO TO 86 00004890 85 LEFT=NCL-5 00004900 MCL=5 00004910 86 IR(1)=1 00004920 24 IR(2)=1 00004930 ICL=MCL 00004940 NREM=NPER 00004950 JI=1 00004960 C DETERMINE AND PRINT COLUMN NUMBERS 00004970 12 DO 8 I=JI,ICL 00004980 8 NC1(I)=NC 00004990 21 IF (NREM.EQ.0.OR.NC.EQ.N) GO TO 9 00005000 NC=NC+1 00005010 JI=ICL+1 00005020 IF(NREM-NCL)10,10,11 00005030 10 ICL=NREM+ICL 00005040 LEFT=NCL-NREM 00005050 NREM=0 00005060 GO TO 12 00005070 11 ICL=NCL+ICL 00005080 NREM=NREM-NCL 00005090 GO TO 12 00005100 9 JI=ICL 00005110 IF(INT.EQ.0) WRITE(6,600)(NC1(I),I=1,JI) 00005120 IF(INT.NE.0) WRITE(6,700)(NC1(I),I=1,JI) 00005130 35 ICL=MCL 00005140 22 NREM=NPER 00005150 IF(MN.EQ.0) GO TO 81 00005160 MN=0 00005170 GO TO 82 00005180 81 IF(MOD((IR(1)-1),40).EQ.0.AND.MOD((IR(2)-1),40).EQ.0) GO TO 17 00005190 82 NR=IR(1) 00005200 IR(1)=IR(1)+1 00005210 JI=1 00005220 NM=0 00005230 C PRINT MATRIX VALUES 00005240 16 DO 7 I=JI,ICL 00005250 IF(NR.GT.M) GO TO 18 00005260 GO TO (30,31,32,33,34),I 00005270 30 WRITE(6,400) NR,A(NR,NC1(I)) 00005280 NM=1 00005290 GO TO 7 00005300 31 IF(NM.NE.0) WRITE(6,500) NR,A(NR,NC1(I)) 00005310 IF(NM.EQ.0) WRITE(6,410) NR,A(NR,NC1(I)) 00005320 GO TO 7 00005330 32 WRITE(6,510) NR,A(NR,NC1(I)) 00005340 GO TO 7 00005350 33 WRITE(6,520) NR,A(NR,NC1(I)) 00005360 GO TO 7 00005370 34 WRITE(6,530) NR,A(NR,NC1(I)) 00005380 7 NR=NR+40 00005390 GO TO 28 00005400 18 IF(NC1(ICL).EQ.N.AND.MOD(NR,40).EQ.0) RETURN 00005410 28 IF(NREM.EQ.0.OR.NC1(ICL).EQ.N) GO TO 35 00005420 JI=ICL+1 00005430 IF(NREM-NCL) 14,14,15 00005440 14 ICL=NREM+ICL 00005450 NREM=0 00005460 IF (NCL.GT.2) GO TO 26 00005470 NR=IR(2)-1 00005480 GO TO 16 00005490 26 NR=IR(2) 00005500 IR(2)=IR(2)+1 00005510 GO TO 16 00005520 15 ICL=NCL+ICL 00005530 NREM=NREM-NCL 00005540 NR=IR(2) 00005550 IR(2)=IR(2)+1 00005560 GO TO 16 00005570 17 INT=1 00005580 NPER=0 00005590 IF(LEFT.EQ.0) GO TO 20 00005600 MCL=LEFT 00005610 IR(1)=NR-39 00005620 MN=1 00005630 IF(MCL.LE.5) GO TO 19 00005640 MCL=5 00005650 LEFT=LEFT-5 00005660 GO TO 24 00005670 19 NPER=5-LEFT 00005680 GO TO 24 00005690 20 NC=NC+1 00005700 MN=1 00005710 GO TO 25 00005720 100 FORMAT('0',14X,I3,4(18X,I3)) 00005730 200 FORMAT('1',14X,I3,4(18X,I3)) 00005740 300 FORMAT(' ',I3,5(3X,D19.12)) 00005750 400 FORMAT(' ',I3,1X,D19.12) 00005760 410 FORMAT(' ',23X,I3,1X,D19.12) 00005770 500 FORMAT('+',23X,I3,1X,D19.12) 00005780 510 FORMAT('+',46X,I3,1X,D19.12) 00005790 520 FORMAT('+',69X,I3,1X,D19.12) 00005800 530 FORMAT('+',92X,I3,1X,D19.12) 00005810 600 FORMAT('0',11X,I3,4(20X,I3)) 00005820 700 FORMAT('1',11X,I3,4(20X,I3)) 00005830 END 00005840