function pix2coord, hdr0, pix, DEBUG=debug, SYSTEM=system, SPACECRAFT=spacecraft, SPHERICAL=spherical, $
    	SOLAR_RADII=solar_radii
;+
; $Id: pix2coord.pro,v 1.10 2009/10/01 20:05:59 nathan Exp $
;
; Project   : STEREO SECCHI
;                   
; Name      : pix2coord
;               
; Purpose   : return coordinate of pixel(s) of a SECCHI image header in chosen system
;               
; Explanation :	This routine converts the native solar coordinates of a FITS
;               file using the World Coordinate System into another solar
;               coordinate system.  The supported coordinate systems are:
;
;                       GEI     Geocentric Equatorial Inertial = RA/Dec
;                       GEO     Geographic
;                       GSE     Geocentric Solar Ecliptic
;                       MAG     Geomagnetic
;                       GSM     Geocentric Solar Magnetospheric
;
;                       Helioprojective-Cartesian       (HPC)
;                       Helioprojective-Radial          (HPR)
;                       Heliocentric-Cartesian          (HCC)
;                       Heliocentric-Radial             (HCR)
;                       Stonyhurst-Heliographic         (HG)
;                       Carrington-Heliographic         (CAR)
;   	    	    	Heliocentric Aries Ecliptic 	(HAE)
;   	    	    	HAE Spherical  	    	    	(HAE, /SPHER)
;
;               The helioprojective coordinates are angular position as seen by
;               the observer, while heliocentric coordinates are the
;               equivalents in spatial units (e.g. kilometers). Spherical coordinates
;   	    	are computed assuming a constant radius, specified as a multiple of 
;   	    	the mean solar radius.
;               
; Use       : IDL> result= pix2coord( hdr, [34,123.3])
;    
; Inputs    : 
;   	hdr 	FITS header structure OR WCS structure
;
; OPtional inputs:
;   	pix 	FLTARR  [2,n1,n2] array of IDL pixel coordinates; if not present, coordinates
;   	    	of entire array for input header will be returned.  	
;               
; Outputs   : 
;   	coords	DBLARR dimension depends on SYSTEM:
;   	    projective/spherical: [[pos angle, elongation] or [longitude, latitude],nx,ny]
;   	    cartesian: [[x,y,z],nx,ny]
;
; Keywords  : 
;   	SYSTEM=  STRING designation of desired system; if not specified, returns
;   	    	 helio-projective radial (HPR)
;   	/DEBUG	Print pixel location of sun center for input image and some other info
;   	SPACECRAFT= String specifying STEREO-A or B. Required for HAE cartesian and input is not a FITS header.
;   	/SPHERICAL  For HAE, return spherical coordinates at SOLAR_RADII distance from sun center
;   	SOLAR_RADII=	Specify solar radii to compute HAE spherical coordinate; 
;   	    Currently works only for 1.
;               
; Calls     : FITSHEAD2WCS, WCS_GET_COORD, WCS_CONV_FIND_DSUN, WCS_CONVERT_FROM_COORD,
;   	      CONVERT_STEREO_COORD	
;
; Common    : 
;               
; Restrictions: HAE coordinates require use ofthe Icy/CSPICE package; see convert_stereo_coord.pro
;               
; Side effects: 
;               
; Category    : Astrometry, coordinates, wcs
;               
; Written     : N.Rich, I2, 6/2008
;               
; $Log: pix2coord.pro,v $
; Revision 1.10  2009/10/01 20:05:59  nathan
; add geocentric systems (only tested gei)
;
; Revision 1.9  2009/09/29 17:25:45  nathan
; if SOHO, load_soho_spice
;
; Revision 1.8  2009/09/22 18:37:23  nathan
; fix typo in input
;
; Revision 1.7  2009/09/09 21:37:18  nathan
; minor changes to accomodate non-secchi headers and strarr headers
;
; Revision 1.6  2009/08/14 19:25:13  nathan
; change computation of HAE /spherical as recommended by W.Thompson
;
; Revision 1.5  2009/08/12 21:12:00  nathan
; update keyword info
;
; Revision 1.4  2009/08/12 20:58:57  nathan
; Complete re-write to accomodate multiple coordinate systems; utilizes
; WCS routines instead of local hpc2hpr.pro.
;
; Revision 1.3  2008/06/25 18:23:45  nathan
; use hpc2hpr.pro
;
; Revision 1.2  2008/06/24 15:17:56  nathan
; fix comments
;
; Revision 1.1  2008/06/23 17:27:30  nathan
; computes r-theta of input
;
;-            
IF keyword_set(SOLAR_RADII) THEN BEGIN
    IF solar_radii NE 1. THEN BEGIN
    	message,'Currently only SOLAR_RADII=1. is allowed; continuing...',/info
	wait,5
    ENDIF
ENDIF
sc=''
IF datatype(hdr0) NE 'STC' THEN hdr=fitshead2struct(hdr0) ELSE hdr=hdr0
IF tag_exist(hdr,'COORD_TYPE') THEN wcs=hdr ELSE BEGIN
    wcs=fitshead2wcs(hdr)
ENDELSE
IF keyword_set(SPACECRAFT) THEN sc=strupcase(spacecraft) ELSE $
IF tag_exist(hdr,'OBSRVTRY') THEN sc=hdr.obsrvtry ELSE $
IF tag_exist(hdr,'TELESCOP') THEN sc=hdr.telescop ELSE sc=''
IF tag_exist(wcs.time,'OBSERV_AVG') THEN dobs=wcs.time.observ_avg ELSE dobs=wcs.time.observ_date
darr=0
final=''    ; set this if there is an intermediary step between HPC and final system
docarr=0
doradial=0

IF keyword_set(DEBUG) THEN debugon=1 ELSE debugon=0
;	rpix=rsec/wcs.cdelt[0]
;	IF cunit EQ 'deg' THEN unitf=3600. ELSE unitf=1.

cunit=wcs.cunit[0]

xy=wcs_get_coord(wcs,pix)   ; HPC coordinates

;coords = hpc2hpr(xy, DEGREES=(cunit EQ 'deg'))	

WCS_CONV_FIND_DSUN, DSUN, dRSUN, WCS=WCS

IF (debugon) THEN BEGIN
    help,xy,dsun,drsun
    ; get sun center
    cxy=wcs_get_pixel(wcs,[0,0])
    cx=cxy[0]
    cy=cxy[1]
    print,'Sun center should be:',cx,cy,' pixels'
    
ENDIF


IF keyword_set(SYSTEM) THEN BEGIN
    sysout=strupcase(strmid(system,0,3))
;--First check for earth-centric systems and do them

    w=where(sysout EQ ['GEI','GEO','GSE','MAG','GSM'],isearthcentered)
    IF (isearthcentered) THEN BEGIN
	IF sc EQ '' THEN BEGIN
	    message,'SPACECRAFT= must be set for '+sysout+' coordinate system.',/info
    	    return,-1
	ENDIF
	if cunit EQ 'arcsec' THEN divi=3600. ELSE divi=1.
	coords=xy/divi
    	message,'Calling SPICE/convert_stereo_lonlat',/info
    	IF sc EQ 'SOHO' THEN load_soho_spice
    	convert_stereo_lonlat,dobs,coords,'HPC',sysout,/degrees,/pos_long,spacecraft=sc
	return,coords
    ENDIF
    CASE sysout OF
    	'HAE' : BEGIN
    	    IF keyword_set(SPHERICAL) THEN sysout='HG' ELSE sysout='HCC'
	    final='HAE'
    	    END
	'HPR' : 
	'HCR' :
	'HPC' : return,xy
	'HCC' : 
	'HG'  :
	'CAR' : BEGIN
	    sysout='HG'
	    docarr=1
	    END
	
        else: message, 'Coordinate type ' + sysout + ' not supported'
    endcase
ENDIF ELSE BEGIN
    sysout='HPR' ; default
    message,'Returning HPR coordinates.',/info
ENDELSE


WCS_CONVERT_FROM_COORD,wcs,xy,sysout,x,y,z  ,/zero_center, /pos_long, $
    arcseconds=(wcs.cunit[0] eq 'arcsec'), CARRINGTON=docarr
IF sysout EQ 'HG' THEN undefine,z

IF debugon THEN help,sysout,x,y,z,xy
sz=size(xy)

; xy is arr[2] or arr[2,i,j]

IF sz[0] EQ 1 AND datatype(z) EQ 'UND' THEN coords=[x,y]
IF sz[0] EQ 1 AND datatype(z) NE 'UND' THEN coords=[x,y,z]
IF sz[0] EQ 3 and datatype(z) EQ 'UND' THEN BEGIN
    coords=dblarr(2,sz[2],sz[3])
    coords[0,*,*]=x
    coords[1,*,*]=y
ENDIF
IF sz[0] EQ 3 and datatype(z) NE 'UND' THEN BEGIN
    coords=dblarr(3,sz[2],sz[3])
    coords[0,*,*]=x
    coords[1,*,*]=y
    coords[2,*,*]=z
ENDIF


IF final NE '' THEN BEGIN
    ; need to convert coordinates using SPICE with HAE
    IF keyword_set(SPHERICAL) THEN Begin
   	convert_stereo_lonlat,dobs,coords,'HEEQ', 'HAE',/degrees,/pos_long,spacecraft=sc
    Endif ELSE BEGIN
	IF sc EQ '' THEN BEGIN
	    message,'SPACECRAFT= must be set for '+final+' coordinate system.',/info
    	    return,-1
	ENDIF
    	message,'Calling SPICE/convert_stereo_coord',/info
	;--I don't think this is ever called if Z is undefined!
	coords[0,*,*]=z
	coords[1,*,*]=x
	coords[2,*,*]=y
    	IF sc EQ 'SOHO' THEN load_soho_spice
	convert_stereo_coord, dobs,coords,'HGRTN',final,spacecraft=sc
    ENDELSE

    ; (STUB)
    IF (0) THEN BEGIN
    	message,'Now converting from cartesian to spherical coordinates.',/info
	x=reform(coords[0,*,*])
    	y=reform(coords[1,*,*])
    	z=reform(coords[2,*,*])
	rho=1.
	S=sqrt(x^2+y^2)
	lat=90.-(180./!pi)*acos(z/drsun)
	xltz=where(x LT 0)
	xgez=where(X GE 0)
	lon=theta
	lon[xltz]=(180./!pi)*asin(y[xltz]/S[xltz])
	lon[xgez]=180.-(180./!pi)*asin(y[xgez]/S[xgez])
	lon=lon+90.
	pltz=where(lon lt 0,npltz)
	IF npltz GT 0 THEN lon[pltz]=lon[pltz]+360.
	
    	coords=dblarr(2,sz[2],sz[3])
	coords[0,*,*]=lat
	coords[1,*,*]=lon
	; lon=phi, lat=theta
    ENDIF
    
ENDIF

return,coords
END