C program dcov2.f / revised by P.laidlaw/Jan 2003/ force 2.08 compiler C PROGRAM DCOV.FOR C DEC 1990 / J.MOLTNER C DETERMINE CO-ORDINATES OF CONTOURS C LINEAR INTERPOLATION OF HAATS C SPLINE ERP'S TO DETERMINE PATTERN C DOC METHOD FOR F50,50 CURVES C OUTPUT DATA FILE:DCOV.OUT C STANDARD INPUT FILES: DH5050.DAT (DISTANCES AND HEIGHTS) C : L5050.DAT (LOW BAND 5050 CURVES) C : H5050.DAT (HIGH BAND F5050 CURVES) C : U5050.DAT (UHF F5050 CURVES) C * INPUT FILE FORMAT * City and channel * Band(M,L,H or U)/Power in Kw(use K) or dbK(use D) * Contour to plot( in dBu)/Contour extention in KM * Lat Coordinates like 44 6 15 * Long Coords like 81 36 16 * Num of extra known hats and erps (n) * azimuth1/HAT1/ERP1 * azimuth2/HAT2/ERP3 * azimuthn/HATn/ERPn * azimuth start/azimuth stop/increment * REAL DD(50),HH(21),V1(50),V2(50),EE(50,21) real base_alt,towerheight,Eval REAL AZ1(100),ERP1(100),HAT1(100),ERP2(720),HAT2(720),DIST(720) CHARACTER DATFIL*20,PWRFLG,BND,STN*20,outfile*20,blank*10 character ext*10 WRITE(6,*)'Program DCOV Computes FM/TV coverage from 5050 curves' write(6,*) WRITE(6,200) 200 format($,'Enter the Input DATA FILE name => ' ) READ(5,5) DATFIL write(6,201) 201 format($,'Enter the output MAPINFO filenaame =>') read(5,5) outfile OPEN(1,FILE=DATFIL) OPEN(4,FILE= 'DCOV.OUT') OPEN(12,FILE='DCOV.err') OPEN(13,FILE='DCOV.COR') OPEN(14,FILE='DCOV.UTM') OPEN(15,FILE= outfile) ! file for mapinfo flagerror=0 C devices C 1 a file entered by the user for data input C 2 DH5050 curves C 4 FILE='DCOV.OUT') C 5 Terminal entry C 6 Terminal output C 12 FILE='DCOV.PLT') C 13 FILE='DCOV.COR') C 14 FILE='DCOV.UTM') C 15 FILE='DCOV.LLL') READ(1,5)STN 5 FORMAT(A20) READ(1,1) BND,PWRFLG 1 FORMAT(A1,1X,A1) WRITE(6,*) ' STATION => ',STN C ------------------------------------------------------- READ(1,*)CONTR,CNEX Read(1,*)base_alt,towerheight WRITE(4,*) ' STATION => ',STN WRITE(6,*) ' BAND => ',BND,' CONTOUR => ',CONTR WRITE(4,*) ' BAND => ',BND WRITE(4,*) ' CONTOUR => ',CONTR WRITE(4,*) ' CONTOUR EXTENSION => ',CNEX READ(1,*) ID,IM,IS WRITE(4,*) ' LAT => ',ID,IM,IS TLAT=ID+IM/60.0+IS/3600.0 READ(1,*) ID,IM,IS WRITE(4,*) ' LNG => ',ID,IM,IS WRITE(4,3) 3 FORMAT(/,3X,'AZ',4X,'HAAT',7X,'ERP',6X,'DIST',6X, *'LAT',18X,'LNG',17X,'EASTING',2X,'NORTHING',/) TLNG=ID+IM/60.0+IS/3600.0 IF (BND.EQ.'M') GOTO 21 C F5050 DISTANCES AND HEIGHTS NH=15 ND=40 OPEN(2,FILE='DH5050.DAT') READ(2,*) (DD(J),J=1,ND) READ(2,*) (HH(J),J=1,NH) write(6,*)'band (BND)=',BND C LOWBAND IF (BND.EQ.'L') THEN OPEN(3,FILE='L5050.DAT') IF (CONTR.NE.68.0.AND.CONTR.NE.47.0) THEN WRITE(6,*) ' ***NON-STANDARD CONTOURS***' END IF END IF C HIGHBAND IF (BND.EQ.'H') THEN OPEN(3,FILE='H5050.DAT') IF (CONTR.NE.71.0.AND.CONTR.NE.56.0) THEN WRITE(6,*) ' ***NON-STANDARD CONTOURS***' END IF END IF C UHF IF (BND.EQ.'U') THEN OPEN(3,FILE='U5050.DAT') IF (CONTR.NE.74.0.AND.CONTR.NE.64.0) THEN WRITE(6,*) ' ***NON-STANDARD CONTOURS***' END IF END IF C F5050 FIELDS DO 20 J=1,ND * write(6,*)J READ(3,*) (EE(J,K),K=1,NH) * write(6,*)J 20 CONTINUE C INPUT KNOWN HATS (m) AND ERPS 21 READ(1,*) NAZ EPMIN=1000000.0 EPMAX=0.0 DO 50 J=1,NAZ READ(1,*) AZ1(J),HAT1(J),ERP1(J) IF (PWRFLG.EQ.'D') ERP1(J)=10.0**(ERP1(J)/10.0) IF (ERP1(J).GT.EPMAX) EPMAX=ERP1(J) IF (ERP1(J).LT.EPMIN) EPMIN=ERP1(J) 50 CONTINUE C ADD EXTRA POINT AT END FOR INTERPOLATION/SPLINE NNAZ=NAZ+1 AZ1(NNAZ)=AZ1(1)+360.0 HAT1(NNAZ)=HAT1(1) ERP1(NNAZ)=ERP1(1) * ----------------------------------------------- * used for terminal entry od assesment contour * WRITE(6,*) ' ENTER AZ START/STOP/INC => ' * READ(5,*) AZSRT,AZSTP,AZINC * ----------------------------------------------- read(1,*)AZSRT,AZSTP,AZINC IF (AZSTP.LE.AZSRT) AZSTP=AZSTP+360.0 NAZI=INT(AZSTP-AZSRT)/AZINC+1 C INTERPOLATE HATS AND SPLINE ERPS watch for low HATTS DO 60 J=1,NAZI AZ=(J-1)*AZINC+AZSRT IF (AZ.GE.360.0) AZ=AZ-360.0 HAT2(J)=F3(NNAZ,AZ1,HAT1,AZ) * write(6,*)hat1(j),hat2(j) * IF (HAT2(J).LT.0.0) HAT2(J)=0.0 ERP2(J)=F1(NNAZ,AZ1,ERP1,AZ) IF (ERP2(J).LT.EPMIN) ERP2(J)=EPMIN IF (ERP2(J).GT.EPMAX) ERP2(J)=EPMAX * write(12,*)j,nnaz,az,erp1(j),erp2(j) 60 CONTINUE WRITE(13,*) NAZI WRITE(14,*) NAZI write(6,*)' No. Azimuth HAAT Distance' ******Start main loop here ****************************************** DO 70 L=1,NAZI AZ=(L-1)*AZINC+AZSRT IF (AZ.GE.360.0) AZ=AZ-360.0 Eval=CONTR-10.0*ALOG10(ERP2(L)) ! scale to the input power level IF (BND.EQ.'M') THEN C CALCULATE MMDS CONTOUR DIST DIST(L)=F2(HAT2(L),Eval) ELSE * ******************************************************************** C CALCULATE TV/FM CONT DIST C DETERMINE VECTOR OF FIELDS AT HAT(L) * ******************************************************************** * start here if(hat2(L) .ge. 30) THEN DO 80 J=1,ND ! 40 points for distances DO 81 K=1,NH ! 15 points along the x axix , haat ! th e5050.dat tables contain values of field for 1 KW ! x axis is haat, y axis is distance V1(K)=EE(J,K) 81 CONTINUE 80 V2(J)=F3(NH,HH,V1,HAT2(L)) DIST(L)=F4(ND,V2,DD,Eval) else ! when we are off the curves !step 1 get d1 from haat=towerheight DO 800 J=1,ND DO 810 K=1,NH 810 V1(K)=EE(J,K) 800 V2(J)=F3(NH,HH,V1,towerheight) DIST(L)=F4(ND,V2,DD,Eval) d1=dist(L) ! step 2 get d2 from haat=hat2+towerheight DO 806 J=1,ND DO 816 K=1,NH 816 V1(K)=EE(J,K) 806 V2(J)=F3(NH,HH,V1,(HAT2(L)+towerheight)) DIST(L)=F4(ND,V2,DD,(Eval+6)) d2=dist(l) !step 3 calculate the actual distance and reset Hat2 and eval DIST(L)=2*d2-d1 end if ! for the case when we are off the curves END IF DIST(L)=DIST(L)+CNEX ! add the contour extention write(12,302)l,az,Hat2(l),dist(l) 302 format(i5,3f10.2) C DETERMINE COORDINATES OF POINT CALL COORD(TLAT,TLNG,AZ,DIST(L),CLAT,CLNG) ND1=INT(CLAT) NM1=INT((CLAT-ND1)*60.0) NS1=NINT(((CLAT-ND1)*60.0-NM1)*60.0) C NS1 is never used RS1=((CLAT-ND1)*60.0-NM1)*60.0 ND2=INT(CLNG) NM2=INT((CLNG-ND2)*60.0) NS2=NINT(((CLNG-ND2)*60.0-NM2)*60.0) RS2=((CLNG-ND2)*60.0-NM2)*60.0 CLNG=-CLNG CALL UTM (1,CLAT,CLNG,NZONE,X,Y) C WRITE(6,2) AZ,HAT2(L),ERP2(L),DIST(L),ND1,NM1,RS1,ND2,NM2,RS2,X,Y WRITE(4,2) AZ,HAT2(L),ERP2(L),DIST(L),ND1,NM1,RS1,ND2,NM2,RS2,X,Y WRITE(13,4) ND1,NM1,RS1,ND2,NM2,RS2,NZONE,X,Y WRITE(14,*) NZONE,X,Y XXX=-ND2-NM2/60.0-RS2/3600.0 YYY=ND1+NM1/60.0+RS1/3600.0 WRITE(15,6)XXX,YYY 70 CONTINUE 2 FORMAT(2F8.2,F10.3,F8.2,2I7,F7.1,2I7,F7.1,2F10.2) 4 FORMAT(2I7,F7.1,2I7,F7.1,I7,2F10.3) 6 FORMAT(F12.6,',',F12.6) WRITE(6,*) ' OUTPUT IN "dcov.out" ' WRITE(6,*) ' AZ/DIST IN "dcov.plt" ' WRITE(6,*) ' COORDS IN "dcov.cor" ' WRITE(6,*) ' ZONE,EAST,NORTH (Km) IN "dcov.utm" ' if (flagerror .eq.1)then write(6,*)'Data is off the 5050 curves,estimates apply' end if END C*********************************************************************** C*********************************************************************** C*********** End of main program *************************************** C*********************************************************************** C*********************************************************************** FUNCTION GCDIS(FLAT,FLON,TLAT,TLON) REAL*8 X,PI,Z PI=4.0*DATAN(1.0D0) Z=180.0/PI X=DSIN(FLAT/Z)*DSIN(TLAT/Z)+DCOS(FLAT/Z)*DCOS(TLAT/Z) &*DCOS((TLON-FLON)/Z) GCDIS=SNGL(6370.14*DACOS(X)) RETURN END C**************************************************** FUNCTION GCANG(FLAT,FLON,TLAT,TLON) REAL*8 X,Y,PI,Z,D,AA PI=4.0*DATAN(1.0D0) Z=180.0/PI X=DSIN(FLAT/Z)*DSIN(TLAT/Z)+DCOS(FLAT/Z)*DCOS(TLAT/Z) &*DCOS((TLON-FLON)/Z) D=DACOS(X) Y=DSIN(TLAT/Z)-DSIN(FLAT/Z)*DCOS(D) X=DSIN(D)*DCOS(FLAT/Z) Y=Y/X AA=Z*DACOS(Y) IF (DSIN((TLON-FLON)/Z).GE.0.D0) AA=360.0-AA GCANG=SNGL(AA) RETURN END C******************************* SUBROUTINE COORD(TLAT,TLNG,AZ,D,CLAT,CLNG) C DETERMINE COORD'S OF PT AT D KM AND AZ DEG FROM TLAT,TLNG C ITERATE UNTIL DOC DIST FORMULA IS SATISFIED D1=D JCNT=0 10 JCNT=JCNT+1 C FIND COORDS USING GREAT CIRCLE CALL FINDP(TLAT,TLNG,AZ,D1,CLAT,CLNG) C CHECK DOC DISTANCE D2=DOCDIS(TLAT,TLNG,CLAT,CLNG) C EXIT IF DOC DIST WITHIN 1 METRE OR 100 ITERATIONS DONE IF (ABS(D2-D).LT.0.001) GOTO 20 IF (JCNT.GT.100) THEN WRITE(6,*) AZ,' DEG DIST NOT CONVERGING' WRITE(6,*) ' DIST =',D WRITE(6,*) ' DOC DIST =',D2 WRITE(6,*) ' EXITING ' GOTO 20 END IF D1=D1+(D-D2)/2.0 GOTO 10 20 RETURN END C********************************* FUNCTION DOCDIS(TLAT,TLNG,CLAT,CLNG) A=(TLAT+CLAT)/2.0 A=A/57.29577951 DLAT=111.108-0.566*COS(2.0*A) DLNG=111.391*COS(A)-0.095*COS(3.0*A) B=DLAT*(TLAT-CLAT) C=DLNG*(TLNG-CLNG) DOCDIS=SQRT(B*B+C*C) RETURN END C********************************* FUNCTION F2(H,Eval) C CALCULATE MMDS CONTOUR DISTANCE C FREE SPACE D1=10.0**((74.77-E)/20.0) C GRAZING D2=10.0**((79.0-Eval+10.0*ALOG10(H))/40.0) C LOS HORIZON D3=3.54*SQRT(H) IF (D1.LT.D2) THEN F2=D1 RETURN ELSE IF (D2.LT.D3) THEN F2=D2 RETURN END IF F2=D3 END IF RETURN END C****************************************** FUNCTION F3(N,XX,YY,XB) C GIVEN N DATA PAIRS XX,YY FIND Y-VALUE AT XB C USING LINEAR INTERPOLATION C XX VALUES ASCENDING REAL XX(50),YY(50) NN=N-1 IF (XB.LT.XX(1)) THEN WRITE(6,*) XB,' BELOW RANGE' F3=YY(1) RETURN END IF IF (XB.GT.XX(N)) THEN WRITE(6,*) XB,' ABOVE RANGE' F3=YY(N) RETURN END IF DO 10 J=1,NN IF (XB.GE.XX(J).AND.XB.LE.XX(J+1)) THEN F3=(YY(J+1)-YY(J))/(XX(J+1)-XX(J))*(XB-XX(J))+YY(J) ELSE CONTINUE END IF 10 CONTINUE RETURN END C******************************************* FUNCTION F4(N,XX,YY,XB) C GIVEN N DATA PAIRS XX,YY FIND Y-VALUE AT XB C USING LINEAR INTERPOLATION FOR X C LOG INTERPOLATION FOR Y C XX VALUES DESCENDING REAL XX(50),YY(50) NN=N-1 IF (XB.GT.XX(1)) THEN WRITE(6,*) XB,' ABOVE RANGE' F3=YY(1) C F3 is never used RETURN END IF IF (XB.LT.XX(N)) THEN WRITE(6,*) XB,' BELOW RANGE' F3=YY(N) RETURN END IF DO 10 J=1,NN IF (XB.LE.XX(J).AND.XB.GE.XX(J+1)) THEN R=(XB-XX(J))/(XX(J+1)-XX(J)) F4=YY(J)*10.0**((ALOG10(YY(J+1))-ALOG10(YY(J)))*R) ELSE CONTINUE END IF 10 CONTINUE RETURN END C******************************************* FUNCTION F1(N,XX,YY,XB) C GIVEN N DATA PAIRS XX,YY FIND Y-VALUE AT XB C USING PERIODIC SPLINE REAL XX(50),YY(50),C(4,50) CALL SPLNPER(N,XX,YY,C) XB1=XB IF (XB.LT.XX(1)) XB1=XX(N)-(XX(1)-XB) IF (XB.GT.XX(N)) XB1=XX(1)+(XB-XX(N)) F1=EPCUB(N,XB1,XX,C) RETURN END C********************************* SUBROUTINE SPLNPER(N,X,Y,C) DIMENSION X(50),Y(50),C(4,50),A(50,50),B(50),WORK(50) INTEGER IPVT(50) N1=N-1 DO 30 J=2,N1 DO 5 K=1,N A(J,K)=0.0 5 CONTINUE JP=J+1 JM=J-1 H=X(JP)-X(J) HM=X(J)-X(JM) B(J)=(Y(JP)-Y(J))/H-(Y(J)-Y(JM))/HM A(J,JM)=HM A(J,J)=2.0*(HM+H) A(J,JP)=H 30 CONTINUE DO 40 K=1,N A(1,K)=0.0 A(N,K)=0.0 40 CONTINUE H1=(X(2)-X(1)) HN1=X(N)-X(N-1) A(1,1)=2.0*H1 A(1,2)=H1 A(1,N-1)=HN1 A(1,N)=2.0*HN1 B(1)=(Y(2)-Y(1))/H1-(Y(N)-Y(N-1))/HN1 A(N,1)=1.0 A(N,N)=-1.0 B(N)=0.0 CALL DECOMP(50,N,A,COND,IPVT,WORK) CALL SOLVE(50,N,A,B,IPVT) DO 20 J=1,N-1 JP=J+1 H=X(JP)-X(J) C(1,J)=Y(J) C(2,J)=(Y(JP)-Y(J))/H-H*(B(JP)+2.0*B(J)) C(3,J)=3.0*B(J) C(4,J)=(B(JP)-B(J))/H 20 CONTINUE RETURN END C********************************************** SUBROUTINE DECOMP(NDIM,N,A,COND,IPVT,WORK) INTEGER NDIM,N REAL A(NDIM,N),COND,WORK(N) INTEGER IPVT(N) IPVT(N)=1 IF (N.EQ.1) GOTO 80 NM1=N-1 ANORM=0.0 DO 10 J=1,N T=0.0 DO 5 I=1,N T=T+ABS(A(I,J)) 5 CONTINUE IF (T.GT.ANORM) ANORM=T 10 CONTINUE DO 35 K=1,NM1 KP1=K+1 M=K DO 15 I=KP1,N IF (ABS(A(I,K)).GT.ABS(A(M,K))) M=I 15 CONTINUE IPVT(K)=M IF (M.NE.K) IPVT(N)=-IPVT(N) T=A(M,K) A(M,K)=A(K,K) A(K,K)=T IF (T.EQ.0.0) GOTO 35 DO 20 I=KP1,N A(I,K)=-A(I,K)/T 20 CONTINUE DO 30 J=KP1,N T=A(M,J) A(M,J)=A(K,J) A(K,J)=T IF (T.EQ.0.0) GOTO 30 DO 25 I=KP1,N A(I,J)=A(I,J)+A(I,K)*T 25 CONTINUE 30 CONTINUE 35 CONTINUE DO 50 K=1,N T=0.0 IF (K.EQ.1) GOTO 45 KM1=K-1 DO 40 I=1,KM1 T=T+A(I,K)*WORK(I) 40 CONTINUE 45 EK=1.0 IF (T.LT.0.0) EK=-1.0 IF (A(K,K).EQ.0.0) GOTO 90 WORK(K)=-(EK+T)/A(K,K) 50 CONTINUE DO 60 KB=1,NM1 K=N-KB T=0.0 KP1=K+1 DO 55 I=KP1,N T=T+A(I,K)*WORK(K) 55 CONTINUE WORK(K)=T M=IPVT(K) IF (M.EQ.K) GOTO 60 T=WORK(M) WORK(M)=WORK(K) WORK(K)=T 60 CONTINUE YNORM=0.0 DO 65 I=1,N YNORM=YNORM+ABS(WORK(I)) 65 CONTINUE CALL SOLVE(NDIM,N,A,WORK,IPVT) ZNORM=0.0 DO 70 I=1,N ZNORM=ZNORM+ABS(WORK(I)) 70 CONTINUE COND=ANORM*ZNORM/YNORM IF (COND.LT.1.0) COND=1.0 RETURN 80 COND=1.0 IF (A(1,1).NE.0.0) RETURN 90 COND=1.0E32 RETURN END C****************************************************** SUBROUTINE SOLVE(NDIM,N,A,B,IPVT) REAL A(NDIM,N),B(N) INTEGER IPVT(N) IF (N.EQ.1) GOTO 50 NM1=N-1 DO 20 K=1,NM1 KP1=K+1 M=IPVT(K) T=B(M) B(M)=B(K) B(K)=T DO 10 I=KP1,N B(I)=B(I)+A(I,K)*T 10 CONTINUE 20 CONTINUE DO 40 KB=1,NM1 KM1=N-KB K=KM1+1 B(K)=B(K)/A(K,K) T=-B(K) DO 30 I=1,KM1 B(I)=B(I)+A(I,K)*T 30 CONTINUE 40 CONTINUE 50 B(1)=B(1)/A(1,1) RETURN END C************************************** FUNCTION EPCUB(N,XBAR,X,C) DIMENSION X(N),C(4,N) IF (N .GT. 1) GO TO 5 EPCUB = C(2,1)*(XBAR - X(1)) + C(1,1) RETURN 5 ILOW = 1 IUP = N 10 IMID = (IUP+ILOW)/2 IF (IUP-ILOW .LE. 1) GO TO 30 IF (XBAR .LE. X(IMID)) GO TO 20 ILOW = IMID GO TO 10 20 IUP = IMID GO TO 10 30 DX = XBAR - X(ILOW) EPCUB = C(1,ILOW) + DX*(C(2,ILOW) + DX*(C(3,ILOW) * + DX*C(4,ILOW))) RETURN END C************************************************* SUBROUTINE SOLVT (TH1,TH2,PHI,S,A,B) * To solve a spherical triangle where two sides TH1 and TH2 * and the included angle PHI are known. * On return, S is the third side and A is the angle (azimuth) * between S and TH1. * B is the angle between S and TH2. * The unknown angles are found in terms of tangents * rather than sines or cosines in order to avoid the indeterminancies * that occur close to +- 1. * Also the angles are placed in the correct quadrant, -pi to pi. * All quantities in this subroutine are in units of radians. * Reference: Astronomie General, by Andre Danjon, pp 16-18. COSTH1 = COS(TH1) COSTH2 = COS(TH2) SINTH1 = SIN(TH1) SINTH2 = SIN(TH2) COSPH = COS(PHI) SINPH = SIN(PHI) * Do not allow SINPH to be zero, since this leads to 0/0 later. IF (ABS(SINPH) .LT. 1.E-10) SINPH = 1.E-10 X = SINTH1*COSTH2/SINTH2 - COSTH1*COSPH A = ATAN2 (SINPH,X) SINA = SINPH / SQRT (SINPH*SINPH + X*X) * (sin(A) is not used because its relative precision suffers near pi) X = SINTH2*COSTH1/SINTH1 - COSTH2*COSPH B = ATAN2 (SINPH,X) SINS = SINPH * SINTH2/SINA COSS = COSTH1*COSTH2 + SINTH1*SINTH2*COSPH S = ATAN2 (SINS,COSS) RETURN END * ------------------------------------------------------------------- SUBROUTINE FINDP (GLAT1,GLONG1,A,S,GLAT,GLONG) * Define a great circle path passing through point GLAT1,GLONG1 * with azimuthal angle A. * This subroutine finds the point along the path that is a given * distance S from the starting point. * Latitudes and longitudes are in degrees, but * A is in degrees, and S is in km. PARAMETER (DR = .0174532925199433, RD = 57.2957795130823) THETA1 = DR * (90.-GLAT1) * We know 2 sides (S and THETA1), and the included angle (A). * Therefore use subroutine SOLVT. aa=a*dr ss=s/6370.14 CALL SOLVT (SS,THETA1,AA,THETA2,B,PHI) GLAT = 90. - RD*THETA2 GLONG = GLONG1 - RD*PHI RETURN END C******************************************************************* SUBROUTINE UTM (INSTR,GLAT,GLONG,NZONE,X,Y) DATA PI,DR,RD/3.14159265358979,.0174532925199433,57.2957795130823/ DATA A,REQ,ESQ / 6378.2064, 6335.0345, .006768658/ DATA E0,E2,E4 / 1.00510892, .00255988, .00000272/ DATA E20,E40,E24 / .00254687, .00000270, .00000378/ IF (INSTR .GT. 1) GO TO 10 PHI = GLAT*DR SINPH = SIN(PHI) COSPH = COS(PHI) ESINSQ = ESQ * SINPH * SINPH RT = 1. - .5*ESINSQ - .125*ESINSQ*ESINSQ RNU = A / RT RHO = REQ / (RT*RT*RT) SIN2PH = 2. * SINPH * COSPH COS2PH = COSPH*COSPH - SINPH*SINPH SIN4PH = 2. * SIN2PH * COS2PH YM = REQ * (E0*PHI - E2*SIN2PH + E4*SIN4PH) ZONE = AMOD ((GLONG+180.)/6., 60.) IF (ZONE .LT. 0.) ZONE = ZONE + 60. NZONE = ZONE + 1. C LONGITUDE DIFFERENCE FROM CENTRAL MERIDIAN MERID = 6*NZONE - 183 DIFF = GLONG - FLOAT(MERID) NCIRC = (DIFF + SIGN(180.,DIFF)) / 360. MERID = MERID + 360*NCIRC DIFF = GLONG - FLOAT(MERID) W = DIFF * DR W2 = W*W W3 = W*W2 C SOME FUNCTIONS OF LATITUDE COS2 = COSPH*COSPH COS3 = COSPH*COS2 T2 = SINPH*SINPH / COS2 C C NORTHING: Y = YM + W2 *.5*RNU * SINPH*COSPH * (1. + (W2/12.)*COS2*(5.-T2)) C (THE LAST TERM W2*W2/12 ETC. IS AN APPROX TO THE VERY SMALL W**4 TERM) C C EASTING: X = W*RNU*COSPH + W3*(RNU/6.)*COS3*(RNU/RHO-T2) C SCALE FACTOR X = X * 0.9996 + 500. Y = Y * 0.9996 C FALSE NORTHING FOR SOUTHERN HEMISPHERE IF (Y .LT. 0.) Y = -(Y + 10000.) C RETURN C C--------------------------------------------------------------------- C C C TRANSFORM UTM COORDINATES TO LATITUDE AND LONGITUDE C C FALSE NORTHING FOR SOUTHERN HEMISPHERE 10 YY = Y IF (Y .LT. 0.) YY = -Y - 10000. C C SCALE FACTOR XX = (X - 500.) / .9996 YY = YY / .9996 C C FIND CENTRAL MERIDIAN ZLONG = 6.*FLOAT(NZONE) + 177. IF (ZLONG .GT. 180.) ZLONG = ZLONG - 360. C C FIND LATITUDE CORESPONDING TO Y AS A MERIDIAN LENGTH C ZERO ORDER APPROXIMATION: PHI0 = YY / (REQ*E0) C FINAL VALUE: PHI1 = PHI0 + E20*SIN(2.*PHI0) + E24*SIN(4.*PHI0) C SINPH = SIN(PHI1) COSPH = COS(PHI1) ESINSQ = ESQ * SINPH * SINPH RT = 1. - .5*ESINSQ - .125*ESINSQ*ESINSQ RNU = A / RT RNU3 = RNU*RNU*RNU RHO = REQ / (RT*RT*RT) T = SINPH / COSPH T2 = T*T X2 = XX*XX X3 = XX*X2 X4 = X2*X2 C LATITUDE CALCULATION C PHI = PHI1 - X2*T/(2.*RHO*RNU) + (X4*T/(24.*RNU3*RHO))*(5.+ 3.*T2) C C LONGITUDE CALCULATION W = XX/(COSPH*RNU) - (X3/(6.*COSPH*RNU3))*(RNU/RHO + 2.*T2) C C CHANGE TO DEGREES GLAT = PHI * RD GLONG = ZLONG + W*RD RETURN END