; $Id: map_proj_info.pro,v 1.1 2004/06/16 20:43:03 chris Exp $
;
; Copyright (c) 2004-2004, Research Systems, Inc.  All rights reserved.
;   Unauthorized reproduction prohibited.

; Undocumented helper routines needed for
; MAP_SET, MAP_GRID, MAP_CONTINENTS.
;
; CT, RSI, June 2004: Split off from map_set.pro

;------------------------------------------------------------------------
; Return information about maps in general, or projection iproj.
; Iproj = projection number, Input.  If CURRENT is set, then the
;   current map projection is returned in Iproj.
; CURRENT = if set, use current proj number, and return that index in iproj.
; UV_RANGE =  UV limits of the given projection, [umin, vmin, umax,
;   vmax].  Also keyword UVRANGE is accepted here for bkwds compat.
; SCALE = (output) the number of meters to one unit of UV space.  The
;     radius of the earth, from Clarke's 1866 ellipsoid is used (6378206).
; NAME = the name of the projection.
; CYLINDRICAL = 1 if the projection is cylindrical or pseudo-cylindrical.
; CONIC = 1 if the projection is conic.
; PROJ_NAMES = string array containing the names of all the
;   available projections.
; UV_LIMITS = returned rectangle of the current map in UV coordinates,
;   [Umin, Vmin, Umax, Vmax].
; LL_LIMITS = returned rectangle of the current map in latitude/longitude
;   coordinates, [Lonmin, Latmin, Lonmax, Latmax].  Note, the
;   range may not always be available, and the min and max may
;   both be set to 0.
;
; Issues for a new projection:
; UVrange, proj_scale, is_cyl, is_azimuthal, circle, name, horizon, clipping
;
PRO map_proj_info, iproj, CURRENT=current, $
       UV_RANGE=r, UVRANGE=rdummy, SCALE=meters, NAME=name, CYLINDRICAL=cyl, $
       PROJ_NAMES = proj_names_p, CIRCLE=is_circular, CONIC=is_conic, $
       AZIMUTHAL=is_azimuthal, UV_LIMITS = uv_limits, $
       LL_LIMITS = ll_limits

    COMPILE_OPT hidden


    common map_set_common, proj_scale, proj_is_cyl, proj_is_conic, $
      proj_names, circle, proj_is_azim, user_defined

    on_error, 2
    if n_elements(proj_scale) eq 0 then begin ;Init constant arrays?
        p2 = 2 * !pi
        s8 = 2 * sqrt(8)
        ro = 2 * 2.628

    ;      St Or LC La Gn  Az Sa  Cy  Me  Mo  Si  Ai  Ha Al Tm  Mi Ro Le Go
    ;       1  2  3  4  5   6  7   8   9  10  11  12  13 14 15  16 17 18 19
      proj_scale =  $
        [0, 4, 2, 4, 4, 4, p2, 1, p2, p2, s8, p2, p2, s8, 4, 0, p2,ro, 0,p2]
      proj_is_cyl = $
        [0, 0, 0, 0, 0, 0,  0, 0,  1,  1,  1,  1,  1,  1, 0, 1,  1, 1, 0, 1]
      proj_is_azim = $
        [0, 1, 1, 0, 1, 1,  1, 0,  0,  0,  0,  0,  0,  0, 0, 0,  0, 0, 0, 0]
      proj_is_conic = $
        [0, 0, 0, 1, 0, 0,  0, 0,  0,  0,  0,  0,  0,  0, 1, 0,  0, 0, 1, 0]
      circle = $
        [0, 1, 1, 0, 1, 1,  1, 1,  0,  0,  1,  0,  1,  1, 0, 0, 0,  0, 0, 0]
      proj_names =  $
        ["", "Stereographic", "Orthographic", "LambertConic", "LambertAzimuthal", $
         "Gnomic", "AzimuthalEquidistant", "Satellite", "Cylindrical", $
         "Mercator", "Mollweide", "Sinusoidal", "Aitoff", "HammerAitoff", $
         "AlbersEqualAreaConic", "TransverseMercator", $
         "MillerCylindrical", "Robinson", "LambertConicEllipsoid", $
         "GoodesHomolosine"]

      user_defined = bytarr(n_elements(proj_names))
    endif

    max_proj = n_elements(proj_names)-1
    if arg_present(proj_names_p) then proj_names_p = proj_names
    if KEYWORD_SET(current) then iproj = !map.projection ;Get current projection

    if ARG_PRESENT(uv_limits) then uv_limits = !map.uv_box
    if ARG_PRESENT(ll_limits) then ll_limits = !map.ll_box

    if n_elements(iproj) eq 0 then return ;Rest of params are proj specific

    if iproj lt 1 or iproj gt max_proj then $
      message, 'Projection number must be within range of 1 to'+ strtrim(max_proj,2)

    ones = [-1., -1., 1., 1.]
    r_earth = 6378206.4d0 ;Earth equatorial radius, meters, ; Clarke 1866 ellipsoid

    cyl = proj_is_cyl[iproj]
    name = proj_names[iproj]
    meters = proj_scale[iproj] * r_earth ;Scale [= 0 if 1:1]
    is_circular = circle[iproj]
    is_conic = proj_is_conic[iproj]
    is_azimuthal = proj_is_azim[iproj]

    if arg_present(r) or arg_present(rdummy) or arg_present(meters) then begin
        case iproj of
        1: r = 2 * ones             ;Stereo
        2: r = ones                 ;Ortho
        3: r = 2 * ones             ;Lambert Conic for sphere.
        4: r = 2 * ones             ;lambert
        5: r = 2 * ones             ;gnomic
        6: r =  !DPI * ones         ;azimuth
        7: if !map.p[1] ne 0 then begin ;Complicated satellite?
            map_satellite_limit, 60, xr, yr
            if n_elements(xr) eq 0 then r = [0,0,0,0] $
            else r = [min(xr, max=maxxr), min(yr, max=maxyr), maxxr, maxyr]
            meters = (r[2]-r[0]) * r_earth
        endif else begin
            st = !map.p[0] > 1.0   ;For satellite
            r = sqrt((st-1)/(st+1)) * ones ;Simple case satellite
            meters = (r[2]-r[0]) * r_earth
        endelse
        8: r = [-1., -.5, 1., .5] * !DPI ;Cylind equidistant
        9: r = [-!DPI, -2.43, !DPI, 2.43] ;Mercator  (+/- 80 degs of lat) ~ 98%
        10: r = sqrt(2.0) * [-2.,-1.,2.,1.] ;Mollweide
        11: r = [-1., -.5, 1., .5] * !DPI ;sinusoidal
        12: r = [-1., -.5, 1., .5] * !DPI ;Aitoff
        13: r = sqrt(2.0) * [-2., -1., 2., 1.] ; Hammer-aitoff
        14: r = 2 * ones            ;Albers Conic
        15: r = 2*!pi * r_earth * ones ;Transverse Mercator, Units = meters.
        16: r = [-!DPI, -1.832, !DPI, 1.832] ;Miller Cylind  (+/- 80 degs of lat)
        17: r = [-2.628, -1.349, 2.628, 1.349] ;Robinson
        18: r = 2 * !pi * r_earth * ones ;Lambert Conic ellipsoid, Units = meters
        19: r = [-!DPI, -1.37, !DPI, 1.37] ;Goodes Homolosine
    endcase

    rdummy = r                      ;For backwards compatibility
    if meters ne 0 then meters = meters / (r[2]-r[0]) $ ;Get scale in UV units
    else meters = 1.0               ;UV Units are meters.
    endif

end