;+
; PURPOSE: transform from helioprojective cartesian to helioprojective radial. See W.T.Thompson, A&A Dec 15, 2005, section 4.
;
; INPUTS:
;   coord : helioprojective cartesian coordinates, in degrees.
;           either two element array : [westward_angle,northward_angle]
;           or array of [0:1,*], with [0,*] the westward angle
;                                     [1,*] the northward angle
;
; OUTPUTS: helioprojective radial coordinates [polar_angle,elongation], in the same format than the input.
;
;-

function heliocart2heliorad,coord

if n_elements(coord) eq 2 then begin
    alpha=coord[0]*!dtor
    delta=coord[1]*!dtor    
endif else begin
    alpha=coord[0,*,*]*!dtor
    delta=coord[1,*,*]*!dtor
endelse


alphaP=0.
deltaP=0.

thetaP=!pi/2
phiP=0.


phi=phiP+atan(-cos(delta)*sin(alpha-alphaP),sin(delta)*cos(deltaP)-cos(deltaP)*sin(deltaP)*cos(alpha-alphaP))
theta=asin(sin(delta)*sin(deltaP)+cos(delta)*cos(deltaP)*cos(alpha-alphaP))

if phi lt 0 then phi+=2*!pi


return,[phi,!pi/2-theta]*!radeg
end