PDS_VERSION_ID                   = PDS3

RECORD_TYPE                      = STREAM

OBJECT                            = DATA_SET_MAP_PROJECTION
 DATA_SET_ID                      = "MRO-M-SHARAD-5-3D-V1.0"

 OBJECT                           = DATA_SET_MAP_PROJECTION_INFO
  MAP_PROJECTION_TYPE             = "POLAR STEREOGRAPHIC"
  MAP_PROJECTION_DESC             = "

In this polar stereographic projection [SNYDER1987], the projection is
centered on the north or south pole. Lines of longitude extend
radially from the center and parallels of latitude are concentric
circles around the center. In the north, Longitude 0 extends straight
down from the center and Longitude 90 East extends to the right.
In the south, Longitude 0 extends straight up from the center, and
Longitude 90 East extends to the right. 

SHARAD observations input to the 3D processing are provided with along-track
information (see GEOM files in the SHARAD US PDS RDR archive) that includes
east longitudes and planetocentric latitudes. The 3D processing procedures
are carried out in cartesian space, so the locations for the polar data were
first transformed to a polar stereographic projection. The IDL procedures
MAP_PROJ_INIT and MAP_PROJ_FORWARD were used for this transformation, 
specifying the GCTP polar stereographic projection as well as Mars semimajor 
and semiminor axes given as A and B below. Given the planetocentric latitudes,
this ellipsoidal projection was not an ideal choice. Nevertheless, accurate
longitudes and planetocentric latitudes can be obtained by inverting the
projection using this same ellipsoid (e.g., using IDL's MAP_PROJ_INVERSE).

Thus, the transformation from line and sample coordinates to planetocentric 
latitude and longitude is 

   X = (I - N/2 - 0.5)*MAP_SCALE
   Y = (J - N/2 - 0.5)*MAP_SCALE

   R = SQRT(X^2 + Y^2)
 
   A = 3396000
   B = 3376000
   E = SQRT(1-B^2/A^2)

   S = R*SQRT(((1+E)^(1+E))*((1-E)^(1-E)))/(2*A)

   LON = (180/PI)*ATAN2(X,-Y)
   For the northern hemisphere:
   LAT =  90 - (360/PI)*ATAN(S*((1-E*SIN(LAT*PI/180))/
          (1+E*SIN(LAT*PI/180)))^(E/2))
   For the southern hemisphere:
   LAT = -90 + (360/PI)*ATAN(S*((1-E*SIN(LAT*PI/180))/
          (1+E*SIN(LAT*PI/180)))^(E/2))

where

   I = sample coordinate
   J = line coordinate
   X,Y = I,J converted to Cartesian coordinate system with 0,0 at center
   MAP_SCALE = the map resolution in meters per pixel
   N = the number of lines or samples per line in the (square) array
   R = radius from center of map to pixel I,J in meters
   A = Mars equatorial radius (semimajor axis) in meters
   B = Mars polar radius (semiminor axis) in meters
   E = eccentricity of the Mars ellipsoid
   LON = east longitude of pixel I,J in degrees
   LAT = planetocentric latitude of pixel I,J in degrees

Values for MAP_SCALE and N can be found in the data volume label as the
keywords MAP_SCALE and AXIS_ITEMS.

The solutions for latitude above require either iteration or recasting 
using a series. See [SNYDER1987] for details.

The transformation from longitude (LON) and planetocentric latitude (LAT) to 
SAMPLE and LINE indices (I,J) is

   T = SQRT((1-SIN(LAT*PI/180))/(1+SIN(LAT*PI/180))*
       ((1+E*SIN(LAT*PI/180))/(1-E*SIN(LAT*PI/180)))^E)
   P = 2*A*T/SQRT((1+E)^(1+E)*(1-E)^(1-E))
   X = P*SIN(LON*PI/180)
   Y = -P*COS(LON*PI/180)
   I = NINT(X/MAP_SCALE + N/2 + 0.5)
   J = NINT(Y/MAP_SCALE + N/2 + 0.5)

NINT is a Fortran intrinsic rounding function, identical to the C math.h
library function rint.  In the case of half-integer values, the NINT
function rounds to the nearest even whole number.
"

  ROTATIONAL_ELEMENT_DESC         = "See [SEIDELMANNETAL2002]."

  OBJECT                          = DS_MAP_PROJECTION_REF_INFO
   REFERENCE_KEY_ID               = "SEIDELMANNETAL2002"
  END_OBJECT                      = DS_MAP_PROJECTION_REF_INFO

  OBJECT                          = DS_MAP_PROJECTION_REF_INFO
   REFERENCE_KEY_ID               = "SNYDER1987"
  END_OBJECT                      = DS_MAP_PROJECTION_REF_INFO

 END_OBJECT                       = DATA_SET_MAP_PROJECTION_INFO
END_OBJECT                        = DATA_SET_MAP_PROJECTION

END