C********************************************************************* C PROGRAM TO CONVERT FROM UTM COORDINATES (IN METERS) TO LATITUDE C AND LONGITUDE. SUBROUTINES WERE PROVIDED BY GARTH D. NEWTON AND C ARE DESCRIBED IN U.S. GEOLOGICAL SURVEY BULLETIN 1642, 1985. C C NOTE: THE REQUIRED USE OF DOUBLE PRECISION C WRITTEN BY WRDANSKIN@DCASCR C LAST REVISION 1/6/87 C********************************************************************* PROGRAM MAIN IMPLICIT DOUBLE PRECISION (A-Z) DIMENSION LAT(3),LON(3) C INITIALIZE DATA DEPENDENT UPON THE CHOICE OF ORIGIN, E.G. SETS ZONE LAT(1)=37. LAT(2)=0. LAT(3)=0. LON(1)=118. LON(2)=0. LON(3)=0. C CONVERT TO DECIMAL DEGREES CALL DEGREE(LAT,ALAT) CALL DEGREE(LON,ALON) C SET GLOBAL PARAMETERS CALL DATA(ALAT,ALON) C READ DATA IN UTM COORDINATES C USE A DO-LOOP OF NCOUNT X,Y PAIRS READ(5,*)NCOUNT WRITE(6,900) DO 100 I=1,NCOUNT READ(5,*)XUTM,YUTM X0=XUTM Y0=YUTM IF(XUTM.EQ.0.)GOTO 200 C CALCULATE DECIMAL LAT/LONG FROM UTM COORDINATES CALL INVERSE(ALAT,ALON,XUTM,YUTM) C CONVERT DECIMAL LAT/LONG TO DEG-MIN-SEC LAT/LONG CALL DECIMAL(LAT,ALAT) CALL DECIMAL(LON,ALON) C PRINT RESULTS WRITE(6,901)X0,Y0,ALAT,ALON,LAT,LON 100 CONTINUE 200 CONTINUE C FORMAT STATEMENTS 900 FORMAT(' XUTM YUTM LAT LONG LAT', -' LONG') 901 FORMAT(2F10.0,2F10.4,1X,3I3,1X,3I3) STOP END C******************************************************************** C CONVERT DEG-MIN-SEC TO DEGREES C SUBROUTINE DEGREE(A,B) IMPLICIT DOUBLE PRECISION (A-Z) C C--------------------------------------------------------------------- DIMENSION A(3) B=A(1)+A(2)/60.+A(3)/3600. RETURN END C********************************************************************* C CONVERT DEGREES TO DEG-MIN-SEC C SUBROUTINE DECIMAL(A,B) IMPLICIT DOUBLE PRECISION (A-Z) C C--------------------------------------------------------------------- DIMENSION A(3) A(1)=AINT(B+.5/3600.) A(2)=AINT((B-A(1))*60.+.5/60.) A(3)=AINT((B-A(1)-A(2)/60.)*3600.+.5) RETURN END C*********************************************************************** C COMPUTE THE LATITUDE AND LONGITUDE GIVEN THE C RECTANGULAR COORDINATES X AND Y C C1 CHANGED AUTMI TO INVERSE 01/09/85 SUBROUTINE INVERSE(LAT,LON,X,Y) IMPLICIT DOUBLE PRECISION (A-Z) C C--------------------------------------------------------------------- COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0 C X=(X-500000.)/M0 Y=Y/M0 C C----- ITERATE TO COMPUTE MERIDIONAL DISTANCE C MAKE INITIAL GUESS PHI=SINV(Y,PHI0) OLD=Y C BEGIN ITERATIONS DO 10 I=1,10 PHI=SINV(Y,PHI) C IF OLD VALUE EQUALS NEW VALUE STOP IF(OLD.EQ.Y0) GOTO 11 OLD=Y0 10 CONTINUE 11 S=DSIN(PHI) C=DCOS(PHI) T=S/C T2=T*T N2=EPS*C*C N=A/SQRT(1.D0-E2*S*S) AA0=X/N AA1=AA0**3*(1.D0+2.D0*T2+N2)/6.D0 AA2=AA0**5*(5.D0+28.D0*T2+24.D0*T2*T2+6.D0*N2+8.D0*T2*N2) #/120.D0 C C COMPUTE LONGITUDE (RADIANS) CC LON=LAM0-(1.D0/C)*(AA0-AA1+AA2) C T3=T*T2 T4=T*T3 T5=T*T4 T6=T*T5 N4=N2*N2 N6=N4*N2 N8=N6*N2 AA1=T*(1.D0+N2)*AA0*AA0 AA2=T*(1.D0+N2)*AA0**4*(5.D0+3.D0*T2+N2-4.D0*N4-9.D0*N2*T2) AA3=T*(1.D0+N2)*AA0**6*(61.D0+90.D0*T2+46.D0*N2+45.D0*T4 # -252.D0*T2*N2 #-3.D0*N4+100*N6-66*T2*N4-90.D0*T4*N2+88.D0*N8 #+225.D0*T4*N4+84.D0*T2*N6-192.D0*T2*N8) AA4=T*(1.D0+N2)*AA0**8*(1385.D0+3633.D0*T2+4095.D0*T4+1575.D0*T6) C C COMPUTE LATITUDE (RADIANS) C LAT=PHI-AA1/2.D0+AA2/24.D0-AA3/720.D0+AA4/40320.D0 C C CONVERT TO DEGREES LAT=LAT /D2R LON=LON /D2R C RETURN END C********************************************************************* C INITIALIZE CONSTANTS RELATED TO THE LATITUDE AND LONGITUDE C OF THE AREA BEING MAPPED C C1 CHANGE FROM UTM TO DATA 01/09/85 SUBROUTINE DATA(BASLAT,LAM) IMPLICIT DOUBLE PRECISION (A-Z) C C--------------------------------------------------------------------- COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0 A=6378206.4 B=6356584.8 M0=0.9996 D2R=.0174533 E2=(A*A-B*B)/(A*A) EPS=(A*A-B*B)/(B*B) E4=E2*E2 E6=E4*E2 E8=E6*E2 AE=A*(1.0-E2) A0=AE*(1.+(3./4.)*E2+(45./64.)*E4+(175./256.)*E6 1 +(11025./16384.)*E8) A1=-AE/2.*((3./4.)*E2+(15./16)*E4+(525./512.)*E6 1 +(2205./2048.)*E8) A2=AE/4.*((15./64.)*E4+(105./256.)*E6+(2205./4096.)*E8) ERR=-AE/6.*((35./512.)*E6+(315./2048.)*E8) C PRINT *,'A-CONSTANTS ',A0,A1*2,A2*2,ERR*2 PHI0=BASLAT*D2R IZONE=AINT((180.-LAM)/6.+1.) LAM0=(183.-(6.*FLOAT(IZONE)))*D2R C PRINT *,'ZONE NUMBER IS ',IZONE C PRINT *,'CENTRAL MERIDIAN IS ',LAM0/D2R RETURN END C********************************************************************* C CALCULATE THE DISTANCE FROM THE EQUATOR ALONG THE MERIDION FUNCTION SS(PHI) IMPLICIT DOUBLE PRECISION (A-Z) C-------------------------------------------------------------------- COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0 LAT=PHI*D2R SS=(A0*LAT + A1*DSIN(2.D0*LAT) + A2*DSIN(4.D0*LAT) + * ERR*DSIN(6.D0*LAT)) RETURN END C*********************************************************************** C COMPUTE THE RECTANGULAR COORDINATES X AND Y FROM THE C LATITUDE AND LONGITUDE C C1 CHANGED FROM UTMF TO XYTRAN 01/09/85 SUBROUTINE XYTRAN(LAT,LON,X,Y) IMPLICIT DOUBLE PRECISION (A-Z) C C--------------------------------------------------------------------- COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0 D=LAM0-LON*D2R D2=D*D D3=D2*D D4=D3*D D5=D4*D D6=D5*D D7=D6*D D8=D7*D C S=DSIN(LAT*D2R) C=DCOS(LAT*D2R) T=S/C T2=T*T T4=T2*T2 T6=T4*T2 C2=C*C C3=C2*C C4=C3*C C5=C4*C C6=C5*C C7=C6*C C N2=EPS*C*C N4=N2*N2 N=A/SQRT(1.-E2*S*S) AA0=N*C AA1=N*C3*(1.-T2+N2) AA2=N*C5*(5.-18.*T2+T4+14.*N2 1-58.*T2*N2) AA3=N*S*C7*(61.-479.*T2+179.*T4-T6) B0=SS(LAT) B1=N*S*C B2=N*S*C3*(5.-T2+9.*N2+4.*N4) B3=N*S*C5*(61.-58.*T2+T4+270.*N2 1-330.*T2*N2) B4=N*S*C7*(1385.-3111.*T2+543.*T4-T6) C X=M0*(AA0*D+AA1*D3/6.+AA2*D5/120.+AA3*D7/5040.)+500000. Y=M0*(B0+B1*D2/2.+B2*D4/24.+B3*D6/720.+B4*D8/40320) RETURN END C********************************************************************* C FUNCTION TO COMPUTE THE INVERSE OF FUNCTION SS. IT FINDS THE C LATITUDE OF A POINT A DISTANCE Y FROM THE EQUATOR. PHI IS A C INITIAL GUESS AT THE LATITUDE OR A POINT NEAR Y C FUNCTION SINV(Y,PHI) IMPLICIT DOUBLE PRECISION (A-Z) C C--------------------------------------------------------------------- COMMON /DATA1/ D2R,A,B,E2,EPS,A0,A1,A2,LAM0,ERR,M0,PHI0,Y0 Y0=SS(PHI/D2R) DY=Y-Y0 C=A*(1.-E2) D=(1.-E2*DSIN(PHI)*DSIN(PHI)) C1=-3.*E2/(C*C) P1=D**(1.5)/C P2=C1*D*D*DCOS(PHI) SINV=(PHI+DY*P1+DY*DY*P2) RETURN END