;+
; $Id: euvi_heliographic.pro,v 1.3 2011/08/26 22:02:59 nathan Exp $
; Project     :	STEREO - SECCHI
;
; Name        :	EUVI_HELIOGRAPHIC()
;
; Purpose     :	Returns heliographic coordinates for EUVI images.
;
; Category    :	STEREO, SECCHI, Coordinates
;
; Explanation :	This routine uses the information in an EUVI FITS header to
;               calculate the heliographic longitude and latitude.
;
; Syntax      :	COORD = EUVI_HELIOGRAPHIC(HEADER  [, PIXELS] [, /CARRINGTON] )
;
; Examples    :	COORD = EUVI_HELIOGRAPHIC(HEADER)
;               LONGITUDE = COORD[0,*,*]
;               LATITUDE  = COORD[1,*,*]
;
; Inputs      :	HEADER  = The FITS header structure.-OR- WCS structure
;
; Opt. Inputs :	PIXELS  = If passed, then contains an array of IDL pixel
;                         locations to return coordinates from.  Otherwise,
;                         coordinates are returned for all the pixels in the
;                         array. The first dimension must be 2.
;
; Outputs     :	The result of the function is an array containing the longitude
;               and latitude in degrees for each pixel.  Points off the limb
;               will be marked as Not-a-Number (NaN).
;
; Opt. Outputs:	None.
;
; Keywords    :	CARRINGTON = If set, then Carrington coordinates are return.
;                            The default is to return Stonyhurst coordinates.
;
; Calls       :	FITSHEAD2WCS, WCS_GET_COORD, FLAG_MISSING
;
; Common      :	None.
;
; Restrictions:	None.
;
; Side effects:	None.
;
; Prev. Hist. :	None.
;
; History     :	Version 1, 17-Mar-2008, William Thompson, GSFC
;
; Contact     :	WTHOMPSON
;-
; $Log: euvi_heliographic.pro,v $
; Revision 1.3  2011/08/26 22:02:59  nathan
; use solar_b0 not hglt_obs to find b0; still will not work for AIA non-carrington
;
; Revision 1.2  2008/03/21 14:06:59  nathan
; fixed case where none are missing; allow wcs header input
;
; Revision 1.1  2008/03/17 19:32:19  nathan
; *** empty log message ***
;
function euvi_heliographic, header, pixels, carrington=carrington, WCS=wcs
;
;  Generate a WCS structure from the header, and extract the coordinates.  Save
;  the dimensions for later.
;
IF tag_exist(header,'WCSNAME') THEN wcs=header ELSE wcs = fitshead2wcs(header)
szp=size(pixels)
IF szp[1] NE 2 THEN coord = wcs_get_coord(wcs) ELSE coord = wcs_get_coord(wcs,pixels)
sz = size(coord)
dim = sz[1:sz[0]]
coord = reform(coord, 2, n_elements(coord)/2, /overwrite)
;
;  Set up some constants, based on the header.
;
deg2rad = !dpi / 180.d0                 ;degrees to radians
sec2rad = deg2rad / 3600.d0             ;arcseconds to radians
rsun = 6.95508d8                        ;meters
dsun = wcs.position.dsun_obs
b0 = wcs.position.solar_b0 * deg2rad
cos_b0 = cos(b0)
sin_b0 = sin(b0)
if keyword_set(carrington) then phi0 = wcs.position.crln_obs else $
  phi0 = wcs.position.hgln_obs
;
;  Calculate the distance of each point from the observer, assuming a distance
;  of one solar radius from sun center.  First, calculate the part inside the
;  square root--points where this is negative are off the limb and are flagged
;  as missing.
;
theta_x = reform(coord[0,*]) * sec2rad
theta_y = reform(coord[1,*]) * sec2rad
cos_x = cos(theta_x)  &  sin_x = sin(theta_x)
cos_y = cos(theta_y)  &  sin_y = sin(theta_y)
ds = dsun^2 * (cos_y^2 * cos_x^2 - 1) + rsun^2
offlimb = where(ds lt 0, noff)
d = dsun * cos_y * cos_x - sqrt(abs(ds))

IF noff GT 0 THEN flag_missing, d, offlimb
;
;  Convert from helioprojective to heliocentric coordinates using equations 15
;  from Thompson (2006).
;
x = d * cos_y * sin_x
y = d * sin_y
z = dsun - d * cos_y * cos_x
;
;  Convert from heliocentric-cartesian to heliographic, using equations 12 from
;  Thompson (2006).  Make sure that the longitude goes from 0 to 360.
;
theta = asin((y*cos_b0 + z*sin_b0) / rsun) / deg2rad
phi = (atan(x, z*cos_b0 - y*sin_b0) / deg2rad  +  phi0) mod 360.d0
w = where(phi lt 0, count)
if count gt 0 then phi[w] = phi[w] + 360.d0
;
;  Store in the coordinates array, and restore the original dimensions.
;
;stop
coord[0,*] = phi
coord[1,*] = theta
coord = reform(coord, dim, /overwrite)
return, coord
;
end