FUNCTION EarthSky3DLoc, UT  , $
        cv2carrington=cv2carrington, $
        nra         = nra       , $
        nde         = nde       , $
        nrr         = nrr       , $
        drr         = drr       , $
        equator     = equator   , $
        degrees     = degrees   , $
        ra          = ra        , $
        dec         = dec       , $
        zero_point  = zero_point, $
        zero_phase  = zero_phase, $
        zero_shift  = zero_shift, $
        rr_earth    = rr_earth  , $
        ut_earth    = tt        , $
        pa_earth    = pa_earth  , $
        elo_earth   = elo_earth , $
        elo_sun     = elo_sun   , $
        center      = center    , $
        _extra      = _extra

    @compile_opt.pro                            ; On error, return to caller
InitVar, equator      , /key
InitVar, cv2carrington, /key

rpm = ToRadians(degrees=degrees)
pi  = !dpi/rpm
pi2 = pi/2

InitVar, zero_shift, 0.0

; Geocentric ecliptic location of Sun in ecliptic and equatorial coordinates

InitVar, center, jpl_body(/earth,/string)

SunEq = big_eph(UT      , $
    center  = jpl_body(/sun,/string), $
    body    = center    , $
    /to_equatorial      , $
    /to_sphere          , $
    /precess            , $
    degrees = degrees   , $
    /onebody            , $
    /silent             )

SunEq[0] = AngleRange(SunEq[0]+pi,degrees=degrees)
SunEq[1] = -SunEq[1]

;SunEq = big_eph(UT     , $
;   body    = jpl_body(/sun,/string), $
;   center  = center    , $
;   /to_equatorial      , $
;   /to_sphere          , $
;   /precess            , $
;   degrees = degrees   , $
;   /onebody            , $
;   /silent             )

SunEcl = CvSky(UT,from_equatorial=SunEq,/to_ecliptic,degrees=degrees,/silent)

InitVar, nra, 360/5     ; # points in ra direction
InitVar, nde, 180/5     ; # points in decl direction

InitVar, nrr, 20        ; # steps along line of sight
InitVar, drr, 0.1       ; Step size along line of sight (AU)

skygrid = nra*nde NE 0

zero_point = AngleRange(([SunEcl[0],SunEq[0]])[equator]+zero_shift,/degrees)
zero_phase = zero_point ; Longitude/ra Sun

CASE skygrid OF

0: BEGIN
    SyncArgs, ra, dec
    nlos = n_elements(ra)

    ;ra = zero_phase+ra                         ; Center on Sun
    rr = drr*gridgen(nrr, /open)                ; Distances along line of sight (AU)

    ; Set up collection of lines of sight by combining every ra with every declination
    ; and with every distance along line of sight.

    R = dblarr(3,nlos*nrr, /nozero)

    R[0,*] = SuperArray(ra ,nrr ,/trail)
    R[1,*] = SuperArray(dec,nrr ,/trail)
    R[2,*] = SuperArray(rr ,nlos,/lead )
END

1: BEGIN
    nra = long(nra)
    nde = long(nde)
    nrr = long(nrr)

    ; ra covers 360 deg in right ascension or ecliptic longitude relative to the
    ; Sun, going from 180 deg west of the Sun to 180 deg east of the Sun
    ; dec covers declinations or latitudes going from -90 to +90
    ; RR covers geocentric distances from 0 to nrr*drr

    ra = gridgen(nra, /open, range=[-pi,pi])    ; Right ascensions/ecliptic longitudes [-180,+180]

    ; ra are right ascensions or ecliptic longitudes

    ra += zero_phase                            ; Center on Sun
    dec = gridgen(nde, /open, range=[-pi2,pi2]) ; Latitudes/declinations
    rr  = drr*gridgen(nrr, /open)               ; Distances along line of sight (AU)

    ; Set up collection of lines of sight by combining every RA with every declination
    ; and with every distance along line of sight.

    R = dblarr(3,nra*nde*nrr, /nozero)

    R[0,*] = SuperArray(ra,nde*nrr,/trail)
    R[1,*] = SuperArray(SuperArray(dec,nra,/lead),nrr,/trail)
    R[2,*] = SuperArray(rr,nra*nde,/lead)
END

ENDCASE

; Convert ra and declinations to geocentric ecliptic coordinates

IF equator THEN R = CvSky(UT,from_equatorial=R,/to_ecliptic,degrees=degrees,/silent)

; Change from geocentric to heliocentric perspective
; !!! oldelo is angle Sun-Observer-Point on los (i.e. the elongation)
;   newelo is angle Observer-Sun-Point on los
; These angles are needed to compute the angle between the radial direction
;   and the direction perpendicular to the line of sight. This angle in turn
;   is needed in the integration for the IPS velocity.

R[0,*] -= SunEcl[0]                         ; Ecliptic longitude relative to Sun
R[2,*] /= SunEcl[2]                         ; Normalize distance to Sun-Earth distance

; Convert from geocentric to heliocentric perspective

R = CvPointOnLos(R,oldelo=elo_earth,oldphase=pa_earth,newelo=elo_sun,degrees=degrees)

R[2,*] *= SunEcl[2]                         ; Convert distance back to AU
R[0,*] += SunEcl[0]+pi                      ; Heliocentric ecliptic longitude
                                            ; Convert from heliocentric ecliptic to heliographic coordinates
R = CvSky(UT,from_ecliptic=R,_extra=_extra,degrees=degrees,/silent,euler_info=euler_info)

IF cv2carrington THEN BEGIN
    to_heliographic = strpos(euler_info,'to heliographic') NE -1
    IF NOT to_heliographic THEN message, 'invalid use of "/cv2carrington" keyword: no heliographic coordinates'
    R[0,*] = Carrington(UT, near_longitude=R[0,*], degrees=degrees, /get_variable)
ENDIF

CASE skygrid OF
0: BEGIN
    R = reform(R,3,nlos,nrr, /overwrite)
    elo_earth = reform(elo_earth, nlos, nrr, /overwrite)
    elo_sun   = reform(elo_sun  , nlos, nrr, /overwrite)
    pa_earth  = reform(pa_earth , nlos, nrr, /overwrite)
END
1: BEGIN
    R = reform(R,3,nra,nde,nrr, /overwrite)
    elo_earth = reform(elo_earth, nra,nde,nrr, /overwrite)
    elo_sun   = reform(elo_sun  , nra,nde,nrr, /overwrite)
    pa_earth  = reform(pa_earth , nra,nde,nrr, /overwrite)
END
ENDCASE

; R now contains all 3D locations (in heliographic coordinates) where values
; need to obtained by interpolation on the input 3D arrays.

; Location Earth in heliographic coordinates

rr_earth = [SunEcl[0]+pi,-SunEcl[1],SunEcl[2]]
rr_earth = CvSky(UT,from_ecliptic=rr_earth,_extra=_extra,degrees=degrees,/silent)

IF cv2carrington THEN   $
    rr_earth[0] = Carrington(UT, near_longitude=rr_earth[0], degrees=degrees, /get_variable)

TT = UT

RETURN, R  &  END