"	SUBROUTINE MAPROJ (N,LON,LAT,X,Y)JC_Title MAPROJ  Conversion from geographic lat-long to mapping projection.C_ARGUMENTSHC Arguments take different meanings with the change in the control flag.IC When the mode is greater than zero, this is the most commonly used modeIC for this routine.  The arguments below are best described at this mode.  *	INTEGER*4 N	   ![I] Process control flag.H 	REAL*4 LAT(*)	   ![I] Latitude to be converted to mapping coordinates.7	REAL*4 LON(*)	   ![I] Longitude array to be converted.*	REAL*4 X(*)	   ![O] Mapping X coordinate.*	REAL*4 Y(*)	   ![O] Mapping Y coordinate.GC_DESCRIPTION  The subroutine MAPROJ converts from geographic latitude-LC longitude to any one of 6 mapping projections.  For a complete descriptionHC of the algorithms used here refer to: Map Projections Used by the U.S.FC Geological Survey, Geological Survey Bulletin 1532, Snyder, John, P.JC Call once to set up, then repeatedly with arrays of positions. Coded to ?C allow relatively easy change of input and output array types.@C Arguments are dependent upon the sign of the first argument N.CC If N is:CC	-  Initialize a projection.#CC	0  Report on current projection.9CC	+  Convert a set of input lon-lat to cartographic X-Y.GCC IF N is less than 0 then initialization mode, the following applies.4CC     Input is negative of the projection control. <CC	- 1 then Mercator.			0,0 is the origin of the projection./CC	- 2 then Lambert conformal(2 parallels)	pole%CC	- 3 then polar-sterographic.		pole$CC	- 4 then simple cylindrical.		0,0'CC	- 5 then sinusoidal equal-area.		0,0=CC	- 6 then point perspective ="view".	sub-space craft point.,CC	- 7 then no projection (x=lon,y=lat).	0,0CC LON ignored.CC LAT ignored.%CC X(8)  input projection parameters.8CC  X(1) = Body major(equitorial) radius. Defaults to 1.CC  X(2) = Body eccentricity.UCC  X(3) = Input angular units of lon and lat, in radians.  Degress =.0174533 defaultPCC  X(4) = Output units, consistant with body radius and scale. inch=2.54e-5 km.3CC  X(5) = Inverse of projection scale. e.g. 25.e6.?CC	for point-perspective, this is scale in output units/radian.PCC  X(6) to X(8): Meaning varies with projection: ignored unless specified here.CC  for projections 1,2,4,5.?CC  X(6) = Central longitude in degrees for Lambert projection.HCC  X(7) = First standard parallel in degrees, sign sets the hemisphere.JCC  X(8) = Second standard parallel in degrees, the sign is ignored for a "CC  polar sterographic projection:PCC  X(7) = Latitude for defined scale degrees, sign sets hemisphere, default=90.%CC  for point-perspective projection:LCC  X(6:7:8) = (x:y:z) coordinates of observer, in same units as body radiusCC  Y ignored.0CC N=0; Projection parameters are reported back:+CC N	 I   [I] 0, Indicate report requested.PCC LON(1) R*4 [O] Line type code for meridians:  1= straight, 2= circle, 3=otherPCC LAT(1) R*4 [O] Line type code for parallels:  1= straight, 2= circle, 3=otherQCC X(8)   R*4 [O] Reports back the projection parameters, in same order as input.1CC Y(1)   R*4 [O] 4 character name of projection.CC N>0; conversion: ,CC  N	  I   [I] Number of points to convert.1CC  LON(*) R*4 [I] West longitude in input units.,CC  LAT(*) R*4 [I] Latitude  in input units.ACC  X(*)   R*4 [O] Position relative to map origin, to the right.7CC  Y(*)   R*4 [O] Position relative to map origin, up.CCC	if a point to be converted is outside the range of a projection,&CC	a large Y value (1.E9) is returned.KCC    If input and output are the same type, a single area can be used for BCC    both. i.e., caller may use same arrays for input and output.&C_Bugs	!!! = code needs to be checked.9C	<<< = code depends upon type specificaton for argumentsFC_Calls	COCOCS COCOSM COCEMC VNEG VNORM VUNIT VCROSS TRANS3 VSUB VDOT C	ROTATE!C_Key	MAP CARTOGRAPHIC PROJECTION>C_Hist 85feb27	HHKieffer  U.S. Geological Survey, Flagstaff AzEC	85oct16 change argument order in vector calls to use *SP conventionC_END5	PARAMETER (NUMPRJ=7)	!Number of defined projections.	LOGICAL SOUTH*	REAL*4  BXX(3),VXX(3),BVXX(3),XX(3),VB(9)B        REAL*4 PROJ(NUMPRJ)		!4 character name of the projections.	INTEGER*4 LLCODE(2,NUMPRJ) BC......	DATA PI , R45/3.1415926,.785398/ !Pi and pi/4 =45 degrees.(	DATA YERR/1.E9/			 !Invalid point flag.7        DATA LLCODE /1,1, 1,2, 1,2, 1,1, 3,1, 3,3, 1,1/&        DATA PROJ/  'MERC',	! mercator.     1	            'LAMB',	! lambert conformal1     1              'POLA',	! polar stereographic0     1              'SIMP',	! simple cylindrical(     1              'SINU',	! sinusoidal2     1	            'VIEW',     ! point perspective+     1	            'NULL'/	! no projection.HC Define conformal function as tangent of 1/2 the conformal co-latitude.MC this can safely be evaluated at all positive latitudes, including the pole.I	CONF(T,ES) = TAN(R45-0.5*T) * ((1.+ES)/(1.-ES))**HE	! conformal function3	NP = N		! transfer first argumnet to local storageA        IF (NP) 10,5,200	! Branch to: initialize, report, convertB                                                                  OC======================================================== report ==============25	CONTINUE		! output a report of input parameters.
	X(1) = BODYA		X(2) = E	X(3) = RINU	X(4) = UOUT
	X(5) = SCALE	X(6) = VXX(1)	X(7) = VXX(2)	X(8) = VXX(3)	Y(1) = PROJ(KPROJ)C        LON(1)=  DFLOTJ(LLCODE(1,KPROJ))	! <<< return line type forF        LAT(1)=  DFLOTJ(LLCODE(2,KPROJ))	! <<< meridians and latitudes	GOTO 9MC========================================= initiate projection ==============
10	KPROJ= -NPH	IF (KPROJ.GT.NUMPRJ) KPROJ = 4	!!! invalid projection request, use SIMP3	BODYA = X(1)		! body major radius (commonly in km)	E     = X(2)		! eccentricity0	RINU  = X(3)		! input angular units, in radians/	UOUT  = X(4)		! output units (commonly inches)"	SCALE = X(5)		! inverse map scale*	VXX(1)= X(6)		! central meridian   or   X*	VXX(2)= X(7)		! 1st std parallel   or   Y*	VXX(3)= X(8)		! 2nd std parallel   or   ZNC - - - - - all set-up code below here independant of argument transfer method1	PI = ACOS(-1.)		! use full precision of hardware	TWOPI = 2.0*PI
	R45=PI/4.0	RAD = PI/180.		! multiplier, degrees to radians	IF (RINU.EQ.0.)RINU=RAD	IF (UOUT.EQ.0.) UOUT=1..	DEG= RINU/RAD		! multiplier, input to degrees	HE = 0.5*E		! eccentricity/23	HEMI = VXX(2)		! save as sign of allowed latitudes	SOUTH = .FALSE.	IF (HEMI.LT.0.) SOUTH = .TRUE.<	CLON=RAD*VXX(1)		! central longitude in radians (if needed).	S = BODYA/(UOUT*SCALE)	! basic scaling factor/	F=S			! most common scaling factor for radians8	FI=S*RINU		!   "    "       "       "   for input units@	GOTO (199,120,130,199,199,160,199),KPROJ	! branch to projection9120	CONTINUE		! lambert conformal, two standard parallels.	RP1 = RAD*ABS(VXX(2))	! 1st standard parallel.	RP2 = RAD*ABS(VXX(3))	! 2nd standard parallel	ESP1 = E*SIN(RP1)	ESP2 = E*SIN(RP2) 	EM1 = COS(RP1)/SQRT(1.-ESP1**2) 	EM2 = COS(RP2)/SQRT(1.-ESP2**2)	TZ1 = CONF(RP1,ESP1)	TZ2 = CONF(RP2,ESP2)G	EL = ALOG(EM1/EM2) / ALOG(TZ1/TZ2) ! also = arcsin of tangent latitude	F = S * EM1/(EL*TZ1**EL)D!C.....	B0=F*CONF(EL,E*EL)**EL	! radius at origin (tangent) latitude		GOTO 199.130	RP1 = RAD*ABS(VXX(2))	! polar-sterographic	IF (RP1.EQ.0) RPI=0.5*PI5	IF (RP1.GT.0.49*PI) THEN	!Test for P1 near the pole.1		F = -2.*S * SQRT(1.-E**2) * ((1.-E)/(1.+E))**HE		    ELSE-		ESP1 = E*SIN(RP1)		! latitude of true scaleA		EM1 = COS(RP1)/SQRT(1.-ESP1**2)	! follows snyder eq. 17-34, but=		F = -S * EM1 / CONF(RP1,ESP1)	!   this blows up at the pole
	    ENDIF		GOTO 1991160	CALL VNORM  (VXX, VB(7))		! point perspective	CALL COCOCS (VXX,  P,Q,R)	CALL COCOSM (P,Q,  DLAT,DLON)&D		WRITE(4,*)' VIEWER',P,Q,R,DLAT,DLOND		CALL VPRINT('VB7',VB(7),4)7C should include here allowance for x-axis rotation !!!D	CALL VNEG   (VB(7), VB(7))	! construct rotation matrix as transform$	CALL VUNIT  (3, XX)			! unit z axis7	CALL VCROSS (VB(7),XX, VB(1))		! new_x = old_z x new_zD		CALL VPRINT('VB1',VB,4)6	CALL VNORM  (VB(1), VB(1))		! have x unit toward east1	CALL VCROSS (VB(1),VB(7), VB(4))	! y = z cross xD		CALL MPRINT('BV',VB,4)@	CALL TRANS3 (VB, VB)		! vb = rotation matrix from body to view.D		CALL MPRINT('VB',VB,4)	F = 1./SCALE	!!!
199	RETURNCOC================================================ do conversions ==============,200	GOTO (210,220,230,240,250,260,270),KPROJC@C __ check for longitude range; common to several projections___0205	RLON = CLON - RINU*LON(I)	! positive azimuth	IF (RLON.LT.-PI) THEN		RLON = RLON + TWOPI	    ELSEIF (RLON.GT.PI) THEN		RLON = RLON - TWOPI
	    ENDIF!	GOTO (212,222,  9,242,252),KPROJCLC------------------------------------------- mercator section --------------
210	DO I=1,NP	  PLAT=LAT(I) 	  RLAT=RINU*ABS(PLAT)	  ESIN=E*SIN(RLAT)@	  GOTO 205	! rlon = long. east from central meridian in radians212	  X(I)=F * RLON?	  Y(I)=F * ALOG( TAN(R45+0.5*RLAT) * ((1-ESIN)/(1+ESIN))**HE )	  IF (PLAT.LT.0.) Y(I)=-Y(I)	  ENDDO	GOTO 9MC---------------------------------------- lambert conformal section   -------
220	DO I=1,NP	  PLAT = LAT(I)	  IF(PLAT*HEMI.LT.0.) THEN		Y(I)=YERR		    ELSE	  RLAT = RINU * ABS(PLAT)	  ESIN = E*SIN(RLAT)	  B = F*CONF(RLAT,ESIN)**EL@	  GOTO 205	! rlon = long. east from central meridian in radians222	  QLON = EL * RLON	  X(I) = B * SIN(QLON)0	  Y(I) = -B * COS(QLON)		! y origin at the pole=!C..	  Y(I) = B0 - B*COS(QLON)	! y origin at tangent latitude	  IF (SOUTH) Y(I)=-Y(I)
	    ENDIF	  ENDDO	GOTO 9MC---------------------------------------- polar stereographic section -------
230	DO I=1,NP	  IF(LAT(I)*HEMI.LT.0.) THEN		Y(I)=YERR		    ELSE	  PLAT=LAT(I)	  RLAT=RINU*ABS(PLAT)	  ESIN=E*SIN(RLAT)	  B= F * CONF(RLAT,ESIN)	  RLON=RINU*LON(I)	  X(I)=B*SIN(RLON)	  Y(I)=B*COS(RLON)	  IF (SOUTH) Y(I)=-Y(I)
	    ENDIF	  ENDDO	GOTO 9NC---------------------------------------- simple cylindrical section ---------
240	DO I=1,NP@	  GOTO 205	! rlon = long. east from central meridian in radians242	  X(I)=F * RLON	  Y(I)=FI * LAT(I)	  ENDDO	GOTO 9NC---------------------------------------- sinusoidal section -----------------250	DO I=1,NP			@	  GOTO 205	! rlon = long. east from central meridian in radians+252	  X(I)=F * RLON * COS(RINU*ABS(LAT(I)))	  Y(I)=FI * LAT(I)	  ENDDO	GOTO 9NC--------------------------------------- point perspective section -----------
260	DO I=1,NP,C get cartesian coordinates of surface point4	  CALL COCEMC (DEG*LAT(I),DEG*LON(I),BODYA,E,  BXX);	  CALL VSUB (BXX,VXX, BVXX)	! offset vector in body system9	  CALL VDOT (BXX,BVXX, SIGN)	! -cosine of emission angle2	  IF (SIGN.GT.0.) THEN		! point is hidden by body		Y(I)=YERR		    ELSE>		CALL ROTATE (VB,BVXX, XX) ! system with z axis to map center3		B = F/XX(3)		  ! scale back onto projection plane		X(I) = B * XX(1)			Y(I) = B * XX(2)
	    ENDIF	  ENDDO	GOTO 9OC--------------------------------------- null projection section --------------
270	DO I=1,NP	  X(I) = FI * LON(I)	  Y(I) = FI * LAT(I)	ENDDO
      	GOTO 99	RETURN        END