IMPLICIT REAL*8(A-H,O-Z) REAL*8 KO,K,L1,L2,MC,MPAZ,MPAZAP,MPDIST,LSK,L3,LSK2,MC2 REAL*8 VCINIT(4,4)/16*0.D0/ DIMENSION DM(2,2),CM(2,2), CM2(2,2), VCDIR(4,4),VCINV(2,2), + DM2(2,2) 30 FORMAT(2F11.2) 31 FORMAT(2(2I4,F8.3)) 32 FORMAT(2F11.2) 35 FORMAT(I5,T40,I2,I3,F8.4,I4,I3,F8.4) 37 FORMAT(I4,I3,F8.4,F10.4,F5.2,F7.4) 40 FORMAT('1',30X, 'DOUBLE STEREOGRAPHIC PROJECTION',////) 41 FORMAT(10X,'ELLIPSOIDAL PARAMETERS',/) 42 FORMAT(10X,'A=',F15.3,5X,'B=',F15.3,///) 43 FORMAT(10X,'GEOGRAPHIC COORDINATES OF THE ORIGIN',/) 44 FORMAT(10X,'PHI=',2I4,F8.3,5X,'LAMDA=',2I4,F8.3,///) 45 FORMAT(10X,'X,Y COORDINATES OF THE ORIGIN',/) 46 FORMAT(10X,'X=',F11.2,5X,'Y=',F11.2,///) 48 FORMAT(10X,'STATION', 6X,'LATITUDE', 9X,'LONGITUDE',10X,'X',12X,'Y 1',//) 49 FORMAT(12X,I5,3X,2I3,F 9.5,3X,2I3,F 9.5,2F13.3,//) 51 FORMAT(10X,'COVARIANCE MATRIX OF MAPPING PLANE COORDINATES IN METR 1ES SQUARED',/) 53 FORMAT(10X,'COVARIANCE MATRIX OF POINTS 1 AND 2, IN METRES SQUARED 1',/) 54 FORMAT(10X,'COVARIANCE MATRIX OF THE AZIMUTH FROM 1 TO 2 IN SECOND 1S SQUARED AND THE DISTANCE IN METRES SQUARED'//) 55 FORMAT(10X,'AZIMUTH FROM 1 TO 2',5X,'AZIMUTH FROM 2 TO 1',5X,'DIST 1ANCE'//) 56 FORMAT(13X,I4,I3,F8.4,10X,I4,I3,F8.4,7X,F8.3,//) 57 FORMAT(10X,'COVARIANCE MATRIX OF PHI 2 AND LAMDA 2, IN SECONDS SQU +ARED',//) 59 FORMAT('0',10X,'T-T CORRECTION',3X,'MERIDIAN CONVERGENCE',3X,'LINE C SCALE') 60 FORMAT('0',10X,I3,I3,F8.4,6X,I3,I3,F8.4,8X,F10.8) 61 FORMAT('1',20X,'DIRECT PROBLEM') 62 FORMAT('1',20X,'INVERSE PROBLEM') 63 FORMAT('1') 141 FORMAT('0',20X,'OBSERVED GEODETIC DATA') 142 FORMAT('0',10X,'AZIMUTH FROM 1 TO 2',5X,'DISTANCE') 143 FORMAT('0',14X,I3,I3,F6.2,8X,F8.3) 144 FORMAT('0',20X,'COMPUTED DATA') 151 FORMAT('1',10X,'COVARIANCE INFORMATION - DIRECT CASE') 152 FORMAT('1',10X,'COVARIANCE INFORMATION - INVERSE CASE') ICODE=1 PI=3.141592653589793D0 C C READ IN DATA - ELLIPSOIDAL PARAMETERS, MAPPING PLANE PARAMETERS, POSITION C DATA AND VARIANCE COVARIANCE INFORMATION OF POINT 1. C READ 30,A,B READ 31,IPND,IPNM,PNS,ILND,ILNM,OLNS READ 32,XO,YO READ 35,NAME,IPD,IPM,PS,LD,LM,OLS READ *,DM(1,1),DM(1,2),DM(2,2) PRINT 40 PRINT 41 PRINT 42,A,B PRINT 43 PRINT 44,IPND,IPNM,PNS,ILND,ILNM,OLNS PRINT 45 PRINT 46,XO,YO C C CONVERT ANGLES TO RADIANS C CALL DMSRAD(IPND,IPNM,PNS,PHIO) CALL DMSRAD(ILND,ILNM,OLNS,ELAMO) KO=0.999912D0 PRINT 61 CALL DMSRAD(IPD,IPM,PS,P1) CALL DMSRAD(LD,LM,OLS,L1) C C TRANSFORM THE GEOGRAPHIC COORDINATES OF POINT 1 TO THE MAPPING PLANE C COORDINATES. C CALL STGINL(PHIO,ELAMO,A,B,R,C1,C2,E,CHIO,SLAMO) CALL ELTSP(P1 , L1 ,E,A,C1,C2,R,CHI,SLAM) CALL SPTPL(CHI,SLAM,XO,YO,KO,CHIO,SLAMO,R,X,Y) X1=X Y1=Y PRINT 48 PRINT 49,NAME,IPD,IPM,PS,LD,LM,OLS,X,Y DSLON=SLAM-SLAMO C C READ IN THE OBSERVATIONAL DATA AND THEIR STANDARD DEVIATIONS. C READ 37,IGD12,IGM12,GS12,GDIST,SDAZ,SDD PRINT 141 PRINT 142 PRINT 143,IGD12,IGM12,GS12,GDIST C C CONVERT THE OBSERVED AZIMUTH TO RADIANS. C CALL DMSRAD(IGD12,IGM12,GS12,GAZ) C C COMPUTE THE MERIDIAN CONVERGENCE ,T-T CORRECTION AND LINE SCALE FACTOR. C CALL PTSKMC(X1,Y1,XO,YO,PHIO,ELAMO,A,B,PSK,MC,KO) MPAZAP=GAZ-MC X2A=X1+DSIN(MPAZAP)*GDIST Y2A=Y1+DCOS(MPAZAP)*GDIST CALL TTLSK(XO,YO,X1,Y1,X2A,Y2A,P1,A,B,KO,TT,LSK) CALL RADMS(TT,ITTD,ITTM,STT) CALL RADMS(MC,IMCD,IMCM,SMC) PRINT 59 PRINT 60,ITTD,ITTM,STT,IMCD,IMCM,SMC,LSK C C SOLVE FOR X2 AND Y2. C CALL MPDIR(X1,Y1,GDIST,GAZ,MC,LSK,TT,MPDIST,MPAZ,X2,Y2) C C PRINT OUT THE COMPUTED DATA. C CALL PLTSP(X2,Y2,XO,YO,KO,R,CHIO,SLAMO,CHI2,SLAM2) CALL SPTEL(CHI2,SLAM2,C1,C2,E,P2,L2) CALL RADMS(P2,IPD2,IPM2,PS2) CALL RADMS(L2,LD2,LM2,SL2) N2=2 PRINT 144 PRINT 48 PRINT 49,N2,IPD2,IPM2,PS2,LD2,LM2,SL2,X2,Y2 C C FORM VARIANCE COVARIANCE MATRIX OF POINT 1 AND THE OBSEVED DATA AND C CONVERT TO THE PROPER UNITS. C DM(2,1)=DM(1,2) DO 14 L=1,2 DO 14 J=1,2 DP2=DM(L,J)*PI**2/180.D0**2 DM(L,J)=DP2/1.296D7 14 CONTINUE CALL STGERP(A,B,P1 ,PHIO,CHI,CHIO,DSLON,CM,DM,ICODE,KO,C1,C2,R) DO 10 I=1,2 DO 10 J=1,2 VCINIT(I,J)=CM(I,J) 10 CONTINUE VCINIT(3,3)=(2.96D-2)**2 VCINIT(4,4)=((5.D0/36.D2)*PI/18.D1)**2 C C COMPUTE THE VARIANCE COVARIANCE MATRIX OF POINTS 1 AND 2. C PRINT 151 CALL ERRMPD(MPAZ,MPDIST,VCINIT,VCDIR) PRINT 51 CALL MOUTD(CM,2,2,2) PRINT 53 CALL MOUTD(VCDIR,4,4,4) DO 11 I=1,2 DO 11 J=1,2 CM2(I,J)=VCDIR(I+2,J+2) 11 CONTINUE ICODE=-1 DSLON=SLAM2-SLAMO CALL STGERP(A,B,P2,PHIO,CHI,CHIO,DSLON,CM2,DM2,ICODE,KO,C1,C2,R) DO 12 I=1,2 DO 12 J=1,2 DM2(I,J)=DM2(I,J)*180.D0**2/PI**2*36.D2**2 12 CONTINUE PRINT 57 CALL MOUTD(DM2,2,2,2) C C ***** INVERSE CASE ***** C PRINT 62 C C SOLVE FOR THE MAPPING PLANE AZIMUTH AND DISTANCE. C CALL TTLSK(XO,YO,X2,Y2,X1,Y1,P2,A,B,KO,TT2,LSK2) CALL PTSKMC(X2,Y2,XO,YO,PHIO,ELAMO,A,B,PSK2,MC2,KO) CALL RADMS(TT2,ITTD,ITTM,STT) CALL RADMS(MC2,IMCD,IMCM,SMC) PRINT 59 PRINT 60,ITTD,ITTM,STT,IMCD,IMCM,SMC,LSK2 CALL MPINV(X1,Y1,X2,Y2,TT,MC,LSK,TT2,MC2,GDIST,GAZ12,GAZ21) CALL RADMS(GAZ12,IGD1,IGM1,SG1) CALL RADMS(GAZ21,IGD2,IGM2,SG2) C C PRINT OUT THE COMPUTED INFORMATION. C PRINT 141 PRINT 48 PRINT 49,NAME,IPD,IPM,PS,LD,LM,OLS,X,Y PRINT 49,N2,IPD2,IPM2,PS2,LD2,LM2,SL2,X2,Y2 PRINT 144 PRINT 55 PRINT 56,IGD1,IGM1,SG1,IGD2,IGM2,SG2,GDIST C C COMPUTE THE VARIANCE COVARIANCE MATRIX OF THE COMPUTED MAPPING PLANE AZIMUTH C AND DISTANCE. C PRINT152 CALL ERRMPI(MPDIST,X1,Y1,X2,Y2,VCDIR,VCINV) PRINT 54 C C CONVERT TO THE PROPER UNITS. C VCINV(1,2)=VCINV(1,2)*180.D0/PI*3600.D0 VCINV(2,1)=VCINV(1,2) VCINV(1,1)=VCINV(1,1)*180.D0**2/PI**2*3600.0D0**2 CALL MOUTD(VCINV,2,2,2) PRINT 63 STOP END SUBROUTINE PTSKMC(X,Y,XO,YO,PHIO,ELAMO,A,B,PSK,MC,KO) C C THIS SUBROUTINE WILL COMPUTE THE POINT SCALE FACTOR AND THE MERIDIAN C CONVERGENCE FOR THE DOUBLE STEREOGRAPHIC PROJECTION. C C INPUT: C X,Y - COORDINATES OF THE POINT (IN METRES). C XO,YO - COORDINATES OF THE ORIGIN (IN METRES). C PHIO,ELAMO - PHI AND LAMDA OF THE ORIGIN (IN RADIANS). C A,B - SEMI-MAJOR AND SEMI-MINOR AXES OF THE ELLIPSOID (IN METRES) C C OUTPUT: C PSK - POINT SCALE FACTOR. C MC - MERIDIAN CONVERGENCE (IN RADIANS). C WRITTEN BY G.BOWIE, JAN. 1978. C IMPLICIT REAL*8(A-Z) CALL STGINL(PHIO,ELAMO,A,B,R,C1,C2,E,CHIO,SLAMO) CALL PLTSP(X,Y,XO,YO,KO,R,CHIO,SLAMO,CHI,SLAM) CALL SPTEL(CHI,SLAM,C1,C2,E,PHI,ELAM) CC=DCOS(CHI) CP=DCOS(PHI) CCO=DCOS(CHIO) CDL=DCOS(SLAM-SLAMO) SC=DSIN(CHI) SP=DSIN(PHI) SDL=DSIN(SLAM-SLAMO) SCO=DSIN(CHIO) N=A/DSQRT(1.D0-E**2*SP**2) PSK=2.D0*C1*R*CC/(N*CP*(1.D0+SCO*SC+CCO*CC*CDL)) NUM=SDL*(SC+SCO) DEN=CC*CCO+(1.D0+SC*SCO)*CDL MC=DATAN2(NUM,DEN) RETURN END SUBROUTINE TTLSK(XO,YO,X1,Y1,X2,Y2,PHIO,A,B,KO,TT,LSK) C C THIS SUBROUTINE WILL COMPUTE THE T-T CORRECTION AND THE C LINE SCALE FACTOR FOR THE DOUBLE STEREOGRAPHIC PROJECTION. C C INPUT: C XO,YO - COORDINATES OF THE ORIGIN (IN METRES). C X1,Y1 - COORDINATES OFTHE KNOWN POINT (IN METRES). C X2,Y2 - APPROXIMATE COORDINATES OF THE UNKNOWN POINT (IN METRES). C PHIO - PHI OF THE ORIGIN (IN RADIANS). C A,B - SEMI-MAJOR AND SEMI-MINOR AXES OF THE ELLIPSOID (IN METRES) C KO - SCALE FACTOR AT THE ORIGIN. C C OUTPUT: C TT - T-T CORRECTION (IN RADIANS). C LSK - LINE SCALE FACTOR. C C WRITTEN BY G. BOWIE, JAN. 1977. C IMPLICIT REAL*8(A-Z) PI=3.141592653589793D0 X1=X1-XO Y1=Y1-YO X2=X2-XO Y2=Y2-YO E2=(A*A-B*B)/(A*A) R=A*DSQRT(1.D0-E2)/(1.D0-E2*DSIN(PHIO)**2) NUM=X2*Y1-X1*Y2 DEN=X1*X2+Y1*Y2+(KO*R*2.D0)**2 TT=DATAN2(NUM,DEN) TT=DABS(TT) AZ01=DATAN2(X1,Y1) AZ10=DATAN2(-X1,-Y1) AZ12=DATAN2(X2-X1,Y2-Y1) IF(AZ01.LT.0.0) AZ01=AZ01+PI*2.D0 IF(AZ10.LT.0.0) AZ10=AZ10+PI*2.D0 IF(AZ12.LT.0.0) AZ12=AZ12+PI*2.D0 IF(AZ10.LT.AZ01) GO TO 3 IF(AZ12.GT.AZ01) GO TO 1 GO TO 2 1 IF(AZ12.LT.AZ10) GO TO 5 2 TT=-TT GO TO 5 3 IF(AZ12.GT.AZ10) GO TO 4 GO TO 5 4 IF(AZ12.GT.AZ01) GO TO 5 GO TO 2 5 CONTINUE TT=-TT Y1=Y1+YO X1=X1+XO Y2=Y2+YO X2=X2+XO DEN2=4.D0*KO*(R**2) K1=KO+((X1-XO)**2+(Y1-YO)**2)/DEN2 K2=KO+((X2-XO)**2+(Y2-YO)**2)/DEN2 XM=(X1+X2)/2.D0 YM=(Y1+Y2)/2.D0 KM=KO+((XM-XO)**2+(YM-YO)**2)/DEN2 LSK=1.D0/((1.D0/6.D0)*(1.D0/K1+4.D0/KM+1.D0/K2)) RETURN END SUBROUTINE STGINL (PHIO,ELAMO,A,B,R,C1,C2,E,CHIO,SLAMO) C C THIS ROUTINE COMPUTES THE INITIAL VALUES TO BE USED IN C THE STEREOGRAPHIC DOUBLE PROJECTION SUBROUTINES. C C INPUT: C C PHIO - ELLIPSOID LATITUDE OF THE ORIGIN OF THE PROJECTION C IN RADIANS. C ELAMO - ELLIPSOIDAL LNGITUDE (POSITIVE EAST OF GREENWICH) C OF THE ORIGIN OF THE PROJECTION IN RADIANS. C A,B - SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE C ELLIPSOID, IN METRES. C C OUTPUT: C C R - RADIUS OF THE CONFORMAL SPHERE, IN METRES. C C1 - CONSTANT USED IN THE TRANSFORMATIONS BETWEEN THE C ELLIPSOID AND THE CONFORMAL SPHERE. C C2 - CONSTANT FOR THE SAME USE AS C1. C E - FIRST ECCENTRICITY OF THE ELLIPSOID. C CHIO - SPHERICAL LATITUDE OF THE ORIGIN OF THE PROJECTION, C IN RADIANS. C SLAMO - SPERICAL LONGITUDE OF THE ORIGIN OF THE PROJECTION C IN RADIANS. C C WRITTEN BY R.R.STEEVES,JULY,1977. C IMPLICIT REAL*8(A-H,O-Z) E2=(A*A-B*B)/(A*A) E=DSQRT(E2) SP=DSIN(PHIO) R=A*DSQRT(1.D0-E2)/(1.D0-E2*SP**2) C1=DSQRT(1.D0+E2 /(1.D0-E2)*DCOS(PHIO)**4) CHIO=DARSIN(SP/C1) SLAMO=C1*ELAMO PI=3.141592653589793D0 C2=DTAN(PI/4.D0+CHIO/2.D0)/(DTAN(PI/4.D0+PHIO/2.D0)*((1.D0-E*SP)/ 1 (1.D0+E*SP))**(E/2.D0))**C1 RETURN END SUBROUTINE ELTSP(PHI,ELAM,E,A,C1,C2,R,CHI,SLAM) C THIS ROUTINE TRANSFORMS ELLIPSOIDAL COORDINATES PHI,ELAM TO C SPHERICAL (CONFORMAL SPHERE) COORDINATES CHI,SLAM AND COMPUTES C THE CORRESPONDING POINT SCALE FACTOR ESK (ELLIPSOID TO SPHERE). C THE POINT SCALE FACTOR AT THE ORIGIN OF THIS CONFORMAL PROJECTION C IS UNITY. C C INPUT: C C PHI - ELLIPSOIDAL LATITUDE OF THE POINT, IN RADIANS. C ELAM - ELLIPSOIDAL LONGITUDE OF THE POINT, IN RADIANS. C (POSITIVE EAST OF GREENWICH). C E - FIRST ECCENTRICITY OF THE ELLIPSOID (COMPUTED IN C SUBROUTINE STGINL). C A - SEMI-MAJOR AXES OF THE REFERENCE ELLIPSOID. C C1 - CONSTANT COMPUTED IN STGINL. C C2 - CONSTANT COMPUTED IN STGINL. C R - RADIUS OF THE CONFORMAL SPHERE (COMPUTED IN STGINL). C C OUTPUT: C C CHI - SPHERICAL LATITUDE OF THE POINT, IN RADIANS. C SLAM - SPHERICAL LONGITUDE OF THE POINT, IN RADIANS. C C WRITTEN BY R.R.STEEVES,JULY,1977. C IMPLICIT REAL*8(A-H,O-Z) SP=DSIN(PHI) PI4=3.141592653589793D0/4.D0 CHI=DATAN(C2*(DTAN(PI4+PHI/2.D0)*((1.D0-E*SP)/(1.D0+E*SP))**(E/2.D 1 0))**C1) CHI=2.D0*(CHI-PI4) SLAM=C1*ELAM RN=A/DSQRT(1.D0-E**2*SP**2) ESK=C1*R*DCOS(CHI)/RN/DCOS(PHI) RETURN END SUBROUTINE SPTEL(CHI,SLAM,C1,C2,E,PHI,ELAM) C C THIS ROUTINE TRANSFORMS SPHERICAL (CONFORMAL SPHERE) COORDINATES C CHI,SLAM TO ELLIPSOIDAL COORDINATES PHI,ELAM USING A NEWTON- C RAPHSON ITERATION. C C INPUT: C C CHI - SPHERICAL LATITUDE OF THE POINT, IN RADIANS. C SLAM - SPHERICAL LONGITUDE OF THE POINT, IN RADIANS. C E - FIRST ECCENTRICITY OF THE ELLIPSOID (COMPUTED IN C SUBROUTINE STGINL). C C1 - CONSTANT COMPUTED IN STGINL. C C2 - CONSTANT COMPUTED IN STGINL. C C OUTPUT: C C PHI - ELLIPSOIDAL LATITUDE OF THE POINT, IN RADIANS. C ELAM - ELLIPSOIDAL LONGITUDE OF THE POINT, IN RADIANS. C C WRITTEN BY R.R.STEEVES,JULY,1977. C IMPLICIT REAL*8(A-H,O-Z) PI4=3.141592653589793D0/4.D0 PHI=CHI 1 ESP=E*DSIN(PHI) P2=((1.D0-ESP)/(1.D0+ESP))**(E/2.D0) P1=DTAN(PI4+PHI/2.D0) F=C2*(P1*P2)**C1-DTAN(PI4+CHI/2.D0) FP=C1*C2*(P1*P2)**(C1-1.D0)*P2*(1.D0/2.D0/DCOS(PI4+PHI/2.D0)**2-E* 1 *2*DCOS (PHI)/(1.D0-ESP**2)*DTAN(PI4+PHI/2.D0)) DPHI=F/FP PHI=PHI-DPHI IF(DABS(DPHI).GT.1.D-11) GO TO 1 ELAM=SLAM/C1 RETURN END SUBROUTINE SPTPL(CHI,SLAM,XO,YO,KO,CHIO,SLAMO,R,X,Y) C C THIS ROUTINE TRANSFORMS SPHERICAL COORDINATES CHI,SLAM TO C STEREOGRAPHIC GRID COORDINATES X,Y. C C INPUT: C C CHI - SPHERICAL LATITUDE OF THE POINT, IN RADIANS. C SLAM - SPHERICAL LONGITUDE OF THE POINT, IN RADIANS. C (POSITIVE EAST OF GREENWICH) C XO - FALSE EASTING OF THE ORIGIN OF THE PROJECTION C YO - FALSE NORTHING OF THE ORIGIN OF THE PROJECTION. C KO - POINT SCALE FACTOR AT THE ORIGIN OF THE PROJECTION. C (FROM SPHERE TO PLANE) C CHIO - SPHERICAL LATITUDE OF THE ORIGIN, IN RADIANS. C SLAMO - SPHERICAL LONGITUDE OF THE ORIGIN, IN RADIANS. C R - RADIUS OF THE SPHERE. C C OUTPUT: C C X - STEREOGRAPHIC GRID EASTING. C Y - STEREOGRAPHIC GRID NORTHING. C K - POINT SCALE FACTOR AT THE POINT, GOING FROM THE SPHERE C TO THE PLANE. C C - MERIDIAN CONVERGENCE AT THE POINT, IN RADIANS. C C WRITTEN BY R.R.STEEVES,JULY,1977. C IMPLICIT REAL*8(A-H,O-Z) REAL*8 KO,K RO=2.D0*KO*R DLAM=SLAM-SLAMO CC=DCOS(CHI) SC=DSIN(CHI) CCO=DCOS(CHIO) SCO=DSIN(CHIO) SDL=DSIN(DLAM) CDL=DCOS(DLAM) DEN=1.D0+SCO*SC+CCO*CC*CDL X=XO+RO*CC*SDL/DEN Y=YO+RO*(SC*CCO-CC*SCO*CDL)/DEN K=2.D0*KO/DEN C=DATAN((SDL*(SC+SCO))/(CC*CCO+(1.D0+SC*SCO)*CDL)) RETURN END SUBROUTINE PLTSP (X,Y,XO,YO,KO,R,CHIO,SLAMO,CHI,SLAM) C C THIS ROUTINE TRANSFORMS STEREOGRAPHIC GRID COORDINATES X,Y TO C SPHERICAL COORDINATES CHI,SLAM. C C INPUT: C C X - STEREOGRAPHIC GRID EASTING. C Y - STEREOGRAPHIC GRID NORTHING. C XO - FALSE EASTING. C YO - FALSE NORTHING. C KO - POINT SCALE FACTOR AT THE ORIGIN (FROM SPHERE TO C PLANE). C R - RADIUS OF THE SPHERE. C CHIO - SPHERICAL LATITUDE OF THE ORIGIN, IN RADIANS. C SLAMO - SPHERICAL LONGITUDE OF THE ORIGIN, IN RADIANS. C (POSITIVE EAST OF GREENWICH.) C C OUTPUT: C C CHI - SPHERICAL LATITUDE OF THE POINT, IN RADIANS. C SLAM - SPHERICAL LONGITUDE OF THE POINT, IN RADIANS. C C WRITTEN BY R.R.STEEVES,JULY,1977. C IMPLICIT REAL*8(A-H,O-Z) REAL*8 KO,K YY=(Y-YO)/KO XX=(X-XO)/KO S=DSQRT(XX**2+YY**2) DEL=2.D0*DATAN(S/2.D0/R) CB=1.D0 IF(S.GT.1.D-50)CB=XX/S SB=0.D0 IF(S.GT.1.D-50)SB=YY/S CD=DCOS(DEL) SD=DSIN(DEL) CHI=DARSIN(DSIN(CHIO)*CD+DCOS(CHIO)*SD*SB) SLAM=SLAMO+DARSIN(CB*SD/DCOS(CHI)) RETURN END SUBROUTINE STGERP(A,B,PHI,PHIO,CHI,CHIO,DSLON,CM,DM,ICODE,KO,C1,C2 1 ,R) C C C THIS ROUTINE COMPUTES THE TWO BY TWO COVARIANCE MATRIX OF THE C STEREOGRAPHIC (DOUBLE PROJECTION) PLANE COORDINATES GIVEN THE C CORRESPONDING TWO BY TWO COVARIANCE MATRIX OF THE ELLIPSOIDAL C COORDINATES OF THE POINT; AND CONVERSELY. C C INPUT: C C A,B - SEMI-MAJOR AND SEMI-MINOR AXES OF THE REFERENCE C ELLIPSOID, IN METRES. C PHI - ELLIPSOIDAL LATITUDE OF THE POINT, IN RADIANS. C PHIO - ELLIPSOIDAL LATITUDE OF THE ORIGIN OF THE C PROJECTION, IN RADIANS. C CHI - SPHERICAL LATITUDE OF THE POINT, IN RADIANS. C CHIO - SPHERICAL LATITUDE OF THE ORIGIN OF THE PROJECTION, C IN RADIANS. C DSLON - SPHERICAL LONGITUDE OF THE POINT MINUS THE C SPHERICAL LONGITUDE OF THE ORIGIN OF THE C PROJECTION, IN RADIANS. (LONGITUDE IS MEASURED C POSITIVE EAST). C CM - COVARIANCE MATRIX (2X2) OF THE PLANE STEREOGRAPHIC C GRID COORDINATES, IN METRES SQUARED. (SEE ICODE C BELOW). C DM - COVARIANCE MATRIX (2X2) OF THE ELLIPSOIDAL C COORDINATES, IN RADIANSSQUARED. (SEE ICODE). C ICODE - MAY BE EITHER 1 OR -1. C - WHEN ICODE IS 1 CM WILL BE COMPUTED GIVEN DM. C (DM IS INPUT AND CM IS OUTPUT). C - WHEN ICODE IS -1 DM WILL BE COMPUTED GIVEN CM. C KO - POINT SCALE FACTOR AT THE ORIGIN. C C1 - A CONSTANT COMPUTED IN STGINL. C C2 - A CONSTANT COMPUTED IN STGINL. C R - RADIUS OF THE CONFORMAL SPHERE,(COMPUTED IN STGINL). C C OUTPUT: C C CM OR DM DEPENDING ON THE VALUE OF ICODE; (1 OR -1 ) C C THIS ROUTINE UTILIZES RIGOROUS EQUATIONS FOR THE ERROR C PROPAGATION; APPROXIMATIONS MAY BE MORE EFFICIENT AND YEILD C SUFFICIENT ACCURACY IN THE RESULTS. C C WRITTEN BY R.R.STEEVES,JULY,1977. C IMPLICIT REAL* 8 (A-H,J-Z) DIMENSION CM(2,2),DM(2,2) E2=(A*A-B*B)/(A*A) E=DSQRT(E2) ESP=E*DSIN(PHI) RR=3.141592653589793D0/4.D0+PHI/2.D0 Q=((1.D0-ESP)/(1.D0+ESP))**(E/2.D0) R2=2.D0*KO*R SC=DSIN(CHI) SCO=DSIN(CHIO) CC=DCOS(CHI) CCO=DCOS(CHIO) CDS=DCOS(DSLON) M=R2*DSIN (DSLON)*(SC+SCO) N=(1.D0+SC*SCO+CC*CCO*CDS)**2 P=R2*(CDS*(1.D0+SC*SCO)+CC*CCO) TR=DTAN(RR) U=C1*C2*(Q*TR)**(C1-1.D0)*Q*(1.D0/2.D0/DCOS(RR)**2-E2*DCOS(PHI)*TR 1/(1.D0-E2*DSIN(PHI)**2))/(1.D0+C2**2*(Q*TR)**(2.D0*C1)) U=U*2.D0 CNC=C1/N*CC J11=-M/N*U J12=CNC*P J21=P/N*U J22=CNC*M IF(ICODE.LT.0)GOTO1 CM(1,1)=J11**2*DM(1,1)+2.D0*J11*J12*DM(1,2)+J12**2*DM(2,2) CM(1,2)=J11*J21*DM(1,1)+J12*J21*DM(1,2)+J11*J22*DM(1,2)+J12*J22*DM 1 (2,2) CM(2,1)=CM(1,2) CM(2,2)=J21**2*DM(1,1)+2.D0*J21*J22*DM(1,2)+J22**2*DM(2,2) GO TO 10 1 CONTINUE DET=J11*J22-J21*J12 J1=J11 J11=J22/DET J12=-J12/DET J21=-J21/DET J22=J1/DET DM(1,1)=J11**2*CM(1,1)+2.D0*J11*J12*CM(1,2)+J12**2*CM(2,2) DM(1,2)=J11*J21*CM(1,1)+J12*J21*CM(1,2)+J11*J22*CM(1,2)+J12*J22*CM 1 (2,2) DM(2,1)=DM(1,2) DM(2,2)=J21**2*CM(1,1)+2.D0*J21*J22*CM(1,2)+J22**2*CM(2,2) 10 RETURN END