pro helio_trans2,p,p0,l0,b0,l,b,r,x,y,z,rsun
;transforms heliographic coordinates [l,b,r] into cartesian coordinates [x,y,z] 
;	p              = position angle (N-->E) 
;	p0	       = rotation (roll) angle of image (N-->E)
;	l0,b0	       = heliografic longitude and latitude of disk center.
;	l[n],b[n]      = array of [l,b] coordinates
;spatial scales of R,X,Y,Z are specified by (photospheric) solar radius RSUN
;	r[n]	       = distance from Sun center (in same units as RSUN) 
;	x[n],y[n],z[n] = array of [x,y,z] coordinates (in same units as RSUN)
;	rsun	       = radius of sun (in arbitrary units)
;M.J.Aschwanden, March 1998, UMd

if (n_params(0) le 10) then rsun=1.	;length units in units of solar radii (old default)

;conversion deg-->rad
q	=(!pi/180.)
p0_	=p0*q
l0_	=l0*q
b0_	=b0*q
p_	=p *q
l_	=l *q
b_	=b *q
;origin of coordinate system (x,y,z)
x1	=0.
y1	=0.
z1	=r/rsun
;rotation by latitude angle 
x2	=x1
y2	=y1*cos(-b_)-z1*sin(-b_)
z2	=y1*sin(-b_)+z1*cos(-b_)
;rotation by longitude angle difference 
x3	=x2*cos(l0_-l_)-z2*sin(l0_-l_)
y3	=y2
z3	=x2*sin(l0_-l_)+z2*cos(l0_-l_)
;rotation by latitude angle 
x4	=x3
y4	=y3*cos(b0_)-z3*sin(b0_)
z4	=y3*sin(b0_)+z3*cos(b0_)
;rotation by position angle 
x 	=x4*cos(p_+p0_)-y4*sin(p_+p0_)
y 	=x4*sin(p_+p0_)+y4*cos(p_+p0_)
z 	=z4
;length unit
x	=x*rsun
y	=y*rsun
z	=z*rsun
end