; docformat = 'rst'
;
; NAME: 
;    areaSubScript
; PURPOSE:
;    Determine the subscript of an area on a map given a projection geometry
;+
; :Description:
;    Determine the subscript of an area on a map given a projection geometry
;
; :Categories:
;    Mapping
; 
; :Params:
;    longitude: in, required, type=array
;      Array of longitudes defining the area
;    latitude: in, required, type=array
;      Array of latitudes defining the area
;
; :Keywords:
;    geometry: in, optional, type=string, default='ECP'
;      The geometry of mapping (ECP | ortho)
;    dimension: in, optional, type=float, default=512
;      The dimension (pixels) of the shortest axis of the map
;
; :Uses:
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in 2008, B. Carry (ESO)::
;      2018 July - B.Carry (OCA) - Complete rewriting
;-
function areaSubScript, longitude, latitude, geometry=geometry, dimension=dimension
  COMPILE_OPT hidden, idl2
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--I-1-- Check Syntax  ------------------------------------------------------------------------;
  if not keyword_set(longitude) or not keyword_set(latitude) then begin
    message, /ioError, 'You must specify both longitude and latitude'
  endif


 ;--I-2-- Test Input Coordinates  --------------------------------------------------------------;
  nLon=n_elements(longitude)
  nLat=n_elements(latitude)
  if nLon ne nLat then begin
    message, /ioError, 'The dimension of longitude and latitude must agree: '
  endif


 ;--I-3-- Default Parameters  ------------------------------------------------------------------;
  if not keyword_set(dimension) then vertAxis=512 else vertAxis=dimension
  if not keyword_set(geometry)  then geometry='ECP'




 ;-Output
  set_plot, 'Z'



;---------------------------------;
;- Orthographic Polar Projection -;
;---------------------------------;
  if geometry eq 'ortho' then begin

   ;-Map size
    horiAxis = vertAxis
    rayon = vertAxis /2.

   ;-Choice of the Pole (North-South)
    signLat = mean(latitude)/abs(mean(latitude))
    if signLat ge 0 then refLat=90 else refLat=-90

   ;-Setting the projection parameters
    map_set, /ORTHOGRAPHIC, $
             /ISOTROPIC, $
             refLat, 0, 0, $
             xMargin=0, yMargin=0, $
             /NOBORDER, $
             position=[0,0,horiAxis, vertAxis], $
             ZVALUE=0

  endif

 

;--------------------------------------;
;- Equidistant Cylindrical Projection -;
;--------------------------------------; 
  if geometry eq 'ECP' then begin

   ;-Map size
    horiAxis = 2*vertAxis

   ;-Limits
;    limitArr=[-90,0,90,360] ;-No Limit --no working very well
;    limitArr=[-85,0,85,360] ;-Without Poles
    limitArr=[-89,0,89,360] ;-Without Poles


   ;-Setting the projection parameters
    map_set, /CYLINDRICAL, $
             /ISOTROPIC, $
             0, 180, 0, $
             limit=limitARR, $
             xMargin=0, yMargin=0, $
             /NOBORDER, $
             position=[0,0,horiAxis, vertAxis], $
             ZVALUE=0

  endif

;-------------------;
;- Filling the Map -;
;-------------------;

 ;-Choosing the resolution
  device, SET_RESOLUTION=[horiAxis, vertAxis]

 ;-Plotting the Area with the given projection
  polyfill, longitude, latitude, $
            clip=[0,0,horiAxis, vertAxis], Z=0

 ;-Retrieving the position
  image=tvrd()
  area = where( image ne 0 )



 ;-Restoring the device and returning the subscript
  set_plot, 'X'
  return, area

end