/* ll2utm: convert lat long to UTM coordinates. *( * -x long or X longitude or UTM X [default is stdin] * -y lat or Y latitude or UTM Y [default is stdout] * * -z zone UTM zone [default is 11] * -s spheroid (1-19) [default is 1] * -ll2utm convert lat/long to UTM [default is yes] * -utm2ll convert UTM back to lat/long * * * 15Dec94 David Okaya Initial installation. */ #include #include #include #include #include #define TRUE 1 #define FALSE 0 FILE *fpin, *fpout; char *malloc(); static unsigned char *buf, *work, *bh, *trchd; char infile[80]; char outfile[80]; int zone; /* UTM zone */ int ispher; /* Spheroid for projection (1-19) */ int ll2utm_flag; /* 1=TRUE=LL->UTM, 0=FALSE=UTM->LL */ float xin, yin; /* lat/long or UTM x,y in */ main(c,v) int c; char *v[]; { int i; float xout,yout; char vv[50][80]; set_defaults(); for(i=0; i=2) /* decode possible user choices */ {i=0; /* i=command line choices; 0='usc_segylmo'*/ while(i=c){errflag=10;} else {ispher= atoi(vv[i]); if(ispher<0 || ispher>20)errflag=10; } break; case 2: /* -z UM zone */ ++i; if(i>=c){errflag=10;} else {zone=atoi(vv[i]); if(zone< -999999 || zone>999999)errflag=10; } break; case 3: /* -x */ ++i; if(i>=c){errflag=10;} else {xwork=atof(vv[i]); xin=xwork; if(xin< -99999999. || xin>99999999.)errflag=10; } break; case 4: /* -y */ ++i; if(i>=c){errflag=10;} else {xwork=atof(vv[i]); yin=xwork; if(yin< -99999999. || yin>99999999.)errflag=10; } break; case 5: /* -ll2utm */ ll2utm_flag=TRUE; break; case 6: /* -utm2ll */ ll2utm_flag=FALSE; break; default: /* can't interpret (userchoice=0) */ errflag=10; break; } if(errflag==10) {switch(userchoice) {case 0: fprintf(stderr,"Error: Invalid command line option\n"); break; case 1: fprintf(stderr,"Error: -s value out of range\n"); break; case 2: fprintf(stderr,"Error: -z value invalid\n"); break; case 3: fprintf(stderr,"Error: -x value invalid \n"); break; case 4: fprintf(stderr,"Error: -y value invalid \n"); break; } } } /* end of while */ } if(xin==0 && yin==0) {errflag=10;} if(errflag>0 ) {printf("\nll2utm -x -y [-ll2utm -utm2ll -z -s] \n"); printf(" -x long_X : longitude or UTM_X\n"); printf(" -y lat_Y : latitude or UTM_Y\n"); printf(" -ll2utm : convert lat/long to UTM [default is yes]\n"); printf(" -utm2ll : convert UTM to lat/long [not yet implemented]\n"); printf(" -z zone : UTM zone [default is 11]\n"); printf(" -s sphere : spheroid for conversion [default is 1]\n"); printf(" 1=CLARKE 1866; 2=CLARKE 1880;\n"); printf(" 3=BESSEL 1967; 4=NEW INTERNATIONAL 1967;\n"); printf(" 5=INTERNATIONAL 1909 (1924) (HAYFORD);\n"); printf(" 6=WORLD GEODETIC SYSTEM 1972; 7=EVEREST;\n"); printf(" 8=WORLD GEODETIC SYSTEM 1966;\n"); printf(" 9=GEODETIC REFERENCE SYSTEM 1980;\n"); printf(" 10=AIRY; 11=MODIFIED EVEREST; 12=MODIFIED AIRY;\n"); printf(" 13=WALBECK; 14=SOUTH ASIA 1960 (FISCHER);\n"); printf(" 15=AUSTRALIAN NATIONAL-SOUTH AMERICAN 1969;\n"); printf(" 16=KRASSOVSKY; 17=HOUGH; 18=MERCURY 1960 (FISCHER);\n"); printf(" 19=MODIFIED MERCURY 1968 (FISCHER)\n"); printf(" 20=WORLD GEODETIC SYSTEM 1984\n"); exit(1);} } /************************************************************************/ /* set_defaults: set default behavior * #define TRUE 1 #define FALSE 0 */ set_defaults() { zone=11; /* UTM zone: California is zone 11 */ ispher=20; /* Spheroid for projection */ ll2utm_flag=TRUE; /* 1=TRUE=LL->UTM, 0=FALSE=UTM->LL */ } /* ll2utm subroutine: shell in front of conversion * subroutines which were obtained from USGS software * originally provided in FORTRAN. * * direction flag: =1 ll2utm, =0 utm2ll * * * 15Dec94 David Okaya Initial conversion. Direction=0 NOT YET INSTALLED. * 19Dec94 David Okaya Install WGS84 spheroid as spheroid=20. * axis = 6378137. * 1/f = 298.257222101 * where f=flattening = (a-b)/a * so b_axis = a(1-f) * = 6356752.31414 */ #include sub_ll2utm(xin,yin,direction_flag,zone,spher,xout,yout) float xin,yin ; /* Should be x.xxxxx */ int direction_flag, zone, spher ; float *xout, *yout; { double x,y,flat,flon,conv; int izone,isph,idir; izone = zone; isph = spher; if (direction_flag == 1) {idir=0; flat=yin; flon=xin; conv =0; plane_for(&x, &y,flat,flon,conv,izone,idir,isph); *xout = (float) x; *yout = (float) y; } else if(direction_flag == 0) {idir=1; x=yin; y=xin; conv =0; plane_inv(x, y,&flat,&flon,conv,izone,idir,isph); *xout = (float) flat; *yout = (float) flon; } } /* SUBROUTINE PLANE (X,Y,FLAT,FLON,CONV,KKZONE,IDIR,ISPHER) C C GENERALIZED PLANE COORDINATE CONVERSION PACKAGE 7-30-82 FIPS CODE C MODIFIED 6-16-83 IMPLICIT REAL*8 (A-H,O-Z) REAL*8 X,Y,FLAT,FLON,CONV C C C X,Y -- UTM IN METERS C FLAT,FLON -- LATITUDE AND LONGITUDE IN DECIDEGREES C CONV -- CONVERGENCE IN DECIDEGREES C IZONE -- ZONE NUMBER C IDIR -- DIRECTION OF CONVERSION: 1 = INVERSE; 0 = FORWARD C FORWARD IS GEOGRAPHIC TO PLANE C INVERSE IS PLANE COORDINATES TO GEOGRAPHIC C ISPHER -- SPHEROID SELECTION: 1=CLARKE 1866; 2=CLARKE 1880; C 3=BESSEL 1967; 4=NEW INTERNATIONAL 1967; C 5=INTERNATIONAL 1909 (1924) (HAYFORD); C 6=WORLD GEODETIC SYSTEM 1972; 7=EVEREST; C 8=WORLD GEODETIC SYSTEM 1966; C 9=GEODETIC REFERENCE SYSTEM 1980; C 10=AIRY; 11=MODIFIED EVEREST; 12=MODIFIED AIRY; C 13=WALBECK; 14=SOUTH ASIA 1960 (FISCHER); C 15=AUSTRALIAN NATIONAL-SOUTH AMERICAN 1969; C 16=KRASSOVSKY; 17=HOUGH; 18=MERCURY 1960 (FISCHER); C 19=MODIFIED MERCURY 1968 (FISCHER) C NOTES: C IF ISPHER.EQ.0: COMPUTE IN LAMBERT OR TRANSVERSE MERCATOR C IF ISPHER.NE.0: COMPUTE IN UNIVERSAL TRANSVERSE MERCATOR C IZONE IS UTM ZONE NUMBER OR STATE PLANE ZONE DEPENDING ON C SELECTION OF ISPHER. IF UTM AND SOUTHERN HEMISPHERE, IZONE C IS NEGATIVE. IF UTM AND IDIR=0, IZONE MAY BE ZERO C c izone will = 11 for California IZONE=KKZONE SECS = 3600.0D0 C C ALL SUBROUTINE INPUT AND OUTPUT OF ANGULAR DATA IS IN C SECONDS. SO CONVERT DECIDEGREES TO SECONDS. 6001 IF (IDIR.NE.0) GO TO 10 ALAT = FLAT * SECS ALON = FLON * SECS ALON = - ALON 10 XX = X YY = Y C CONVERSION STATEMENT ZONES SPHERE C UTM 400 1-60 NON-ZERO 400 CALL UTM(XX,YY,ALAT,ALON,CONV,IZONE,IDIR,ISPHER) C C AFTER THE CONVERSION -- CHANGE BACK TO DECIDEGREES 500 FLAT = ALAT / SECS IF (KKZONE.EQ.0) KKZONE=IZONE FLON = -ALON / SECS CONV = CONV / SECS IF (IDIR.NE.0) RETURN 600 X = XX Y = YY RETURN 900 WRITE(6,901) IZONE 901 FORMAT(1X,I6,' ZONE NUMBER IS INVALID') RETURN END */ /* forward: from ll to UTM */ plane_for (X,Y,FLAT,FLON,CONV,KKZONE,IDIR,ISPHER) double *X, *Y,FLAT,FLON,CONV; int KKZONE,IDIR,ISPHER; { int IZONE; double SECS,ALAT,ALON,XX,YY; SECS = 3600.0; IZONE = KKZONE; /* ALL SUBROUTINE INPUT AND OUTPUT OF ANGULAR DATA IS IN SECONDS. SO CONVERT DECIDEGREES TO SECONDS. */ if(IDIR==0) /* this is true */ { ALAT = FLAT * SECS; ALON = FLON * SECS; ALON = - ALON; } *X = 0; *Y = 0; XX = *X; YY = *Y; TO_UTM(&XX,&YY,ALAT,ALON,CONV,IZONE,IDIR,ISPHER); /* AFTER THE CONVERSION -- CHANGE BACK TO DECIDEGREES */ /* FLAT = ALAT / SECS if (KKZONE == 0) KKZONE=IZONE; FLON = -ALON / SECS; CONV = CONV / SECS; */ if(IDIR == 0) /* this is true */ { *X = XX; *Y = YY; } } /* reverse: from UTM to ll */ plane_inv (X,Y,FLAT,FLON,CONV,KKZONE,IDIR,ISPHER) double X,Y, *FLAT,*FLON,*CONV; int KKZONE,IDIR,ISPHER; { int IZONE; double SECS,ALAT,ALON,XX,YY,WCONV; SECS = 3600.0; IZONE = KKZONE; /* ALL SUBROUTINE INPUT AND OUTPUT OF ANGULAR DATA IS IN SECONDS. SO CONVERT DECIDEGREES TO SECONDS. */ if(IDIR==0) /* this is FALSE */ { ALAT = *FLAT * SECS; ALON = *FLON * SECS; ALON = - ALON; } XX = X; YY = Y; WCONV = *CONV; FROM_UTM(XX,YY,&ALAT,&ALON,&WCONV,IZONE,IDIR,ISPHER); /* AFTER THE CONVERSION -- CHANGE BACK TO DECIDEGREES */ *FLAT = ALAT / SECS; if (KKZONE == 0) KKZONE=IZONE; *FLON = -ALON / SECS; *CONV = WCONV / SECS; if(IDIR == 0) /* this is FALSE */ { X = XX; Y = YY; } } /* SUBROUTINE UTM(X,Y,SLAT,SLON,CONV,IIZONE,IDIR,ISPHER) C C UNIVERSAL TRANSVERSE MERCATOR CONVERSION C C REWRITTEN 6-16-83 BY J.F. WAANANEN USING FORMULATION C BY J.P. SNYDER IN BULLETIN 1532, PAGES 68-69 C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 AXIS(19),BXIS(19) C DATA AXIS/6378206.4D0,6378249.145D0,6377397.155D0, . 6378157.5D0,6378388.D0,6378135.D0,6377276.3452D0, . 6378145.D0,6378137.D0,6377563.396D0,6377304.063D0, . 6377341.89D0,6376896.0D0,6378155.0D0,6378160.D0, . 6378245.D0,6378270.D0,6378166.D0,6378150.D0/ C DATA BXIS/6356583.8D0,6356514.86955D0,6356078.96284D0, . 6356772.2D0,6356911.94613D0,6356750.519915D0,6356075.4133D0, . 6356759.769356D0,6356752.31414D0,6356256.91D0,6356103.039D0, . 6356036.143D0,6355834.8467D0,6356773.3205D0,6356774.719D0, . 6356863.0188D0,6356794.343479D0,6356784.283666D0, . 6356768.337303D0/ C DATA AK0/0.9996D0/ C RADSEC = 206264.8062470964D0 C A = AXIS(ISPHER) B = BXIS(ISPHER) ES= (A**2-B**2)/A**2 C IZONE = IIZONE IF (IZONE.NE.0) GO TO 10 C C COMPUTE UTM ZONE(IZONE) AND CENTRAL MERIDIAN IN SECONDS FOR C GEODETIC TO UTM CONVERSION WHERE ZONE IS NOT INPUT. C IF (IDIR.NE.0) GO TO 150 IIZ=SLON/2.16D4 IZONE = 30 - IIZ IF (SLON.LT.0.0D0) IZONE=31- IIZ 10 IF(IABS(IZONE).GT.30) GO TO 20 IUTZ = 30-IABS(IZONE) CM=((IUTZ*6.0D0)+3.0D0)*3600.0D0 GO TO 30 20 IUTZ = IABS(IZONE)-30 CM=((IUTZ*6.0D0)-3.0D0)*(-3600.0D0) C 30 IF (IDIR.EQ.0.AND.SLAT.LT.0.0D0) IZONE=-IABS(IZONE) IF (IDIR.EQ.0.AND.IZONE.LT.0) SLAT = -DABS(SLAT) IF (IIZONE.EQ.0) IIZONE=IZONE IF (SLAT.LT.0.0D0.AND.IIZONE.GT.0) IIZONE=-IIZONE IF (IDIR.NE.0) GO TO 100 C-------------------------------------------------------------- C FORWARD COMPUTATION C-------------------------------------------------------------- C PHI = SLAT/RADSEC DLAM = -(SLON-CM)/RADSEC C EPRI = ES/(1.0D0-ES); EN = A/ sqrt(1.00-ES*sin(PHI)*sin(PHI)); T = TAN(PHI)*tan(PHI); C = EPRI*cos(PHI)*cos(PHI); AA = DLAM*cos(PHI); S2 = sin(2.0*PHI); S4 = sin(4.0*PHI); S6 = sin(6.0*PHI); F1 = 1. - (ES/4.) - (3.0*ES*ES/6.41)-(5.D0*ES*ES*ES/2.56D2)) F2 = ((3.D0*ES/8.D0)+(3.D0*ES*ES/3.2D1)+(4.5D1*ES*ES*ES/1.024D3)) F3 = ((1.5D1*ES*ES/2.56D2)+(4.5D1*ES*ES*ES/1.024D3)) F4 = (3.5D1*ES*ES*ES/3.072D3) EM = A*(F1*PHI - F2*S2 + F3*S4 - F4*S6) C XX = AK0*EN*(AA + (1.D0-T+C)*AA**3/6.D0 + (5.D0-1.8D1*T+T*T+ . 7.2D1*C-5.8D1*EPRI)*AA**5/1.2D2) XX = XX + 5.0D5 C YY = AK0*(EM + EN*DTAN(PHI)*((AA*AA/2.D0) + . (5.D0-T+9.D0*C+4.D0*C*C)*AA**4/2.4D1 + . (6.1D1-5.8D1*T+T*T+6.D2*C-3.3D2*EPRI)*AA**6/7.2D2)) IF (IZONE.LT.0 .OR. SLAT.LT.0.0D0) YY = YY + 1.0D7 C GO TO 200 C-------------------------------------------------------------- C INVERSE COMPUTATION C-------------------------------------------------------------- 100 YY = Y IF (IZONE.LT.0) YY = YY - 1.0D7 XX = X - 5.0D5 EM = YY/AK0 UM = EM/(A*(1.D0-(ES/4.D0)-(3.D0*ES*ES/6.4D1) . -(5.D0*ES*ES*ES/2.56D2))) E1 = (1.D0-DSQRT(1.D0-ES))/(1.D0+DSQRT(1.D0-ES)) C PHI1 = UM+((3.D0*E1/2.D0)-(2.7D1*E1**3/3.2D1))*DSIN(2.D0*UM) + . ((2.1D1*E1*E1/1.6D1)-(5.5D1*E1**4/3.2D1))*DSIN(4.D0*UM) + . (1.51D2*E1**3/9.6D1)*DSIN(6.D0*UM) C EN = A/DSQRT(1.0D0-ES*DSIN(PHI1)**2) T = DTAN(PHI1)**2 EPRI = ES/(1.D0-ES) C = EPRI*DCOS(PHI1)**2 R = (A*(1.D0-ES))/((1.D0-ES*DSIN(PHI1)**2)**1.5D0) D = XX/(EN*AK0) C PHI = PHI1 - (EN*DTAN(PHI1)/R) * ((D*D/2.D0) .-(5.D0+3.D0*T+10.D0*C-4.D0*C*C-9.D0*EPRI)*D**4/2.4D1+(6.1D1+ . 9.D1*T+2.98D2*C+4.5D1*T*T-2.52D2*EPRI-3.D0*C*C)*D**6/7.2D2) C ALAM = (CM/RADSEC)-(D-(1.D0+2.D0*T+C)*D**3/6.D0 + . (5.D0-2.D0*C+2.8D1*T-3.D0*C*C+8.D0*EPRI+2.4D1*T*T)* . D**5/1.2D2)/DCOS(PHI1) C SLAT = PHI*RADSEC SLON = ALAM*RADSEC DLAM = -(SLON-CM)/RADSEC C-------------------------------------------------------------- 200 CONTINUE C C TEST TO SEE IF WITHIN DEFINITION OF UTM PROJECTION C C ---- LATITUDE WITHIN 84 DEGREES IF (DABS(SLAT).GT.302400) GO TO 150 C ---- DELTA LONGITUDE WITHIN 0.16 RADIANS OF CENTRAL MERIDIAN IF (DABS(DLAM).GT.1.6D-1) GO TO 150 C IF (IDIR.NE.0) GO TO 210 X = XX Y = YY C C COMPUTE CONVERGENCE C 210 CONV = DLAM*(DSIN(PHI)+1.9587D-12*(DLAM**2) . *DSIN(PHI)*DCOS(PHI*PHI)) CONV = CONV*RADSEC C 300 RETURN C C ERROR MESSAGE C 150 WRITE(6,1) 1 FORMAT(' ************** ERROR UTM *************') GO TO 300 C END */ TO_UTM(X,Y,SLAT,SLON,CONV,IIZONE,IDIR,ISPHER) double *X, *Y, SLAT, SLON, CONV; int IIZONE,IDIR,ISPHER; { static double AXIS[]={6378206.4, 6378249.145, 6377397.155, 6378157.5, 6378388., 6378135., 6377276.3452, 6378145., 6378137., 6377563.396, 6377304.063, 6377341.89, 6376896.0, 6378155.0, 6378160., 6378245., 6378270., 6378166., 6378150., 6378137.}; static double BXIS[]={6356583.8, 6356514.86955, 6356078.96284, 6356772.2, 6356911.94613, 6356750.519915, 6356075.4133, 6356759.769356, 6356752.31414, 6356256.91, 6356103.039, 6356036.143, 6355834.8467, 6356773.3205, 6356774.719, 6356863.0188, 6356794.343479, 6356784.283666, 6356768.337303, 6356752.31414}; int IZONE,IIZ,IUTZ; double AK0,RADSEC,A,B,ES,CM,XX,YY; double PHI,DLAM,EPRI,EN,T,C,AA,S2,S4,S6,F1,F2,F3,F4,EM; AK0=0.9996; RADSEC = 206264.8062470964; A = AXIS[ISPHER-1]; /* -1 for C, not FORTRAN indexing */ B = BXIS[ISPHER-1]; /* -1 for C, not FORTRAN indexing */ ES= (A*A-B*B)/(A*A); IZONE = IIZONE; /* COMPUTE UTM ZONE(IZONE) AND CENTRAL MERIDIAN IN SECONDS FOR */ /* GEODETIC TO UTM CONVERSION WHERE ZONE IS NOT INPUT. */ if (IZONE == 0) { if (IDIR !=0) {printf("ERROR IN UTM_FORWARD\n"); exit(-1);} IIZ= (int) (SLON/(double)21600.) ; IZONE = 30 - IIZ; if (SLON < (double)0.0) IZONE=31- IIZ; } if(abs(IZONE) <= 30) { IUTZ = 30-abs(IZONE); CM= (double) (((IUTZ*6.00)+3.00)*3600.00); } else if(abs(IZONE) > 30) { IUTZ = abs(IZONE)-30; CM=(double)((IUTZ*6.0)-3.0)*(-3600.0); } /* for forward, IDIR=0; for inverse, IDIR=1 */ if (IDIR == 0 && SLAT < 0.) IZONE = -abs(IZONE); if (IDIR == 0 && IZONE < 0) SLAT = -fabs(SLAT); if (IIZONE==0) IIZONE= IZONE; if (SLAT <0.0 && IIZONE > 0) IIZONE= -IIZONE; /*IF (IDIR.NE.0) GO TO 100*/ /* this is explicitly the forward*/ /*-------------------------------------------------------------- C FORWARD COMPUTATION C-------------------------------------------------------------- */ PHI = SLAT/RADSEC; DLAM = -(SLON-CM)/RADSEC; EPRI = ES/(1.0-ES); EN = A/sqrt(1.0-ES*sin(PHI)*sin(PHI)); T = tan(PHI)*tan(PHI); C = EPRI*cos(PHI)*cos(PHI); AA = DLAM*cos(PHI); S2 = sin(2.0*PHI); S4 = sin(4.0*PHI); S6 = sin(6.0*PHI); F1 = 1.0 - (ES/4.0) - (3.0*ES*ES/64.0) - (5.0*ES*ES*ES/256.0); F2 = (3.0*ES/8.0)+(3.0*ES*ES/32.)+(45.*ES*ES*ES/1024.); F3 = (15.*ES*ES/256.)+(45.*ES*ES*ES/1024.); F4 = 35.*ES*ES*ES/3072.; EM = A*(F1*PHI - F2*S2 + F3*S4 - F4*S6); XX = AK0*EN*(AA + (1.0-T+C)*AA*AA*AA/6.0 + (5.0 -18.*T +T*T +72.*C -58.*EPRI)*AA*AA*AA*AA*AA/120.); XX = XX +500000.; YY = EM + EN*tan(PHI)* ((AA*AA/2.0) + (5.0-T+9.0*C+4.0*C*C)*AA*AA*AA*AA/24. + (61.-58.*T+T*T+600.*C-3300.*EPRI)*AA*AA*AA*AA*AA*AA/720.); YY = AK0 * YY; if (IZONE < 0 || SLAT < 0.) YY = YY + 10000000.; /* GO TO 200 C-------------------------------------------------------------- C INVERSE COMPUTATION C-------------------------------------------------------------- 100 YY = Y IF (IZONE.LT.0) YY = YY - 1.0D7 XX = X - 5.0D5 EM = YY/AK0 UM = EM/(A*(1.D0-(ES/4.D0)-(3.D0*ES*ES/6.4D1) . -(5.D0*ES*ES*ES/2.56D2))) E1 = (1.D0-DSQRT(1.D0-ES))/(1.D0+DSQRT(1.D0-ES)) C PHI1 = UM+((3.D0*E1/2.D0)-(2.7D1*E1**3/3.2D1))*DSIN(2.D0*UM) + . ((2.1D1*E1*E1/1.6D1)-(5.5D1*E1**4/3.2D1))*DSIN(4.D0*UM) + . (1.51D2*E1**3/9.6D1)*DSIN(6.D0*UM) C EN = A/DSQRT(1.0D0-ES*DSIN(PHI1)**2) T = DTAN(PHI1)**2 EPRI = ES/(1.D0-ES) C = EPRI*DCOS(PHI1)**2 R = (A*(1.D0-ES))/((1.D0-ES*DSIN(PHI1)**2)**1.5D0); D = XX/(EN*AK0); C PHI = PHI1 - (EN*DTAN(PHI1)/R) * ((D*D/2.D0) -(5.D0+3.D0*T+10.D0*C-4.D0*C*C-9.D0*EPRI)*D**4/2.4D1+(6.1D1+ 9.D1*T+2.98D2*C+4.5D1*T*T-2.52D2*EPRI-3.D0*C*C)*D**6/7.2D2); C ALAM = (CM/RADSEC)-(D-(1.D0+2.D0*T+C)*D**3/6.D0 + (5.D0-2.D0*C+2.8D1*T-3.D0*C*C+8.D0*EPRI+2.4D1*T*T)* D**5/1.2D2)/DCOS(PHI1); C SLAT = PHI*RADSEC; SLON = ALAM*RADSEC; DLAM = -(SLON-CM)/RADSEC; C-------------------------------------------------------------- 200 CONTINUE */ /* TEST TO SEE IF WITHIN DEFINITION OF UTM PROJECTION /* LATITUDE WITHIN 84 DEGREES */ if (fabs(SLAT) > 302400.) {printf("ERROR IN UTM_FORWARD\n"); exit(-1);} /*DELTA LONGITUDE WITHIN 0.16 RADIANS OF CENTRAL MERIDIAN */ if (fabs(DLAM) > 0.16) {printf("ERROR IN UTM_FORWARD\n"); exit(-1);} if (IDIR == 0.) { *X = XX; *Y = YY; } /* COMPUTE CONVERGENCE */ CONV = DLAM*(sin(PHI)+DLAM*DLAM* 0.0000000000019587 *sin(PHI)*cos(PHI*PHI)); CONV = CONV*RADSEC; } FROM_UTM(X,Y,SLAT,SLON,CONV,IIZONE,IDIR,ISPHER) double X, Y, *SLAT, *SLON, *CONV; int IIZONE,IDIR,ISPHER; { *SLAT=0.0; *SLON=0.0; *CONV=0; }