C***********************************************************************00000040 C* *00000050 C* *00000060 C* ************************* *00000070 C* * * *00000080 C* * C O M P A R E * *00000090 C* * * *00000100 C* ************************* *00000110 C* *00000120 C* *00000130 C* PROGRAM 'COMPARE' COMPUTES THE SCALE CHANGE AND ROTATION BETWEEN *00000140 C* TWO SETS OF ADJUSTED GEODETIC COORDINATES OF THE SAME NETWORK. THE *00000150 C* SCALE CHANGE AND ROTATION ARE COMPUTED WITH RESPECT TO THE NETWORK *00000160 C* FIXED POINT (SEE BELOW). THE PROGRAM COMPUTES FOR AN *00000170 C* ENTIRE NETWORK (UP TO 1000 POINTS) AND ANY NUMBER OF SUBNETWORKS *00000180 C* CONTAINED IN THE NETWORK (UP TO 500 POINTS EACH). THE REQUIRED *00000190 C* INPUT IS AS FOLLOWS: *00000200 C* *00000210 C* CARD 1: FORMAT (20A4); JOB DESCRIPTION *00000220 C* *00000230 C* CARD 2: FORMAT(I5,I3,2I2); (1) NUMBER OF STATIONS IN NETWORK *00000240 C* (2) NUMBER OF SUBNETWORKS *00000250 C* (3) UNIT NUMBER FOR 'OLD' COORDINATES *00000260 C* (4) UNIT NUMBER FOR 'NEW' COORDINATES *00000270 C* *00000280 C* -STATION NUMBERS AND ADJUSTED COORDINATES CAN BE READ FROM *00000290 C* AN AUXILLARY DEVICE (DISK OR TAPE) OR FROM CARDS, FOR CARDS *00000300 C* UNIT NUMBER MUST BE 5, FOR AUXILLARY DEVICE UNIT NUMBER IS *00000310 C* ANY VALID FORTRAN UNIT NUMBER WITH THE APPROPRIATE 'DD' *00000320 C* STATEMENT. *00000330 C* *00000340 C* CARD 3: FORMAT(I9,1X,12A4); STATION NUMBER AND DESCRIPTION OF *00000350 C* FIXED POINT *00000360 C* *00000370 C* CARD 4: FORMAT(I9,1X,2I3,F7.0,I4,I3,F7.0); STATION NUMBER, *00000380 C* LATITUDE,LONGITUDE OF NETWORK POINT FOR 'OLD' *00000390 C* COORDINATES *00000400 C* *00000410 C* CARD 5: SAME AS CARD 4 BUT FOR 'NEW' COORDINATES *00000420 C* *00000430 C* NOTE:-CARDS 4 & 5 ARE REQUIRED ONLY IF UNIT NUMBERS ON CARD *00000440 C* NUMBER 2 ARE 5 *00000450 C* -IF CARDS 4 & 5 ARE REQUIRED THEY ARE REPEATED ONCE FOR *00000460 C* EACH NETWORK POINT *00000470 C* -IF INPUT IS NOT FROM CARDS (CARD TYPES 4 & 5) THE FORMAT *00000480 C* TO USE IN CREATING THE DISC (TAPE) DATA SET IS *00000490 C* 'WRITE(IUNIT)STATION NO.,LATITUDE,LONGITUDE' WHERE *00000500 C* STATION NO. IS AN INTEGER,LATITUDE AND LONGITUDE ARE REAL*8*00000510 C* VARIABLES WITH THE COORDINATES IN RADIANS (CREADED BY *00000520 C* PROGRAM 'SOLVE') *00000530 C* *00000540 C* CARD 6: FORMAT(I9); STATION NUMBER OF POINTS IN SUBNETWORKS *00000550 C* *00000560 C* NOTE:-CARD 6 IS REQUIRED ONCE FOR EACH SUBNETWORK POINT *00000570 C* -EACH SUBNETWORK IS SEPERATED BY A BLANK CARD *00000580 C* -IF THERE ARE NO SUBNETWORKS CARD 6 IS NOT REQUIRED *00000590 C* *00000600 C* GENERAL NOTE: PROGRAM COMPUTES 'NEW - OLD' WHICH ARE DISTINGUISHED *00000610 C* ONLY BY THE ORDER IN WHICH THE COORDINATES ARE READ IN. *00000620 C* *00000630 C* *00000640 C* C. CHAMBERLAIN SEPT/1975 *00000650 C* *00000660 C***********************************************************************00000670 IMPLICIT REAL*8(A-H,O-Z) 00000680 REAL*4 S1,S2,S3,S4,ABS 00000690 DIMENSION DJOB(20),DESC(12),ISTNO(1000),ISTNN(1000),PHIO(1000), 00000700 1PHIN(1000),XLONO(1000),XLONN(1000),DPHIR(1000),DPHIM(1000), 00000710 2DXLONR(1000),DXLONM(1000),DSCALE(1000),DAZ(1000),VEC(1000), 00000720 3ISTNS(500),IPOINT(500) 00000730 C 00000740 C CLARKE 1866 ELLIPSOID 00000750 C 00000760 A=6378206.4D0 00000770 B=6356583.8D0 00000780 ESQ=(A*A-B*B)/A/A 00000790 C 00000800 C COUNTER FOR OUTPUT CONTROL 00000810 C 00000820 N=1 00000830 C 00000840 C READ AND PRINT CARDS 1,2,3 00000850 C 00000860 READ 1000,DJOB 00000870 PRINT 1010,DJOB 00000880 READ 1020,NSTN,NSUBN,IUO,IUN 00000890 READ 1030,IFIX,DESC 00000900 PRINT 1040,DESC,IFIX 00000910 C 00000920 C READ IN STATION NUMBERS AND COORDINATES 00000930 C 00000940 DO 50 I=1,NSTN 00000950 IF(MOD(I,26).EQ.1)PRINT 1050 00000960 IF(IUO.EQ.5)READ 1060,ISTNO(I),K,L,S1,KK,LL,S2 00000970 IF(IUO.NE.5)READ(IUO)ISTNO(I),PHIO(I),XLONO(I) 00000980 IF(IUO.EQ.5)GO TO 10 00000990 CALL RADARC(PHIO(I),K,L,S1) 00001000 CALL RADARC(XLONO(I),KK,LL,S2) 00001010 GO TO 20 00001020 10 CALL ARCRAD(K,L,S1,PHIO(I)) 00001030 CALL ARCRAD(KK,LL,S2,XLONO(I)) 00001040 20 IF(IUN.EQ.5)READ 1060,ISTNN(I),J,M,S3,JJ,MM,S4 00001050 IF(IUN.NE.5)READ(IUN)ISTNN(I),PHIN(I),XLONN(I) 00001060 IF(IUN.EQ.5)GO TO 30 00001070 CALL RADARC(PHIN(I),J,M,S3) 00001080 CALL RADARC(XLONN(I),JJ,MM,S4) 00001090 GO TO 40 00001100 30 CALL ARCRAD(J,M,S3,PHIN(I)) 00001110 CALL ARCRAD(JJ,MM,S4,XLONN(I)) 00001120 40 IF(ISTNO(I).NE.ISTNN(I))GO TO 999 00001130 50 PRINT 1070,ISTNO(I),K,L,S1,KK,LL,S2,J,M,S3,JJ,MM,S4 00001140 C 00001150 C FIND SEQUENCE NUMBER AND COORDINATES OF FIXED POINT -- ERROR 00001160 C EXIT IF NOT FOUND 00001170 C 00001180 DO 55 I=1,NSTN 00001190 IF(ISTNO(I).EQ.IFIX)GO TO 56 00001200 55 CONTINUE 00001210 GO TO 996 00001220 56 IFIXS=I 00001230 PHI=PHIO(IFIXS) 00001240 XLON=XLONO(IFIXS) 00001250 C 00001260 C DO REQUIRED COMPUTATIONS 00001270 C 00001280 SUM1=0.0D0 00001290 SUM2=0.0D0 00001300 SUM3=0.0D0 00001310 SUM4=0.0D0 00001320 DO 60 I=1,NSTN 00001330 IF(I.EQ.IFIXS)GO TO 60 00001340 C 00001350 C COORDINATE DIFFERENCES IN RADIANS (NEW - OLD) 00001360 C 00001370 DPHIR(I)=PHIN(I)-PHIO(I) 00001380 DXLONR(I)=XLONN(I)-XLONO(I) 00001390 C 00001400 C COORDINATE DIFFERENCES IN METRES AND LENGTH OF 'MISCLOSURE' VECTOR 00001410 C 00001420 DEN=DSQRT(1.0D0-ESQ*DSIN(PHIO(I))**2) 00001430 RN=A/DEN 00001440 RM=A*(1.0D0-ESQ)/DEN**3 00001450 DPHIM(I)=DPHIR(I)*RM 00001460 DXLONM(I)=DXLONR(I)*RN*DCOS(PHIO(I)) 00001470 VEC(I)=DSQRT(DPHIM(I)**2+DXLONM(I)**2) 00001480 C 00001490 C COMPUTE SCALE AND AZIMUTH CHANGE BY COMPUTING AZIMUTH & DISTANCE 00001500 C BETWEEN THE NETWORK POINTS AND THE 'FIXED POINT' 00001510 C 00001520 CALL VININV(A,B,PHI,XLON,PHIN(I),XLONN(I),AZN,AZ,DISN) 00001530 CALL VININV(A,B,PHI,XLON,PHIO(I),XLONO(I),AZO,AZ,DISO) 00001540 DSCALE(I)=(DISN-DISO)/DISN*1.0D6 00001550 DAZ(I)=AZN-AZO 00001560 SUM1=SUM1+DSCALE(I) 00001570 SUM2=SUM2+DSCALE(I)**2 00001580 SUM3=SUM3+DAZ(I) 00001590 SUM4=SUM4+DAZ(I)**2 00001600 60 CONTINUE 00001610 C 00001620 C PRINT OUT RESULTS--FLAG SCALE CHANGES > 10 PPM AND ROTATIONS > 1 SEC 00001630 C 00001640 DO 70 I=1,NSTN 00001650 IF(MOD(I,26).EQ.1)PRINT 1080 00001660 IF(I.NE.IFIXS)GO TO 65 00001670 PRINT 1085,IFIX 00001680 GO TO 70 00001690 65 CALL RADARC(DAZ(I),J,K,S1) 00001700 CALL RADARC(DXLONR(I),J,K,S3) 00001710 CALL RADARC(DPHIR(I),J,K,S2) 00001720 PRINT 1090,ISTNN(I),S2,S3,DPHIM(I),DXLONM(I),VEC(I),DSCALE(I00001730 1 ),S1 00001740 IF(DABS(DSCALE(I)).GT.10.0D0.OR.ABS(S1).GT.1.0)PRINT 1100 00001750 70 CONTINUE 00001760 C 00001770 C COMPUTE AND PRINT MEAN SCALE AND ROTATION 00001780 C 00001790 DSCM=SUM1/DFLOAT(NSTN-1) 00001800 DAZM=SUM3/DFLOAT(NSTN-1) 00001810 STDS=DSQRT((DFLOAT(NSTN)*SUM2-SUM1**2)/DFLOAT(NSTN)/DFLOAT(NSTN-00001820 1 1)) 00001830 STDA=DSQRT((DFLOAT(NSTN)*SUM4-SUM3**2)/DFLOAT(NSTN)/DFLOAT(NSTN-00001840 1 1)) 00001850 CALL RADARC(DAZM,I,J,S1) 00001860 CALL RADARC(STDA,I,J,S2) 00001870 PRINT 1110,DSCM,STDS,S1,S2 00001880 C 00001890 C COMPUTE ANY SUBNETWORKS 00001900 C 00001910 90 IF(NSUBN.EQ.0)GO TO 998 00001920 NSUBN=NSUBN-1 00001930 C 00001940 C READ SUBNETWORK STATION NUMBERS UNTIL BLANK CARD ENCOUNTERED 00001950 C 00001960 NSTNSN=1 00001970 100 READ 1120,ISTNS(NSTNSN) 00001980 IF(ISTNS(NSTNSN).EQ.0)GO TO 110 00001990 NSTNSN=NSTNSN+1 00002000 GO TO 100 00002010 110 NSTNSN=NSTNSN-1 00002020 PRINT 1130,N 00002030 N=N+1 00002040 C 00002050 C IDENTIFY SUBNETWORK POINTS--ERROR EXIT IF THEY ARE NOT IN NETWORK 00002060 C 00002070 DO 140 I=1,NSTNSN 00002080 DO 120 J=1,NSTN 00002090 IF(ISTNO(J).EQ.ISTNS(I))GO TO 130 00002100 120 CONTINUE 00002110 GO TO 997 00002120 130 IPOINT(I)=J 00002130 140 CONTINUE 00002140 C 00002150 C CHECK IF FIXED POINT IN SUBNETWORK 00002160 C 00002170 IFLAG=0 00002180 IFIXSS=0 00002190 DO 145 I=1,NSTNSN 00002200 IF(ISTNS(I).EQ.IFIX)GO TO 146 00002210 145 CONTINUE 00002220 GO TO 147 00002230 146 IFIXSS=I 00002240 IFLAG=1 00002250 C 00002260 C PRINT OUT COORDINATES OF SUBNETWORK 00002270 C 00002280 147 DO 150 I=1,NSTNSN 00002290 IF(MOD(I,26).EQ.1)PRINT 1050 00002300 NN=IPOINT(I) 00002310 CALL RADARC(PHIO(NN),K,L,S1) 00002320 CALL RADARC(XLONO(NN),KK,LL,S2) 00002330 CALL RADARC(PHIN(NN),J,M,S3) 00002340 CALL RADARC(XLONN(NN),JJ,MM,S4) 00002350 150 PRINT 1070,ISTNS(I),K,L,S1,KK,LL,S2,J,M,S3,JJ,MM,S4 00002360 C 00002370 C PRINT OUT SCALE AND ROTATION CHANGES FOR SUBNETWORK (FLAG LARGE VALUES00002380 C 00002390 SUM1=0.0D0 00002400 SUM2=0.0D0 00002410 SUM3=0.0D0 00002420 SUM4=0.0D0 00002430 DO 160 I=1,NSTNSN 00002440 IF(MOD(I,26).EQ.1)PRINT 1080 00002450 IF(I.NE.IFIXSS)GO TO 155 00002460 PRINT 1085,IFIX 00002470 GO TO 160 00002480 155 NN=IPOINT(I) 00002490 CALL RADARC(DAZ(NN),J,K,S1) 00002500 CALL RADARC(DPHIR(NN),J,K,S2) 00002510 CALL RADARC(DXLONR(NN),J,K,S3) 00002520 PRINT 1090,ISTNS(I),S2,S3,DPHIM(NN),DXLONM(NN),VEC(NN), 00002530 1 DSCALE(NN),S1 00002540 IF(DABS(DSCALE(NN)).GT.10.0D0.OR.ABS(S1).GT.1.)PRINT110000002550 SUM1=SUM1+DSCALE(NN) 00002560 SUM2=SUM2+DSCALE(NN)**2 00002570 SUM3=SUM3+DAZ(NN) 00002580 SUM4=SUM4+DAZ(NN)**2 00002590 160 CONTINUE 00002600 C 00002610 C COMPUTE AND PRINT MEAN SCALE AND ROTATION 00002620 C 00002630 DSCM=SUM1/DFLOAT(NSTNSN-IFLAG) 00002640 DAZM=SUM3/DFLOAT(NSTNSN-IFLAG) 00002650 STDS=DSQRT((DFLOAT(NSTNSN-IFLAG)*SUM2-SUM1**2)/DFLOAT(NSTNSN-IFL00002660 1 AG)/DFLOAT(NSTNSN-IFLAG-1)) 00002670 STDA=DSQRT((DFLOAT(NSTNSN-IFLAG)*SUM4-SUM3**2)/DFLOAT(NSTNSN-IFL00002680 1 AG)/DFLOAT(NSTNSN-IFLAG-1)) 00002690 CALL RADARC(DAZM,I,J,S1) 00002700 CALL RADARC(STDA,I,J,S2) 00002710 PRINT 1110,DSCM,STDS,S1,S2 00002720 C 00002730 C GO COMPUTE NEXT SUBNETWORK 00002740 C 00002750 GO TO 90 00002760 C 00002770 C EXITS -996- FIXED POINT NOT IN NETWORK 00002780 C -997- SUBNETWORK POINT NOT IN NETWORK 00002790 C -998- NORMAL ENDING 00002800 C -999- STATION NUMBERS IN INPUT FILES DO NOT CORRESPOND 00002810 C 00002820 996 PRINT 1135,IFIX 00002830 STOP 00002840 997 PRINT 1140,ISTNS(I) 00002850 STOP 00002860 998 PRINT 1150 00002870 STOP 00002880 999 PRINT 1160 00002890 STOP 00002900 C 00002910 C FORMATS 00002920 C 00002930 1000 FORMAT(20A4) 00002940 1010 FORMAT('1'///20X,20A4////) 00002950 1020 FORMAT(I5,I3,2I2) 00002960 1030 FORMAT(I9,1X,12A4) 00002970 1040 FORMAT('-',35X,'FIXED POINT'/'+',35X,11('_')//26X,12A4///35X,'STAT00002980 1ION NO. ',I9) 00002990 1050 FORMAT('1',44X,'INPUT COORDINATES'/'+',44X,17('_')// 00003000 1 ,2X,'STATION NO.',19X,5H'OLD',33X,5H'NEW'/'+',1X,'______00003010 1_ ___',19X,'_____',33X,'_____'//) 00003020 1060 FORMAT(I9,1X,2I3,F7.0,I4,I3,F7.0) 00003030 1070 FORMAT('0',2X,I9,4X,4(I4,1X,I2,1X,F7.4,4X)) 00003040 1080 FORMAT('1',35X,30HDISPLACEMENTS -- 'NEW' - 'OLD'//1X,'STATION NO.'00003050 1,4X,'DELTA PHI',4X,'DELTA LONG.',4X,'DELTA PHI',4X,'DELTA LONG.', 00003060 24X,'LENGTH',4X,'SCALE CHANGE',4X,'ROTATION'/16X,'(SECONDS)',5X, 00003070 3'(SECONDS)',5X,'(METRES)',6X,'(METRES)',5X,'(METRES)',6X,'(PPM)', 00003080 47X,'(SECONDS)'//) 00003090 1085 FORMAT('0',1X,I9,24X,15('*'),5X,'FIXED POINT',5X,15('*')) 00003100 1090 FORMAT('0',1X,I9,6X,F7.3,F15.3,F14.3,F13.3,F13.3,F13.3,F13.3) 00003110 1100 FORMAT('+','*',83X,'_______',6X,'_______',2X,'**') 00003120 1110 FORMAT('-',10X,'MEAN SCALE CHANGE =',F8.3,3X,'+',2X,F5.1,2X,'PPM'/00003130 1 '+',T42,'_'//11X,'MEAN ROTATION =',F8.3,3X,'+',2X,F5.1,2X,'00003140 2 ARC-SECONDS'/'+',T38,'_') 00003150 1120 FORMAT(I9) 00003160 1130 FORMAT('1'///////////////35X,'SUB NETWORK NUMBER ',I2) 00003170 1135 FORMAT('1',15X,'***** ERROR *****',5X,'FIXED POINT ',I9,' IS NOT I00003180 1N NETWORK'/'1') 00003190 1140 FORMAT('1',15X,'***** ERROR ***** SUBNETWORK POINT ',I10,' IS NOT 00003200 1IN THE MAIN NETWORK'/'1') 00003210 1150 FORMAT('1') 00003220 1160 FORMAT('1',15X,'***** ERROR ***** THE STATION NUMBERS ON THE TWO 00003230 1INPUT FILES DO NOT CORRESPOND'/'1') 00003240 END 00003250 C***********************************************************************00003260 C* *00003270 C* S U B R O U T I N E A R C R A D *00003280 C* *00003290 C***********************************************************************00003300 C 00003310 C 00003320 C SUBROUTINE 'ARCRAD' CONVERTS DEGREES MINUTES AND SECONDS TO 00003330 C RADIANS. FOR NEGATIVE VALUES OF THE ANGLE ONLY THE LEFTMOST NON-ZERO 00003340 C VALUE IS NEGATIVE. (EGS. -30,15,30.0;0,-25,15.5;0,0,-37.2) 00003350 C 00003360 C INPUT: I = DEGREES (INTEGER) 00003370 C J = MINUTES (INTEGER) 00003380 C S = SECONDS (REAL*4) 00003390 C 00003400 C OUTPUT: RADS = ANGLE IN RADIANS (REAL*8) 00003410 C 00003420 C 00003430 SUBROUTINE ARCRAD(I,J,S,RADS) 00003440 REAL*8 RADS,DABS,DFLOAT,DBLE,SIGN,RH0/57.29577951308233D0/ 00003450 C 00003460 C CHECK FOR POSITIVE OR NEGATIVE ANGLES 00003470 C 00003480 SIGN=1.0D0 00003490 IF(I.LT.0.OR.J.LT.0.OR.S.LT.0.0)SIGN=-1.0D0 00003500 C 00003510 C CHECK FOR ANGLE OF ZERO -- IF ZERO SET RADS EXACTLY 0.0D0 00003520 C 00003530 IF(I.EQ.0.AND.J.EQ.0.AND.S.EQ.0.0)GO TO 10 00003540 C 00003550 C COMPUTE RADIAN VALUE 00003560 C 00003570 RADS=((DABS(DFLOAT(I))+DABS(DFLOAT(J)/60.0D0)+DABS(DBLE(S)/3600.0D00003580 10))/RH0)*SIGN 00003590 RETURN 00003600 10 RADS=0.0D0 00003610 RETURN 00003620 END 00003630 C 00003640 C SUBROUTINE 'RADARC' CONVERTS RADIANS TO DEGREES MINUTES AND 00003650 C SECONDS. FOR NEGATIVE ANGLES ONLY THE LEFTMOST NONZERO VALUE IS 00003660 C NEGATIVE (EGS. -50,15,30.5 ; 0,-35,30.0 ; 0,0,-50.5) 00003670 C 00003680 C NOTE: THE 0.0005 VALUE IS TO GUARD AGAINST ROUNDOFF 00003690 C 00003700 C INPUT: A = RADIAN VALUE OF ANGLE (REAL*8) 00003710 C 00003720 C OUTPUT: I = DEGREES (INTEGER) 00003730 C J = MINUTES (INTEGER) 00003740 C S = SECONDS (REAL*4) 00003750 C 00003760 SUBROUTINE RADARC(A,I,J,S) 00003770 DOUBLE PRECISION A,SEC,AD,AJ,RHO, SIGN 00003780 DATA RHO/206264.8062470963D0/ 00003790 C 00003800 C CHECK SIGN OF 'A' -- SET SIGN=-1 IF NEGATIVE AND CONVERT 'A' TO 00003810 C POSITIVE VALUE 00003820 C 00003830 SIGN=1.0D0 00003840 IF(A.LT.0.0)SIGN=-1.0D0 00003850 IF(SIGN.LT.0.0)A=-A 00003860 C 00003870 C CONVERT 'A' TO ARCSECONDS 00003880 C 00003890 SEC=A*RHO+0.0005D0 00003900 C 00003910 C FIND INTEGER DEGREES 00003920 C 00003930 I=SEC/3600.0D0 00003940 AD=I 00003950 C 00003960 C FIND INTEGER MINUTES 00003970 C 00003980 J=SEC/60.0D0-AD*60.0D0 00003990 AJ=J 00004000 C 00004010 C FIND REAL*4 SECONDS 00004020 C 00004030 S=SEC-AD*3600.0D0-AJ*60.0D0-0.0005D0 00004040 C 00004050 C SET LEFTMOST VALUE NEGATIVE IF SIGN=-1 00004060 C 00004070 IF(I.NE.0)GO TO 20 00004080 IF(J.EQ.0)GO TO 10 00004090 J=J*SIGN 00004100 GO TO 30 00004110 10 S=S*SIGN 00004120 GO TO 30 00004130 20 I=I*SIGN 00004140 C 00004150 C CONVERT 'A' BACK TO NEGATIVE IF SIGN=-1 00004160 C 00004170 30 IF(SIGN.LT.0.0)A=-A 00004180 RETURN 00004190 END 00004200 SUBROUTINE VININV(A,B1,PHI1,ONG1,PHI2,ONG2,A12,A21,S12) 00004210 C 00004220 C SUBROUTINE 'VININV' INVERSE SOLUTION OF GEODESICS ON THE ELLIPSOID 00004230 C BETWEEN POINTS 1 AND 2 00004240 C 00004250 C REFERENCE: 'DIRECT AND INVERSE SOLUTIONS OF GEODESICS ON THE 00004260 C ELLIPSOID WITH APPLICATION OF NESTED EQUATIONS' BY T. VINCENTY 00004270 C SURVEY REVIEW # 176 APRIL 1975 PAGES 88 TO 93 00004280 C 00004290 C INPUT: A=SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID 00004300 C B1=SEMI-MINOR AXIS OF REFERENCE ELLIPSOID *OR* RECIPROCAL O00004310 C THE FLATTENING 00004320 C PHI1=GEODETIC LATITUDE OF POINT 1 00004330 C ONG1=GEODETIC LONGITUDE OF POINT 1 00004340 C PHI2=GEODETIC LATITUDE OF POINT 2 00004350 C ONG2=GEODETIC LONGITUDE OF POINT 2 00004360 C OUTPUT: A12=GEODESIC AZIMUTH POINT 1 TO 2 00004370 C A21=GEODESIC AZIMUTH POINT 2 TO 1 00004380 C S12=GEODESIC DISTANCE POINT 1 TO 2 00004390 C 00004400 C NOTES: ALL LENGTH MEASURES ARE IN METRES 00004410 C ALL ANGULAR QUANTITIES ARE IN RADIANS 00004420 C LATITUDES POSITIVE NORTH 00004430 C LONGITUDES POSITIVE EAST 00004440 C AZIMUTHS POSITIVE CLOCKWISE FROM NORTH 00004450 C ITERATION TO 1*10-12 RADIANS 00004460 C ACCURACY OF DISTANCES GREATER THAN 0.01 MM AT ANY DISTANCE 00004470 C ACCURACY OF AZIMUTHS GREATER THAN 1*10-5 ARC-SECONDS 00004480 C ROUTINE CANNOT HANDLE ANTIPODAL SOLUTION - MESSAGE PRINTED 00004490 C AND ROUTINE STOPS IF ANTIPODAL SOLUTION FOUND 00004500 C FOR EQUIVALENT DIRECT SOLUTION SEE SUBROUTINE 'VINDIR' 00004510 C 00004520 C C. CHAMBERLAIN MARCH 1976 00004530 C 00004540 IMPLICIT REAL*8(A-H,O-Z) 00004550 DATA PI/3.141592653589793D0/,SING/1.0D-12/ 00004560 C 00004570 C ELLIPSOID PARAMETERS (ALSO VALID FOR A SPHERE OF RADIUS A) 00004580 C 00004590 IF(B1.EQ.0.0)B1=A 00004600 IF(B1.NE.0.0.AND.B1.LT.1000.0)F=1.0D0/B1 00004610 IF(B1.NE.0.0.AND.B1.LT.1000.0)B=A-A*F 00004620 IF(B1.GT.1000.0)B=B1 00004630 IF(B1.GT.1000.0)F=(A-B)/A 00004640 IF(B.EQ.A)B1=0.0D0 00004650 B2=B*B 00004660 OMF=1.0D0-F 00004670 C 00004680 C REDUCED LATITUDE OF POINTS 1 AND 2 AND THEIR TRIG FUNCTIONS AND 00004690 C COMBINATIONS 00004700 C 00004710 TPI=2.0D0*PI 00004720 POT=PI/2.0D0 00004730 U1=POT 00004740 IF(PHI1.LT.0.0)U1=-U1 00004750 IF(DABS(DABS(PHI1)-POT).GT.SING)U1=DATAN(OMF*DTAN(PHI1)) 00004760 SU1=DSIN(U1) 00004770 CU1=DCOS(U1) 00004780 U2=POT 00004790 IF(PHI2.LT.0.0)U2=-U2 00004800 IF(DABS(DABS(PHI2)-POT).GT.SING)U2=DATAN(OMF*DTAN(PHI2)) 00004810 SU2=DSIN(U2) 00004820 CU2=DCOS(U2) 00004830 CU1SU2=CU1*SU2 00004840 SU1CU2=SU1*CU2 00004850 SU1SU2=SU1*SU2 00004860 CU1CU2=CU1*CU2 00004870 C 00004880 C CHANGE IN LONGITUDE ON AUXILLARY SPHERE ITERATION TO 1*10-12 RADIANS00004890 C (EQUATIONS 13,14,15,16,17,18,10,11) 00004900 C 00004910 ITER=0 00004920 DLON=ONG2-ONG1 00004930 IF(DABS(DLON-PI).LT.SING.AND.DABS(PHI1).LT.SING.AND.DABS(PHI2).LT.00004940 1SING)GO TO 30 00004950 DLAS=DLON 00004960 10 SDLAS=DSIN(DLAS) 00004970 CDLAS=DCOS(DLAS) 00004980 SSIG=DSQRT((CU2*SDLAS)**2+(CU1SU2-SU1CU2*CDLAS)**2) 00004990 CSIG=SU1SU2+CU1CU2*CDLAS 00005000 TSIG=SSIG/CSIG 00005010 SIG=DARCOS(CSIG) 00005020 SALPHA=CU1CU2*SDLAS/SSIG 00005030 CALPH2=1.0D0-SALPHA**2 00005040 CTSM=0.0D0 00005050 IF(DABS(CALPH2).GT.SING )CTSM=CSIG-2.0D0*SU1SU2/CALPH2 00005060 CTSM2=CTSM**2 00005070 C=F/16.0D0*CALPH2*(4.0D0+F*(4.0D0-3.0D0*CALPH2)) 00005080 DLASUP=DLON+(1.0D0-C)*F*SALPHA*(SIG+C*SSIG*(CTSM+C*CSIG*00005090 1 (-1.0D0+2.0D0*CTSM2))) 00005100 IF(DABS(DLASUP-DLAS).LT.SING )GO TO 20 00005110 DLAS=DLASUP 00005120 ITER=ITER+1 00005130 IF(ITER.GT.50)GO TO 30 00005140 GO TO 10 00005150 C 00005160 C GEODESIC DISTANCE (EQUATIONS 3,4,6,19) 00005170 C 00005180 20 USQ=CALPH2*(A*A-B2)/B2 00005190 A3=1.0D0+USQ/16384.0D0*(4096.0D0+USQ*(-768.0D0+USQ*(320.0D0- 00005200 1 175.0D0*USQ))) 00005210 B4=USQ/1024.0D0*(256.0D0+USQ*(-128.0D0+USQ*(74.0D0-47.0D0*USQ)))00005220 DELSIG=B4*SSIG*(CTSM+B4/4.0D0*(CSIG*(-1.0D0+2.0D0*CTSM2)-B4/ 00005230 1 6.0D0*CTSM*(-3.0D0+4.0D0*SSIG**2)*(-3.0D0+4.0D0*CTSM2))) 00005240 S12=B*A3*(SIG-DELSIG) 00005250 C 00005260 C GEODESIC AZIMUTHS 1 TO 2 AND 2 TO 1 (EQUATIONS 20 AND 21) 00005270 C 00005280 A12=DATAN2(CU2*SDLAS,CU1SU2-SU1CU2*CDLAS) 00005290 A21=DATAN2(-CU1*SDLAS,SU1CU2-CU1SU2*CDLAS) 00005300 IF(A12.GT.TPI)A12=A12-TPI 00005310 IF(A12.LT.0.0)A12=A12+TPI 00005320 IF(A21.GT.TPI)A21=A21-TPI 00005330 IF(A21.LT.0.0)A21=A21+TPI 00005340 RETURN 00005350 C 00005360 C ERROR EXIT FOR ANTIPODAL SOLUTION 00005370 C 00005380 30 PRINT 1000 00005390 STOP 00005400 1000 FORMAT('1'//////////15X,'SUBROUTINE VININV WILL NOT HANDLE ANTIPOD00005410 1AL SOLUTION'/'1') 00005420 END 00005430