;Mapping from image space into heliographic coordinates
; mfw -- 99 mar 11
; (see manila folder "Transformation from  image coord's
; to heliographic coord's -- approximate treatment)

pro remap_lct, B0, gx, gy, nx, ny, rsun, mask, flow_map, $
 theta, phi, dtheta, dphi

;mask=1 <-> good points, mask=0 <-> bad points

;initialization
eps=1e-20
cb = cos(B0)
sb = sin(B0)
cx=(nx-1.)/2.
cy=(ny-1.)/2.
s=size(flow_map)
mx=s[2]
my=s[3]

;map from pixel coordinates to "relative" image coordinates
eta=mask*(gx-cx)/rsun
zeta=mask*(gy-cy)/rsun
xi=sqrt(1.-eta^2-zeta^2)

;map to heliographic cartesian coordinates
xip  = -sb*zeta + cb*xi
etap = eta
zetap = cb*zeta + sb*xi

;determine heliographic colatitude, longitude
; for points on the solar surface:
theta = acos(zetap)
ct=zetap
st=sin(theta)
phi=asin(etap/st)
;conventionalize longitude values
sxip = xip/(abs(xip)+eps) ;is there a cleaner way to compute
setap = etap/(abs(etap)+eps) ; the algebraic sign of a number?
phi=sxip*phi+(1.-sxip)*setap*(!pi)/2.

;determine theta and phi components of flow velocities from
; cartesian proper motion components (assuming that the velocities
; are tangent to the solar surface). 
cp=cos(phi)
sp=sin(phi)
a=replicate(0.0,2,2,mx,my)
dtheta=theta ;will become meridional flow velocity component, in radians per ..
dphi=phi ;will become azimuthal flow velocity component, in rad./ ..
a[0,0,*,*]=rsun*ct*sp
a[1,0,*,*]=rsun*st*cp
a[0,1,*,*]= -rsun*(cb*st+sb*ct*cp)
a[1,1,*,*]= rsun*sb*st*sp
for ix=0,mx-1 do begin ; would be nice idl could invert an array
 for iy=0,my-1 do begin ; of arrays
  b=invert(a[*,*,ix,iy]) 
  dtheta[ix,iy]=b[0,0]*flow_map[0,ix,iy]+b[1,0]*flow_map[1,ix,iy]
  dphi[ix,iy]=b[0,1]*flow_map[0,ix,iy]+b[1,1]*flow_map[1,ix,iy]
 endfor
endfor

end