C***********************************************************************00000040 C* *00000050 C* *00000060 C* *************************** *00000070 C* * * *00000080 C* * C O M P A R E 1 * *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* PROGRAM TAKES AS INPUT TWO SETS OF ADJUSTED COORDINATES AND A LIST *00000160 C* OF SEQUENCE NUMBERS THAT SPECIFY ALONG WHICH LINES OBSERVATIONS *00000170 C* HAVE BEEN MADE. THE PROGRAM COMPUTES A SCALE CHANGE AND AZIMUTH *00000180 C* ROTATION ALONG EACH LINE THAT HAS BEEN OBSERVED (EACH LINE ONCE). *00000190 C* *00000200 C* INPUT IS AS FOLLOWS: *00000210 C* *00000220 C* CARD 1: FORMAT (20A4); JOB DESCRIPTION *00000230 C* *00000240 C* CARD 2: FORMAT(I5,3I2) ; (1) NUMBER OF STATIONS IN NETWORK *00000250 C* (2) UNIT NUMBER FOR 'OLD' COORDINATES *00000260 C* (3) UNIT NUMBER FOR 'NEW' COORDINATES *00000270 C* (4) UNIT NUMBER FOR 'FROM-TO' DATA *00000280 C* *00000290 C* -STATION NUMBERS AND ADJUSTED COORDINATES CAN BE READ FROM *00000300 C* AN AUXILLARY DEVICE (DISK OR TAPE) OR FROM CARDS, FOR CARDS *00000310 C* UNIT NUMBER MUST BE 5, FOR AUXILLARY DEVICE UNIT NUMBER IS *00000320 C* ANY VALID FORTRAN UNIT NUMBER WITH THE APPROPRIATE 'DD' *00000330 C* STATEMENT. *00000340 C* -'FROM-TO' DATA CAN ALSO BE READ FROM CARDS OR AN AUXILLARY *00000350 C* DEVICE - SEE COMMENT ABOVE *00000360 C* *00000370 C* CARD 3: 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 4: SAME AS CARD 3 BUT FOR 'NEW' COORDINATES *00000420 C* *00000430 C* NOTE:-CARDS 3 & 4 ARE REQUIRED ONLY IF UNIT NUMBERS ON CARD *00000440 C* NUMBER 2 ARE 5 *00000450 C* -IF CARDS 3 & 4 ARE REQUIRED THEY ARE REPEATED ONCE FOR *00000460 C* EACH NETWORK POINT *00000470 C* -IF INPUT IS NOT FROM CARDS (CARD TYPES 3 & 4) 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 (PROGRAM *00000520 C* 'SOLVE' CREATES THIS DATA SET) *00000530 C* *00000540 C* CARD 5: FORMAT(I5,I3,20(I3)); SEQUENCE NUMBER OF FROM STATION, *00000550 C* NUMBER OF TO STATIONS, SEQUENCE NUMBERS OF TO *00000560 C* STATIONS *00000570 C* NOTE:-IF INPUT FOR 'FROM-TO' IS FROM CARDS CARD 5 IS REPEATED *00000580 C* ONCE FOR EACH FROM STATION *00000590 C* -IF INPUT FOR 'FROM-TO' IS FROM DISC(TAPE) CREATE A *00000600 C* SEQUENTIAL FILE AS FOLLOWS-> WRITE(IUNIT)IFROM,NTO,(ISTO(I)*00000610 C* ,I=1,NTO) *00000620 C* -IF INPUT IS NOT FROM CARDS CARD 5 IS NOT REQUIRED *00000630 C* -IF INPUT IS FROM CARDS THE PROGRAM COMPUTES THOSE LINES *00000640 C* SPECIFIED, IF INPUT IS FROM DISC(TAPE) (OUTPUT FROM FORM3) *00000650 C* THE PROGRAM CHECKS THE NUMBERS AND ELIMINATES REFERENCES *00000660 C* THAT WOULD CAUSE A COMPUTATION FROM J TO I IF I TO J HAS *00000670 C* ALREADY BEEN COMPUTED *00000680 C* *00000690 C* GENERAL NOTE: PROGRAM COMPUTES 'NEW - OLD' WHICH ARE DISTINGUISHED *00000700 C* ONLY BY THE ORDER IN WHICH THE COORDINATES ARE READ IN. *00000710 C* THE PROGRAM IS DIMENSIONED TO HANDLE A 1000 STATION *00000720 C* NETWORK WITH 6000 OBSERVATIONS, *00000730 C* FOR LARGER NETWORKS THE DIMENSIONS OF THE APPROPRIATE DIMENSIONED *00000740 C* VARIABLES MUST BE INCREASED *00000750 C* *00000760 C* *00000770 C* C. CHAMBERLAIN MAY /1976 *00000780 C* *00000790 C***********************************************************************00000800 IMPLICIT REAL*8(A-H,O-Z) 00000810 REAL*4 S1,S2,S3,S4,ABS 00000820 DIMENSION DJOB(20),NTOV(30),ISTNO(1000),ISTNN(1000),PHIO(1000), 00000830 1PHIN(1000),XLONO(1000),XLONN(1000),IFROMV(6000),ITOV(6000) 00000840 C 00000850 C CLARKE 1866 ELLIPSOID 00000860 C 00000870 A=6378206.4D0 00000880 B=6356583.8D0 00000890 ESQ=(A*A-B*B)/A/A 00000900 F=(A-B)/A 00000910 ESQP=ESQ/(1.0D0-ESQ) 00000920 C 00000930 C READ AND PRINT CARDS 1,2 00000940 C 00000950 READ 1000,DJOB 00000960 PRINT 1010,DJOB 00000970 READ 1020,NSTN,IUO,IUN,IUFT 00000980 C 00000990 C READ IN STATION NUMBERS AND COORDINATES 00001000 C 00001010 DO 50 I=1,NSTN 00001020 IF(MOD(I,26).EQ.1)PRINT 1030 00001030 IF(IUO.EQ.5)READ 1040,ISTNO(I),K,L,S1,KK,LL,S2 00001040 IF(IUO.NE.5)READ(IUO)ISTNO(I),PHIO(I),XLONO(I) 00001050 IF(IUO.EQ.5)GO TO 10 00001060 CALL RADARC(PHIO(I),K,L,S1) 00001070 CALL RADARC(XLONO(I),KK,LL,S2) 00001080 GO TO 20 00001090 10 CALL ARCRAD(K,L,S1,PHIO(I)) 00001100 CALL ARCRAD(KK,LL,S2,XLONO(I)) 00001110 20 IF(IUN.EQ.5)READ 1040,ISTNN(I),J,M,S3,JJ,MM,S4 00001120 IF(IUN.NE.5)READ(IUN)ISTNN(I),PHIN(I),XLONN(I) 00001130 IF(IUN.EQ.5)GO TO 30 00001140 CALL RADARC(PHIN(I),J,M,S3) 00001150 CALL RADARC(XLONN(I),JJ,MM,S4) 00001160 GO TO 40 00001170 30 CALL ARCRAD(J,M,S3,PHIN(I)) 00001180 CALL ARCRAD(JJ,MM,S4,XLONN(I)) 00001190 40 IF(ISTNO(I).NE.ISTNN(I))GO TO 999 00001200 50 PRINT 1050,ISTNO(I),K,L,S1,KK,LL,S2,J,M,S3,JJ,MM,S4 00001210 C 00001220 C READ IN 'FROM-TO' DATA AND STORE IN VECTORS 00001230 C 00001240 IC=1 00001250 60 CONTINUE 00001260 IF(IUFT.EQ.5)READ(5,1060,END=130)IFROM,NTO,(NTOV(K),K=1,NTO) 00001270 IF(IUFT.EQ.5)GO TO 70 00001280 READ(IUFT,END=90)IFROM,NTO,(NTOV(K),K=1,NTO) 00001290 70 ICT=IC+NTO-1 00001300 K=1 00001310 DO 80 J=IC,ICT 00001320 IFROMV(J)=IFROM 00001330 ITOV(J)=NTOV(K) 00001340 80 K=K+1 00001350 IC=ICT+1 00001360 GO TO 60 00001370 90 IC=IC-1 00001380 C 00001390 C IF INPUT IS FROM DISK(TAPE) ELIMINATE REFERENCES TO LINES OBSERVED 00001400 C MORE THAN ONCE 00001410 C 00001420 DO 120 I=1,IC 00001430 ICHECK=IFROMV(I) 00001440 JCHECK=ITOV(I) 00001450 IF(ICHECK.EQ.0)GO TO 120 00001460 DO 110 K=1,IC 00001470 IF(K.EQ.I)GO TO 110 00001480 IF(ITOV(K).NE.ICHECK.OR.IFROMV(K).NE.JCHECK)GO TO 100 00001490 IFROMV(K)=0 00001500 100 IF(ITOV(K).NE.JCHECK.OR.IFROMV(K).NE.ICHECK)GO TO 110 00001510 IFROMV(K)=0 00001520 110 CONTINUE 00001530 120 CONTINUE 00001540 130 IC=IC-1 00001550 C 00001560 C COMPUTE CHANGES IN COORDINATES AT EACH STATION 00001570 C 00001580 DO 140 I=1,NSTN 00001590 C 00001600 C COORDINATE DIFFERENCES IN RADIANS (NEW - OLD) 00001610 C 00001620 DPHIR=PHIN(I)-PHIO(I) 00001630 DXLONR=XLONN(I)-XLONO(I) 00001640 C 00001650 C COORDINATE DIFFERENCES IN METRES AND LENGTH OF 'MISCLOSURE' VECTOR 00001660 C 00001670 DEN=DSQRT(1.0D0-ESQ*DSIN(PHIO(I))**2) 00001680 RN=A/DEN 00001690 RM=A*(1.0D0-ESQ)/DEN**3 00001700 DPHIM=DPHIR*RM 00001710 DXLONM=DXLONR*RN*DCOS(PHIO(I)) 00001720 VEC=DSQRT(DPHIM**2+DXLONM**2) 00001730 CALL RADARC(DPHIR ,J,K,S2) 00001740 CALL RADARC(DXLONR ,J,K,S3) 00001750 IF(MOD(I,26).EQ.1)PRINT 1070 00001760 PRINT 1080,ISTNN(I),S2,S3,DPHIM,DXLONM,VEC 00001770 140 CONTINUE 00001780 C 00001790 C COMPUTE CHANGE IN SCALE AND AZIMUTH ALONG EACH OBSERVED LINE 00001800 C FLAG SCALE CHANGES > 10 PPM AND ROTATIONS > 1SEC 00001810 C 00001820 SUM1=0.0D0 00001830 SUM2=0.0D0 00001840 SUM3=0.0D0 00001850 SUM4=0.0D0 00001860 KOUNT=0 00001870 DO 150 I=1,IC 00001880 IF(IFROMV(I).EQ.0)GO TO 150 00001890 KOUNT=KOUNT+1 00001900 J=IFROMV(I) 00001910 K=ITOV(I) 00001920 CALL VININV(A,B,PHIN(J),XLONN(J),PHIN(K),XLONN(K),AZN,AZ,DISN) 00001930 CALL VININV(A,B,PHIO(J),XLONO(J),PHIO(K),XLONO(K),AZO,AZ,DISO) 00001940 DSCALE=(DISN-DISO)/DISN*1.0D6 00001950 DAZ=AZN-AZO 00001960 CALL RADARC(DAZ,KK,LL,S1) 00001970 DLEN=DISN-DISO 00001980 SUM1=SUM1+DSCALE**2 00001990 SUM2=SUM2+DSCALE 00002000 SUM3=SUM3+DAZ**2 00002010 SUM4=SUM4+DAZ 00002020 IF(MOD(KOUNT,26).EQ.1)PRINT 1090 00002030 PRINT 1100,ISTNN(J),ISTNN(K),S1 ,DLEN,DSCALE 00002040 IF(DABS(DSCALE).GT.10.0D0.OR.ABS(S1).GT.1.0)PRINT 1110 00002050 150 CONTINUE 00002060 C 00002070 C COMPUTE MEAN AND STANDARD DEVIATION OF SCALE AND ROTATION CHANGE 00002080 C 00002090 FLN=DFLOAT(KOUNT) 00002100 DSCM=SUM2/FLN 00002110 STS=DSQRT((FLN*SUM1-SUM2**2)/FLN/(FLN-1.0D0)) 00002120 DAZM=SUM4/FLN 00002130 STA=DSQRT((FLN*SUM3-SUM4**2)/FLN/(FLN-1.0D0)) 00002140 CALL RADARC(DAZM,I,J,S1) 00002150 CALL RADARC(STA,I,J,S2) 00002160 PRINT 1120,DSCM,STS,S1,S2,KOUNT 00002170 C 00002180 C EXITS -998- NORMAL ENDING 00002190 C -999- STATION NUMBERS IN INPUT FILES DO NOT CORRESPOND 00002200 C 00002210 998 PRINT 1130 00002220 STOP 00002230 999 PRINT 1140 00002240 STOP 00002250 C 00002260 C FORMATS 00002270 C 00002280 1000 FORMAT(20A4) 00002290 1010 FORMAT('1'///20X,20A4////) 00002300 1020 FORMAT(I5,3I2) 00002310 1030 FORMAT('1',44X,'INPUT COORDINATES'/'+',44X,17('_')// 00002320 1 ,2X,'STATION NO.',19X,5H'OLD',33X,5H'NEW'/'+',1X,'______00002330 1_ ___',19X,'_____',33X,'_____'//) 00002340 1040 FORMAT(I9,1X,2I3,F7.0,I4,I3,F7.0) 00002350 1050 FORMAT('0',2X,I9,4X,4(I4,1X,I2,1X,F7.4,4X)) 00002360 1060 FORMAT(I5,I3,20(I3)) 00002370 1070 FORMAT('1',35X,30HDISPLACEMENTS -- 'NEW' - 'OLD'//11X,'STATION NO.00002380 1',4X,'DELTA PHI',4X,'DELTA LONG.',4X,'DELTA PHI',4X,'DELTA LONG.',00002390 24X,'LENGTH' /26X,'(SECONDS)',5X, 00002400 3'(SECONDS)',5X,'(METRES)',6X,'(METRES)',5X,'(METRES)'//) 00002410 1080 FORMAT('0',11X,I9,6X,F7.3,F15.3,F14.3,F13.3,F13.3) 00002420 1090 FORMAT('1',17X,'LINE',14X,'CHANGE IN',8X,'CHANGE IN LENGTH'/36X, 00002430 1'DIRECTION'/34X,'(ARC-SECONDS) METRES PPM') 00002440 1100 FORMAT('0',5X,I9,5X,I9,9X,F5.1,11X,F5.1,4X,F5.1) 00002450 1110 FORMAT('+','**',69X,'**') 00002460 1120 FORMAT('-',10X,'MEAN SCALE CHANGE =',F8.2,2X,'+',F5.1,2X,'PPM'/'+'00002470 1,39X,'_'//14X,'MEAN DIRECTION =',F8.2,2X,'+',F5.1,2X,'ARC-SECONDS'00002480 2/'+',39X,'_'////10X,'BASED ON', I7,3X,'DETERMINATIONS') 00002490 1130 FORMAT('1') 00002500 1140 FORMAT('1',15X,'***** ERROR ***** THE STATION NUMBERS ON THE TWO 00002510 1INPUT FILES DO NOT CORRESPOND'/'1') 00002520 END 00002530 SUBROUTINE VININV(A,B1,PHI1,ONG1,PHI2,ONG2,A12,A21,S12) 00002540 C 00002550 C SUBROUTINE 'VININV' INVERSE SOLUTION OF GEODESICS ON THE ELLIPSOID 00002560 C BETWEEN POINTS 1 AND 2 00002570 C 00002580 C REFERENCE: 'DIRECT AND INVERSE SOLUTIONS OF GEODESICS ON THE 00002590 C ELLIPSOID WITH APPLICATION OF NESTED EQUATIONS' BY T. VINCENTY 00002600 C SURVEY REVIEW # 176 APRIL 1975 PAGES 88 TO 93 00002610 C 00002620 C INPUT: A=SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID 00002630 C B1=SEMI-MINOR AXIS OF REFERENCE ELLIPSOID *OR* RECIPROCAL O00002640 C THE FLATTENING 00002650 C PHI1=GEODETIC LATITUDE OF POINT 1 00002660 C ONG1=GEODETIC LONGITUDE OF POINT 1 00002670 C PHI2=GEODETIC LATITUDE OF POINT 2 00002680 C ONG2=GEODETIC LONGITUDE OF POINT 2 00002690 C OUTPUT: A12=GEODESIC AZIMUTH POINT 1 TO 2 00002700 C A21=GEODESIC AZIMUTH POINT 2 TO 1 00002710 C S12=GEODESIC DISTANCE POINT 1 TO 2 00002720 C 00002730 C NOTES: ALL LENGTH MEASURES ARE IN METRES 00002740 C ALL ANGULAR QUANTITIES ARE IN RADIANS 00002750 C LATITUDES POSITIVE NORTH 00002760 C LONGITUDES POSITIVE EAST 00002770 C AZIMUTHS POSITIVE CLOCKWISE FROM NORTH 00002780 C ITERATION TO 1*10-12 RADIANS 00002790 C ACCURACY OF DISTANCES GREATER THAN 0.01 MM AT ANY DISTANCE 00002800 C ACCURACY OF AZIMUTHS GREATER THAN 1*10-5 ARC-SECONDS 00002810 C ROUTINE CANNOT HANDLE ANTIPODAL SOLUTION - MESSAGE PRINTED 00002820 C AND ROUTINE STOPS IF ANTIPODAL SOLUTION FOUND 00002830 C FOR EQUIVALENT DIRECT SOLUTION SEE SUBROUTINE 'VINDIR' 00002840 C 00002850 C C. CHAMBERLAIN MARCH 1976 00002860 C 00002870 IMPLICIT REAL*8(A-H,O-Z) 00002880 DATA PI/3.141592653589793D0/,SING/1.0D-12/ 00002890 C 00002900 C ELLIPSOID PARAMETERS (ALSO VALID FOR A SPHERE OF RADIUS A) 00002910 C 00002920 IF(B1.EQ.0.0)B1=A 00002930 IF(B1.NE.0.0.AND.B1.LT.1000.0)F=1.0D0/B1 00002940 IF(B1.NE.0.0.AND.B1.LT.1000.0)B=A-A*F 00002950 IF(B1.GT.1000.0)B=B1 00002960 IF(B1.GT.1000.0)F=(A-B)/A 00002970 IF(B.EQ.A)B1=0.0D0 00002980 B2=B*B 00002990 OMF=1.0D0-F 00003000 C 00003010 C REDUCED LATITUDE OF POINTS 1 AND 2 AND THEIR TRIG FUNCTIONS AND 00003020 C COMBINATIONS 00003030 C 00003040 TPI=2.0D0*PI 00003050 POT=PI/2.0D0 00003060 U1=POT 00003070 IF(PHI1.LT.0.0)U1=-U1 00003080 IF(DABS(DABS(PHI1)-POT).GT.SING)U1=DATAN(OMF*DTAN(PHI1)) 00003090 SU1=DSIN(U1) 00003100 CU1=DCOS(U1) 00003110 U2=POT 00003120 IF(PHI2.LT.0.0)U2=-U2 00003130 IF(DABS(DABS(PHI2)-POT).GT.SING)U2=DATAN(OMF*DTAN(PHI2)) 00003140 SU2=DSIN(U2) 00003150 CU2=DCOS(U2) 00003160 CU1SU2=CU1*SU2 00003170 SU1CU2=SU1*CU2 00003180 SU1SU2=SU1*SU2 00003190 CU1CU2=CU1*CU2 00003200 C 00003210 C CHANGE IN LONGITUDE ON AUXILLARY SPHERE ITERATION TO 1*10-12 RADIANS00003220 C (EQUATIONS 13,14,15,16,17,18,10,11) 00003230 C 00003240 ITER=0 00003250 DLON=ONG2-ONG1 00003260 IF(DABS(DLON-PI).LT.SING.AND.DABS(PHI1).LT.SING.AND.DABS(PHI2).LT.00003270 1SING)GO TO 30 00003280 DLAS=DLON 00003290 10 SDLAS=DSIN(DLAS) 00003300 CDLAS=DCOS(DLAS) 00003310 SSIG=DSQRT((CU2*SDLAS)**2+(CU1SU2-SU1CU2*CDLAS)**2) 00003320 CSIG=SU1SU2+CU1CU2*CDLAS 00003330 TSIG=SSIG/CSIG 00003340 SIG=DARCOS(CSIG) 00003350 SALPHA=CU1CU2*SDLAS/SSIG 00003360 CALPH2=1.0D0-SALPHA**2 00003370 CTSM=0.0D0 00003380 IF(DABS(CALPH2).GT.SING )CTSM=CSIG-2.0D0*SU1SU2/CALPH2 00003390 CTSM2=CTSM**2 00003400 C=F/16.0D0*CALPH2*(4.0D0+F*(4.0D0-3.0D0*CALPH2)) 00003410 DLASUP=DLON+(1.0D0-C)*F*SALPHA*(SIG+C*SSIG*(CTSM+C*CSIG*00003420 1 (-1.0D0+2.0D0*CTSM2))) 00003430 IF(DABS(DLASUP-DLAS).LT.SING )GO TO 20 00003440 DLAS=DLASUP 00003450 ITER=ITER+1 00003460 IF(ITER.GT.50)GO TO 30 00003470 GO TO 10 00003480 C 00003490 C GEODESIC DISTANCE (EQUATIONS 3,4,6,19) 00003500 C 00003510 20 USQ=CALPH2*(A*A-B2)/B2 00003520 A3=1.0D0+USQ/16384.0D0*(4096.0D0+USQ*(-768.0D0+USQ*(320.0D0- 00003530 1 175.0D0*USQ))) 00003540 B4=USQ/1024.0D0*(256.0D0+USQ*(-128.0D0+USQ*(74.0D0-47.0D0*USQ)))00003550 DELSIG=B4*SSIG*(CTSM+B4/4.0D0*(CSIG*(-1.0D0+2.0D0*CTSM2)-B4/ 00003560 1 6.0D0*CTSM*(-3.0D0+4.0D0*SSIG**2)*(-3.0D0+4.0D0*CTSM2))) 00003570 S12=B*A3*(SIG-DELSIG) 00003580 C 00003590 C GEODESIC AZIMUTHS 1 TO 2 AND 2 TO 1 (EQUATIONS 20 AND 21) 00003600 C 00003610 A12=DATAN2(CU2*SDLAS,CU1SU2-SU1CU2*CDLAS) 00003620 A21=DATAN2(-CU1*SDLAS,SU1CU2-CU1SU2*CDLAS) 00003630 IF(A12.GT.TPI)A12=A12-TPI 00003640 IF(A12.LT.0.0)A12=A12+TPI 00003650 IF(A21.GT.TPI)A21=A21-TPI 00003660 IF(A21.LT.0.0)A21=A21+TPI 00003670 RETURN 00003680 C 00003690 C ERROR EXIT FOR ANTIPODAL SOLUTION 00003700 C 00003710 30 PRINT 1000 00003720 STOP 00003730 1000 FORMAT('1'//////////15X,'SUBROUTINE VININV WILL NOT HANDLE ANTIPOD00003740 1AL SOLUTION'/'1') 00003750 END 00003760 C 00003770 C SUBROUTINE 'RADARC' CONVERTS RADIANS TO DEGREES MINUTES AND 00003780 C SECONDS. FOR NEGATIVE ANGLES ONLY THE LEFTMOST NONZERO VALUE IS 00003790 C NEGATIVE (EGS. -50,15,30.5 ; 0,-35,30.0 ; 0,0,-50.5) 00003800 C 00003810 C NOTE: THE 0.0005 VALUE IS TO GUARD AGAINST ROUNDOFF 00003820 C 00003830 C INPUT: A = RADIAN VALUE OF ANGLE (REAL*8) 00003840 C 00003850 C OUTPUT: I = DEGREES (INTEGER) 00003860 C J = MINUTES (INTEGER) 00003870 C S = SECONDS (REAL*4) 00003880 C 00003890 SUBROUTINE RADARC(A,I,J,S) 00003900 DOUBLE PRECISION A,SEC,AD,AJ,RHO, SIGN 00003910 DATA RHO/206264.8062470963D0/ 00003920 C 00003930 C CHECK SIGN OF 'A' -- SET SIGN=-1 IF NEGATIVE AND CONVERT 'A' TO 00003940 C POSITIVE VALUE 00003950 C 00003960 SIGN=1.0D0 00003970 IF(A.LT.0.0)SIGN=-1.0D0 00003980 IF(SIGN.LT.0.0)A=-A 00003990 C 00004000 C CONVERT 'A' TO ARCSECONDS 00004010 C 00004020 SEC=A*RHO+0.0005D0 00004030 C 00004040 C FIND INTEGER DEGREES 00004050 C 00004060 I=SEC/3600.0D0 00004070 AD=I 00004080 C 00004090 C FIND INTEGER MINUTES 00004100 C 00004110 J=SEC/60.0D0-AD*60.0D0 00004120 AJ=J 00004130 C 00004140 C FIND REAL*4 SECONDS 00004150 C 00004160 S=SEC-AD*3600.0D0-AJ*60.0D0-0.0005D0 00004170 C 00004180 C SET LEFTMOST VALUE NEGATIVE IF SIGN=-1 00004190 C 00004200 IF(I.NE.0)GO TO 20 00004210 IF(J.EQ.0)GO TO 10 00004220 J=J*SIGN 00004230 GO TO 30 00004240 10 S=S*SIGN 00004250 GO TO 30 00004260 20 I=I*SIGN 00004270 C 00004280 C CONVERT 'A' BACK TO NEGATIVE IF SIGN=-1 00004290 C 00004300 30 IF(SIGN.LT.0.0)A=-A 00004310 RETURN 00004320 END 00004330 C***********************************************************************00004340 C* *00004350 C* S U B R O U T I N E A R C R A D *00004360 C* *00004370 C***********************************************************************00004380 C 00004390 C 00004400 C SUBROUTINE 'ARCRAD' CONVERTS DEGREES MINUTES AND SECONDS TO 00004410 C RADIANS. FOR NEGATIVE VALUES OF THE ANGLE ONLY THE LEFTMOST NON-ZERO 00004420 C VALUE IS NEGATIVE. (EGS. -30,15,30.0;0,-25,15.5;0,0,-37.2) 00004430 C 00004440 C INPUT: I = DEGREES (INTEGER) 00004450 C J = MINUTES (INTEGER) 00004460 C S = SECONDS (REAL*4) 00004470 C 00004480 C OUTPUT: RADS = ANGLE IN RADIANS (REAL*8) 00004490 C 00004500 C 00004510 SUBROUTINE ARCRAD(I,J,S,RADS) 00004520 REAL*8 RADS,DABS,DFLOAT,DBLE,SIGN,RH0/57.29577951308233D0/ 00004530 C 00004540 C CHECK FOR POSITIVE OR NEGATIVE ANGLES 00004550 C 00004560 SIGN=1.0D0 00004570 IF(I.LT.0.OR.J.LT.0.OR.S.LT.0.0)SIGN=-1.0D0 00004580 C 00004590 C CHECK FOR ANGLE OF ZERO -- IF ZERO SET RADS EXACTLY 0.0D0 00004600 C 00004610 IF(I.EQ.0.AND.J.EQ.0.AND.S.EQ.0.0)GO TO 10 00004620 C 00004630 C COMPUTE RADIAN VALUE 00004640 C 00004650 RADS=((DABS(DFLOAT(I))+DABS(DFLOAT(J)/60.0D0)+DABS(DBLE(S)/3600.0D00004660 10))/RH0)*SIGN 00004670 RETURN 00004680 10 RADS=0.0D0 00004690 RETURN 00004700 END 00004710