function hpc2hpr, cin, DEGREES=degrees, SILENT=silent, OLD=old
;+
; $Id: hpc2hpr.pro,v 1.3 2010/09/30 22:00:34 nathan Exp $
;
; Project   : STEREO SECCHI
;                   
; Name      : hpc2hpr
;               
; Purpose   : convert cartesian coordinates to radial coordinates
;               
; Explanation: 
;               
; Use       : IDL> result= hpc2hpr( [34,123.3])
;    
; Inputs    : 
;   	cin 	DBLARR	[2,n] or [2,nx,ny] (x is longitude and y is latitude)
;
; OPtional inputs:
;               
; Outputs   : 
;   	coords	DBLARR of same dimension as cin in units degrees (1st dimension is position angle, 2nd dim is elongation)
;
; Keywords  : 
;   	/DEGREES     cin is in units degrees; otherwise assumes arcsec. Also 
;   	    	    Does error checking for max/min degrees of 180
;   	/SILENT     Do not print result
;               
; Calls from LASCO : NONE
;
; Common    : 
;               
; Restrictions: 
;               
; Side effects: prints results if just 1 point
;               
; Category    : astrometry wcs coordinates mapping 
;               
; Prev. Hist. : None.
;
; Written     : N.Rich, I2, 6/2008
;               
; $Log: hpc2hpr.pro,v $
; Revision 1.3  2010/09/30 22:00:34  nathan
; added /OLD; reverse order of output to match WCS convention
;
; Revision 1.2  2008/06/26 18:16:50  nathan
; used formula for spherical coords from W.Thompson email
;
; Revision 1.1  2008/06/25 18:20:28  nathan
; from pix2coord.pro
;
;-            

xy=double(cin)
; xy is [2,n] or [2,nx,ny], units of degrees or arcsec
szo=size(xy)

rx=reform(xy[0,*,*])	;
ry=reform(xy[1,*,*])

IF keyword_set(DEGREES) THEN BEGIN
    rx0=rx
    ry0=ry
    wok=where(rx GE -180 and rx LT 180 and ry GE -180 and ry LT 180, nok)
    IF ~keyword_set(SILENT) THEN print,trim(nok),' points are in range.'
    rx=rx[wok]*!pi/180
    ry=ry[wok]*!pi/180
ENDIF ELSE BEGIN
    rx=rx*(!pi/180)/3600
    ry=ry*(!pi/180)/3600
ENDELSE

IF keyword_set(OLD) THEN BEGIN
rsec=sqrt(rx^2+ry^2)*180./!pi

;r = rpix * arcsc/ solrc   	; Rsun
;IF r LE 1.0 and elong EQ 2 THEN BEGIN
;    ;--Display Carrington latitude and longitude
;    cl=euvi_heliographic(wcs,[ev.x,ev.y],/CARRINGTON)
;    r=	    cl[0]   ; longitude
;    angle=  cl[1]   ; latitude
;    pfx='CLN='
;    sfx='!9%!3'
;    pfxh='CLT='
;    fonttype=-1   

;ENDIF ELSE $
;--Compute latitude angle

angle=-atan(-rx,-ry)*180./!pi
ENDIF ELSE BEGIN
;
;print,'% PLANE_PROJECT: spherical'
;radial = acos(cos(rx)*cos(ry))
radial = atan(sqrt(cos(ry)^2*sin(rx)^2 + sin(ry)^2), cos(ry)*cos(rx))
;azim   = atan(sin(ry), cos(ry)*sin(rx))
azim   = atan(cos(ry)*sin(rx), -sin(ry))
;x_proj = radial * cos(azim)
;y_proj = radial * sin(azim)
;rad = minmax(radial)

rsec=radial*180/!pi
angle=azim*180/!pi

;IF (min(angle) lt 0) THEN 
;
ENDELSE
angle=angle+180.
;stop
IF szo[0] EQ 1 THEN BEGIN
    coords=[angle,rsec]
    IF ~keyword_set(SILENT) THEN BEGIN
    	print, coords
	wait,2
    ENDIF
ENDIF ELSE BEGIN
    IF keyword_set(DEGREES) THEN BEGIN
    	; make sure output has same dim as input
	; and mark missing data as NaN
	rsec0=rsec
	angle0=angle
	rsec=float(rx0)
	angle=ry0
	rsec[*]=!values.d_nan
	angle[*]=!values.d_nan
	rsec[wok]=rsec0
	angle[wok]=angle0
    ENDIF
    ;for i=0,szin[2] do print,rsec[i],angle[i]
    IF szo[0] EQ 2 THEN BEGIN
	coords=make_array(szo[1:2],/double)
	coords[1,*]=rsec
	coords[0,*]=angle
    ENDIF ELSE BEGIN
	coords=make_array(szo[1:3],/double)
	coords[1,*,*]=rsec
	coords[0,*,*]=angle
    ENDELSE    	
ENDELSE

return,coords
END