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