; (Lat,Lon) to SSM/I (X,Y) (northern hemisphere only).
; Input:  alat, alon (vectors)
; Output: x, y (vectors)

PRO LLtoXY, alat, alon, x, y

rot = 45.0 ; orientation of +X axis, degrees

nel = n_elements(alat)
x = fltarr(nel)
y = fltarr(nel)
t0 = fltarr(nel)
rho = fltarr(nel)

; Definition of constants

re=6378.273		; Radius of earth (km)  -Hughes Ellipsoid
e2=DOUBLE(0.006693883)	; Eccentricity of earth -Hughes Ellipsoid
e=sqrt(e2)
slat=70.		; Standard parallel

XP = where(alat LT 89.995, num)
if num GT 0 then begin ; for points south of 89.995 N latitude

   t0(XP)=TAN((!DPI/4.)-(alat(XP)/(2.*!RADEG)))/     $
          ((1.-e*SIN(!DTOR*alat(XP)))/               $
           (1.+e*SIN(!DTOR*alat(XP))))^(e/2)

   t1=    TAN((!DPI/4.)-(slat/(2.*!RADEG)))/         $
          ((1.-e*SIN(!DTOR*slat))/                   $
           (1.+e*SIN(!DTOR*slat)))^(e/2)

   cm=COS(!DTOR*slat)/SQRT(1.-e2*(SIN(!DTOR*slat)^2))
   rho(XP)=re*cm*t0(XP)/t1

   x(XP) =  rho(XP)*SIN(!DTOR*(alon(XP)+rot))
   y(XP) = -rho(XP)*COS(!DTOR*(alon(XP)+rot))

endif

NP = where(alat GE 89.995, num)
if num GT 0 then begin ; for points north of 89.995 N latitude
   x(NP) = replicate(0.0, num)
   y(NP) = replicate(0.0, num)
endif

END

;================================================================

; SSM/I (X,Y) to (Lat,Lon) (northern hemisphere only).
; Input:  x, y (vectors)
; Output: alat, alon (vectors)

PRO XYtoLL, x, y, alat, alon

rot = 45.0 ; orientation of +X axis, degrees

; Definition of constants

re=6378.273		; Radius of earth (km)  -Hughes Ellipsoid
e2=DOUBLE(0.006693883)	; Eccentricity of earth -Hughes Ellipsoid
e=sqrt(e2)
slat=70.		; Standard parallel

rho = sqrt(x*x + y*y)
cm  = COS(!DTOR*slat)/SQRT(1.-e2*(SIN(!DTOR*slat)^2))
t   = TAN((!DPI/4.)-(slat/(2.*!RADEG)))/         $
      ((1.-e*SIN(!DTOR*slat))/                   $
       (1.+e*SIN(!DTOR*slat)))^(e/2)
tt  = rho*t/(re*cm)
chi = (!DPI/2.0)-2.0*ATAN(tt)

alat = chi + ((e2/2)+(5*(e2^2)/24)+((e2^3)/12)) * SIN(2*chi) $
              + ((7*(e2^2)/48)+(29*(e2^3)/240)) * SIN(4*chi) $
                               + (7*(e2^3)/120) * SIN(6*chi)
alat = alat * !RADEG
alon = ATAN(x,-y)*!RADEG - rot

neg = where(alon LE -180.0, num)
if num GT 0 then alon(neg) = alon(neg) + 360.

END