;+ ; NAME: ; smei_zld_model ; PURPOSE: ; Implements the current analytic fit to the zodiacal light ; distribution for the SMEI observations ; CATEGORY: ; camera/idl/dust ; CALLING SEQUENCE: FUNCTION smei_zld_model, tt, rr, $ ndim = ndim , $ range = range_ , $ degrees = degrees , $ ecliptic = ecliptic , $ equatorial = equatorial, $ model_1 = model_1 , $ model_2 = model_2 , $ model_3 = model_3 , $ model_4 = model_4 , $ noslab = noslab , $ version = version , $ get_version = get_version,$ silent = silent , $ camera = camera , $ keep_gegenschein = keep_gegenschein ; INPUTS: ; tt array[1]; type: time structure ; UT time at which zodiacal light brightness is required ; rr array[2,*]; type: float or double ; angles defining the lines of sight for which brightness is required ; ; rr[0,*] is ecliptic longitude of los relative to Sun ; rr[1,*] is ecliptic latitude of los ; (i.e. a sun-centered ecliptic coordinate system is assumed) ; ; if /ecliptic is SET: ; rr[0,*] is ecliptic longitude of los ; rr[1,*] is ecliptic latitude of los ; ; if /equatorial is SET: ; rr[0,*] is right ascension of los ; rr[1,*] is declination of los ; ; The units of the angles are set by keyword /degrees ; NOTE: it is allowed to specify rr as an array[3,*] with the ; 3rd element in the first dimension the line of sight distance ; (only the angles are used) ; OPTIONAL INPUT PARAMETERS: ; ndim=ndim array[2]; type: integer; default: none ; if present, a full-sky map is set up in the coordinate system ; implied by setting of /ecliptic and /equatorial ; The lines of sight used are returned in rr ; range=range array[2,2]; type: float; default: none ; range of grid in longitude and latitude. ; Used in conjunction with keyword 'ndim'. ; range[0:1,0]; start and end longitude ; range[0:1,0]; start and end latitude ; These are the outer bin edges, NOT the bin centers. ; /degrees if SET the angles in argument 'rr' and 'range' ; are in degrees; default is radians ; /equatorial if SET then 'rr' is specified in equatorial coordinates ; /ecliptic if SET then 'rr' is specified in ecliptic coordinates ; if neither /equatorial nor /ecliptic is set then 'rr' is ; assumed to be in sun-centered ecliptic coordinates ; version=version scalar; type: float ; version number of zodiacal light model ; OUTPUTS: ; result array[*]; type: double ; the model brightness of the zodiacal dust cloud ; rr array[2,n,m]; type: double ; (only if keyword ndim is used) ; The lines of sight used to create a full-sky map ; ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, big_eph, jpl_body, CvSky, AngleUnits, AngleRange, gridgen ; EXAMPLE: ; To get the zodiacal dust cloud brightness on an equatorial grid with a resolution ; of one degree at the current time: ; ; tt = TimeSystem() ; rr = gridgen([360,180],/open,range=[[0,360],[-90,90]],/multid) ; zz = smei_zld_model(tt,rr,/degrees,/equatorial) ; ; returns an double precision array zz[360,180]. ; PROCEDURE: ; Purely empirical fit ; ; The brightness scale corresponds to camera 2 at nominal gain 1. ; To subtract this model from a skymap the map needs to multiplied by the ; time- and camera-dependent gain factor returned by href=smei_camera_gain=. ; MODIFICATION HISTORY: ; AUG-2006, Jordan Vaughan (jtvaugha@ucsd.edu), Paul Hick (UCSD/CASS) ; Based on Andys zodiac_model.for, but substantially rewritten. ; DEC-2007, Paul Hick (UCSD/CASS) ; Added second zodiacal light model (keyword /model_2), ; based on the model first used for the comet movies. ; The old model is available with keyword /model_1. ; APR-2008, Paul Hick (UCSD/CASS) ; Next iteration on zld model added as /model_3. ; (/model_2 remains the default). ; Added keyword "version" ; NOV-2008, Paul Hick (UCSD/CASS) ; Next iteration on zld model added as /model_4. ; DEC-2008, Paul Hick (UCSD/CASS) ; Model 4 is now the default model (was model 2) ; Bug fix in new model (V4.01) ;- InitVar, ecliptic , /key InitVar, equatorial , /key InitVar, model_1 , /key InitVar, model_2 , /key InitVar, model_3 , /key InitVar, model_4 , /key InitVar, noslab , /key InitVar, silent , 0 InitVar, get_version, /key InitVar, keep_gegenschein, /key model_4 OR= (1-model_1) AND (1-model_2) AND (1-model_3) CASE 1 OF model_1: version = 1.0 model_2: version = 2.0 model_3: version = 3.0 model_4: version = 4.01 ENDCASE IF get_version THEN RETURN, version ; Internally all angle calculations are done in radians pi = !dpi halfpi = pi/2.0d0 dpr = 180.0d0/pi ; Degrees per radian ; Geocentric position Sun at time tt in ecliptic and equatorial coordinates rsun_equ = big_eph( tt , $ body = jpl_body(/sun ,/string) , $ center = jpl_body(/earth,/string) , $ /to_equatorial , $ /to_sphere , $ /precess , $ /onebody , $ /silent ) rsun_ecl = CvSky(tt, from_equatorial=rsun_equ, /to_ecliptic, /silent) CASE IsType(ndim,/defined) OF 0: BEGIN ; Convert to 2D array with first dimension 2 or 3 ; Convert to double precision (this mainly deals with integer input angles) ; Drop the distance if present rrsz = size(rr,/dim) rzod = reform(double(rr),rrsz[0],n_elements(rr)/rrsz[0]) IF rrsz[0] EQ 3 THEN rzod = rzod[0:1,*] ; Convert input angles to radians ; From here on all angles are in radians. ; Where degrees are needed multiply by dpr rzod = AngleUnits(from_angle=rzod, degrees=degrees, /to_radians) END 1: BEGIN CASE IsType(range_, /defined) OF 0: BEGIN ; Full-sky grid in radians CASE 1 OF ecliptic : range = pi*[[ 0.0d0,2.0d0],[-0.5d0,0.5d0]] equatorial: range = pi*[[ 0.0d0,2.0d0],[-0.5d0,0.5d0]] ELSE : range = pi*[[-1.0d0,1.0d0],[-0.5d0,0.5d0]] ENDCASE END 1: range = AngleUnits(from_angle=range_, degrees=degrees, /to_radians) ENDCASE rzod = gridgen(ndim,range=range,/open,/double) rr = reform(rzod,[2,ndim]) InitVar, degrees, /key IF degrees THEN rr *= dpr END ENDCASE ; Convert equatorial input grid to ecliptic coordinates IF equatorial THEN rzod = CvSky(tt, from_equatorial=rzod, /to_ecliptic,/silent) ; Convert ecliptic longitudes to Sun-centered ecliptic longitude IF ecliptic OR equatorial THEN rzod[0,*] -= rsun_ecl[0] ; Solar elongations elo = acos(cos(rzod[0,*])*cos(rzod[1,*])) elo = reform(elo,/overwrite); Loose leading dummy dimension ; Position angle from ecliptic north (counterclockwise) ; Range is [-pi,pi] pa = atan(sin(rzod[0,*]), tan(rzod[1,*])) pa = reform(pa,/overwrite) ; Loose leading dummy dimension power_r = 2.3d0 CASE 1 OF model_1: BEGIN ; Constants used in the parameterization. ; Uses Bernie's original parameterization to generate a ; zodiacal sky brightness. ; This is then subtracted from the selected skymap. A0 = 22.4000d0 ; Constant background A = 157453.6090d0 ; Used in (1/elo)^2 part A1 = 87.7441d0 A2 = -56.6437d0 B = 248322.7810d0 ; Used in (1/elo)^3 part B1 = -2.6228d0 B2 = 1.3625d0 C = 55.0000d0 ; Used for Gegenschein C1 = 25.0000d0 R = 1.5000d0/dpr ; Inclination of zodiacal dust plane to ecliptic ?? Omega = 96.0000d0/dpr ; ?? Ascending node of zodiacal dust plane on ecliptic ?? BrightBand = 13.0d0 ; Used to fit the dust bands OmegaBand = 180.0d0/dpr ; Change 2-16-06 values to reduce zodiac near Sun and eliminate Y. A = 280000.0d0 ; 5-24-06 B = 550000.0d0 ; 5-24-06 Z1 = 0.000230d0 ; 5-25-06 ; Some sort of correction to the position angle ; Must have something to do with the tilt of the zld plane to the ecliptic ; in which case I would expect pa = pa+R*cos((run_ecl[0]-Omega)) pa += R*cos((rsun_equ[0]-Omega)) ; Angular distance to plane of symmetry of zld cloud (-90<=dec<=90) dec = halfpi-abs(pa) ; 90 at north pole, 0 in plane; -90 at south pole tmp = abs(dec) ; Need abs(dec) again for fudge 1, so save it in tmp sin_dec = sin(tmp) ; 0 in plane; 1 at both poles pole sin2_dec = sin_dec*sin_dec ; 0 in plane; 1 at both poles ; ====== FUDGE 1: "fix" cusp at abs(pa) = 90 ; ; sin_dec has a cusp where the plane of symmetry is crossed (abs(pa)=90;sin_dec=0) ; The 2 degrees around pa +/-90 degrees substitutes a quadratic in place of linear ; Values and slope are matched at +/-2 degrees. ; This looked odd within a few degrees of the Sun, when 5 degrees ; was used in place of 2 here. This may be still true. ; Updated on 5-26-2006 ; delta = 0.034899497d0 ; sin(2 deg) delta = 2.0/dpr i = where(tmp LE delta) ; tmp = abs(dec) IF i[0] NE -1 THEN BEGIN tmp = tmp[i]/delta sin_dec[i] = 0.01745d0*(1.0d0+tmp*tmp) ENDIF ; ====== END FUDGE 1 ; The main zodiac model: ; - elongation dependence has (1/elo)^2 and (1/elo)^3 terms ; - position angle dependence has abs(sin(dec)) and sin^2(dec) terms ; ====== FUDGE 2: force 1/elo to zero for elo -> 180.0 deg ; 1/elo with extra term for elo > 90 to get zero at elo=180 tmp = 1.0d0/(elo*dpr) i = where(elo GT halfpi) ; In antisolar hemisphere IF i[0] NE -1 THEN tmp[i] *= sqrt(sin(elo[i])) ; ====== END FUDGE 2 zld = ( A*(1.0d0 + A1*sin_dec + A2*sin2_dec)*tmp + B*(1.0d0 + B1*sin_dec + B2*sin2_dec) )*tmp*tmp ; ====== FUDGE 3: Goose it up near the Sun and ecliptic tmp = elo*dpr-70.0d0 i = where(tmp LT 0.0d0) ; Less then 70 deg elongation IF i[0] NE -1 THEN BEGIN tmp = tmp[i] zld[i] *= 1.0d0+Z1*(1.0d0-sin2_dec[i])^6*tmp*tmp ENDIF ; ====== END FUDGE 3 ; Constant background, indistinguishable from Sidereal zld += A0 ; ===== Add Exponential Gegenschein ; Remove the cusp by using a hyperbola in antisolar elongation. ; NOTE: With characteristic angle at 1 deg, this thing practically DOES have a cusp. tmp = (pi-elo)*dpr zld += C*exp((1.0d0-sqrt(1.0d0+tmp*tmp))/C1) ; ====== FUDGE 4: try to model dust bands ; The full range of possible values for the ecliptic longitude of the Sun relative ; to OmegaBand is divided in two 180 degree sections. ; The correction is BrightBand*abs(cos(BandLng))*sin(elo) ; multiplied by a factor between 0 and 1: ; ; abs(BandLng) > 90: ; factor = 0 for dec < -10 ; factor = dec/20+0.5 for -10 < dec < +10, i.e. increases linearly from 0 to 1 ; factor = 1 for dec > +10 ; abs(BandLng) < 90: ; factor = 1 for dec < -10 ; factor = -dec/20+0.5 for -10 < dec < +10, i.e. decreases linearly from 1 to 0 ; factor = 0 for dec > +10 BandLng = AngleRange(rsun_ecl[0]-OmegaBand,/pi) ; -pi < BandLng <= pi zld += BrightBand*abs(cos(BandLng))*sin(elo)* $ ( ( (1-2*(abs(BandLng) LT halfpi))*(dec/(20.0d0/dpr))+0.5d0 ) > 0.0d0 < 1.0d0 ) ; ====== END FUDGE 4 ; Everything is adjusted for Sun-Earth distance using 2.3 power from Astrophysical Quantities. ; Sun-Earth distance correction: zld map is reduced for > 1 AU and vice versa. zld /= rsun_ecl[2]^power_r END model_2: BEGIN ; This uses the first of Steve Price's two suggested functional forms. ; See Lamy and Perrin (Astr. Ap. 163, 269 (1986)), with a VSF from Ney ; and Merrill (Science 194, 1051 (1976)), which uses a sum of ; cosines up through power 3, all divided by sin(elo)^2.3. cos_elo = cos(elo) sin_elo = sin(elo) sin_elo_power_r = sin_elo^power_r ; Change angular scale zld = 75.0d0+120.0d0*cos_elo sunward_hemisphere = cos_elo GT 0.0d0 antisolar_hemisphere = 1B-sunward_hemisphere sunward = where(sunward_hemisphere, complement=antisolar) tmp = cos_elo*cos_elo IF sunward[0] NE -1 THEN BEGIN zld[sunward] += tmp[sunward]*(-185.0d0+140.0d0*cos_elo[sunward]) zld[sunward] /= sin_elo_power_r[sunward] ENDIF IF antisolar[0] NE -1 THEN BEGIN ; We include an ad hoc 6th power cosine and exponential in the antisolar ; hemisphere to remove the Gegenschein residue. tmp = tmp[antisolar] zld[antisolar] += tmp*(154.0d0+88.0d0*cos_elo[antisolar]+20.0d0*tmp*tmp) zld[antisolar] += 20.0d0*exp(-(pi-elo[antisolar])/(6.0d0/dpr)) ENDIF ; Ecliptic latitude and longitude relative to Sun hlat = reform(rzod[1,*]) abs_hlng = abs(reform(rzod[0,*])) ; Dropoff exponential in sin(rLat), with elo slope ; The cusp is removed by using a hyperbola of rLat ; in the exponential ; This has a flat bottom, pulls the rest of the curve down ; by 1.5 degrees, with the transition at about 1.5 degrees ; This may generate another cusp at +-90 deg., hopefully ; too small to see ... {to fix, multiply by 90/88.5} ; The elongation offset below comes from observed slopes on ; the hlat - elo plot of log10 (Observed/Model) ; An offset range between about 30 to 45 is permissable. hyperbole_width = 1.5d0/dpr ; rzod[1,*] is ecliptic latitude ; Hyperbola: tmp[0] = 0 tmp = hyperbole_width*(sqrt(1.0d0+(hlat/hyperbole_width)^2)-1.0d0) zld /= 10.0d0^( abs(sin(tmp))/(0.007d0*(elo*dpr+40.0d0)) ) ; In the sunward hemisphere, include term propto cos(tmp) zld += sunward_hemisphere*( 30.0d0*(1.0d0/sin_elo_power_r-1.0d0)*cos(tmp) ) IF NOT noslab THEN BEGIN Omega = 70.0d0/dpr ; ?? Ascending node of zodiacal dust plane on ecliptic ?? R1 = 5.0d0/dpr REW = 6.0d0 RRR = 6.0d0 ; Slab brightness?? tmp = rsun_ecl[0]-Omega NS = sin(tmp) ; Max when sun is 90 deg away from node EW = cos(tmp) ; Max when sun is in node ; Blur edge between ecl latitudes -R1 and +R1 ; NS > 0 if Earth is above plane of dust; in antisolar ; hemisphere the dust brightness is biased to the south ; NS > 0: zslab = 1 if hlat < -R1 (southern hemisphere) ; = 0 if hlat > +R1 (northern hemisphere) ; with linear transition between -R1 and +R1 ; NS < 0 if Earth is below plane of dust; in antisolar ; hemisphere the dust brightness is biased to the north ; NS < 0: zslab = 0 if hlat < -R1 (southern hemisphere) ; = 1 if hlat > +R1 (northern hemisphere) ; with linear transition between -R1 and +R1 zslab = ( (1-2*(NS GT 0.0d0))*(hlat/(2*R1)) +0.5d0 ) > 0.0d0 < 1.0d0 ; Close to the Sun the NS asymmetry is less tmp = sunward_hemisphere*sin_elo + antisolar_hemisphere*1.0d0 zslab = RRR*abs(NS)*(1.0d0+tmp*(zslab-1.0d0)) ; Attempt at the EW/NS component, greatest when Earth passes ; through the nodes of the zodiacal dust plane on the ecliptic ; Bend it by R1 degree tilt ; 0 for elo=0 and 180; max for elo=90 ; For abs(EW) close to one (sun in node) the maximum is where ; the dust is farthest from the ecliptic (elo close to 90). R1EW = R1*abs(EW)*sin_elo ; >=0 ; EW > 0 sun closer to ascending node; ; dust brightness is biased to the north/south in the ; east/west hemisphere (pa > 0/pa < 0) ; EW < 0 sun closer to descending node; ; dust brightness is biased to the south/north in the ; east/west hemisphere (pa > 0/pa < 0) ; blurNS is 1.0 on the ecliptic (hlat=0) ; EW > 0: ; pa > 0 (east): blurNS = 1 if hlat > 0 (north) ; blurNS = 0 if hlat < -R1EW (south) ; pa < 0 (west): blurNS = 0 if hlat > +R1EW (north) ; blurNS = 1 if hlat < 0 (south) ; EW < 0: ; pa > 0 (east): blurNS = 0 if hlat > +R1EW (north) ; blurNS = 1 if hlat < 0 (south) ; pa < 0 (west): blurNS = 1 if hlat > 0 (north) ; blurNS = 0 if hlat < -R1EW (south) ; So this correction adds to two opposite quadrants of sky blurNS = ( (1-2*(EW GT 0.0d0))*(1-2*(pa GT 0.0d0))*(hlat/R1EW) +1.0d0 ) > 0.0d0 < 1.0d0 ; Andy's "feathering" of the NS transition at hlng = 0 (through Sun) ; and hlng=180 (through anti-solar direction). ; blurNS LT 1.0d0 identifies the two opposite quadrants with no slab ; It adds an exponential dropoff from 1 (at the NS direction) to zero ; moving into the quadrants with no slab. ;blurEW = (blurNS LT 1.0d0)*(exp(-abs_hlng/(R1/3.0d0))+exp(-(pi-abs_hlng)/(R1/2.0d0))) blurEW = 0.0d0 zslab += REW*abs(EW)*abs(cos_elo)*(blurNS+blurEW) ; Increase slab brightness towards Sun. ; Drop off the slab brightness in antisolar hemisphere zslab *= sunward_hemisphere/sin_elo_power_r + antisolar_hemisphere*sin_elo zld += zslab ENDIF ; If elo=90 and hlat=90 (ecliptic poles), this +DC offset below should ; agree with AQ p329 (22.4), or p330 (24) ; Nov 2 value is 6 {from the 75 constant above top line for zod, ; times 10^(-1/(.007*130))}; but see just below zld += 16.0d0 ; DC offset, a case can be made that this should be "zld+=REW" ; Diminish zodical light a bit near ecliptic and around 90 degrees, ; This probably needs refinement, especially towards the Sun. zld -= 8.0d0*sin_elo*exp(-(hlat*dpr)^2/10.0d0 ) zld /= rsun_ecl[2]^power_r END model_3: BEGIN ; Andys revision of Model II (2008/04/22) ; (see $SMEI/user/abuffington/zld/zodiac_modelb_3.f) ; This uses the first of Steve Price's two suggested functional forms. ; See Lamy and Perrin (Astr. Ap. 163, 269 (1986)), with a VSF from Ney ; and Merrill (Science 194, 1051 (1976)), which uses a sum of ; cosines up through power 3, all divided by sin(elo)^2.3. cos_elo = cos(elo) sin_elo = sin(elo) sin_elo_power_r = sin_elo^power_r ; Change angular scale ;zld = 60.0d0+120.0d0*cos_elo ; @@@@@@@@@@@@ zld = 65.0d0+120.0d0*cos_elo ; @@@@@@@@@@@@ sunward_hemisphere = cos_elo GT 0.0d0 antisolar_hemisphere = 1B-sunward_hemisphere sunward = where(sunward_hemisphere, complement=antisolar) tmp = cos_elo*cos_elo IF sunward[0] NE -1 THEN BEGIN zld[sunward] += tmp[sunward]*(-185.0d0+140.0d0*cos_elo[sunward]) zld[sunward] /= sin_elo_power_r[sunward] ENDIF IF antisolar[0] NE -1 THEN BEGIN ; We include an ad hoc 6th power cosine and exponential in the antisolar ; hemisphere to remove the Gegenschein residue. tmp = tmp[antisolar] ;zld[antisolar] += tmp*(154.0d0+88.0d0*cos_elo[antisolar]+22.0d0*tmp*tmp) ;@@@@@@ zld[antisolar] += tmp*(154.0d0+88.0d0*cos_elo[antisolar]) ;@@@@@@ ENDIF ; Ecliptic latitude and longitude relative to Sun hlat = reform(rzod[1,*]) abs_hlng = abs(reform(rzod[0,*])) ; Dropoff exponential in sin(rLat), with elo slope ; The cusp is removed by using a hyperbola of rLat ; in the exponential ; This has a flat bottom, pulls the rest of the curve down ; by 1.5 degrees, with the transition at about 1.5 degrees ; This may generate another cusp at +-90 deg., hopefully ; too small to see ... {to fix, multiply by 90/88.5} ; The elongation offset below comes from observed slopes on ; the hlat - elo plot of log10 (Observed/Model) ; An offset range between about 30 to 45 is permissable. hyperbole_width = 1.5d0/dpr ; rzod[1,*] is ecliptic latitude ; Hyperbola: tmp[0] = 0 tmp = hyperbole_width*(sqrt(1.0d0+(hlat/hyperbole_width)^2)-1.0d0) ;zld /= 10.0d0^( abs(sin(tmp))/(0.007d0*(elo*dpr+40.0d0)) ) ; @@@@@@@@@@@ zld /= 10.0d0^( abs(sin(tmp))/(0.009d0*(elo*dpr+40.0d0)) ) ; @@@@@@@@@@@ zld -= 8.0d0*(cos(tmp)-1.0d0) ; @@@@@@@@@ ; Small-scale Gegenschein term outside of the Latitude dropoff ; since it is "round in the sky" ;zld += 21.0d0*exp(-(pi-elo)/(6.0d0/dpr)) ; @@@@@@@@@@ ; Partial gaussian contribution sigma = 14 deg. ; (18 April 2008, for camera 1 only) ;zld += 5.0d0*exp(-(hlat*dpr)^2/400.0d0) ; @@@@@@@@@ zld += 6.0d0*exp(-(hlat*dpr)^2/512.0d0) ; @@@@@@@@@ ; In the sunward hemisphere, include term propto cos(tmp) zld += sunward_hemisphere*( 30.0d0*(1.0d0/sin_elo_power_r-1.0d0)*cos(tmp) ) IF NOT noslab THEN BEGIN Omega = 70.0d0/dpr ; Ascending node of zodiacal dust plane on ecliptic ?? R1 = 5.0d0/dpr REW = 6.0d0 RRR = 6.0d0 ; Slab brightness?? tmp = rsun_ecl[0]-Omega NS = sin(tmp) ; Max when sun is 90 deg away from node EW = cos(tmp) ; Max when sun is in node ; Blur edge between ecl latitudes -R1 and +R1 ; NS > 0 if Earth is above plane of dust; in antisolar ; hemisphere the dust brightness is biased to the south ; NS > 0: zslab = 1 if hlat < -R1 (southern hemisphere) ; = 0 if hlat > +R1 (northern hemisphere) ; with linear transition between -R1 and +R1 ; NS < 0 if Earth is below plane of dust; in antisolar ; hemisphere the dust brightness is biased to the north ; NS < 0: zslab = 0 if hlat < -R1 (southern hemisphere) ; = 1 if hlat > +R1 (northern hemisphere) ; with linear transition between -R1 and +R1 zslab = ( (1-2*(NS GT 0.0d0))*(hlat/(2*R1)) +0.5d0 ) > 0.0d0 < 1.0d0 ; Close to the Sun the NS asymmetry is less tmp = sunward_hemisphere*sin_elo + antisolar_hemisphere*1.0d0 zslab = RRR*abs(NS)*(1.0d0+tmp*(zslab-1.0d0)) ; Attempt at the EW/NS component, greatest when Earth passes ; through the nodes of the zodiacal dust plane on the ecliptic ; Bend it by R1 degree tilt ; 0 for elo=0 and 180; max for elo=90 ; For abs(EW) close to one (sun in node) the maximum is where ; the dust is farthest from the ecliptic (elo close to 90). R1EW = R1*abs(EW)*sin_elo ; >=0 ; EW > 0 sun closer to ascending node; ; dust brightness is biased to the north/south in the ; east/west hemisphere (pa > 0/pa < 0) ; EW < 0 sun closer to descending node; ; dust brightness is biased to the south/north in the ; east/west hemisphere (pa > 0/pa < 0) ; blurNS is 1.0 on the ecliptic (hlat=0) ; EW > 0: ; pa > 0 (east): blurNS = 1 if hlat > 0 (north) ; blurNS = 0 if hlat < -R1EW (south) ; pa < 0 (west): blurNS = 0 if hlat > +R1EW (north) ; blurNS = 1 if hlat < 0 (south) ; EW < 0: ; pa > 0 (east): blurNS = 0 if hlat > +R1EW (north) ; blurNS = 1 if hlat < 0 (south) ; pa < 0 (west): blurNS = 1 if hlat > 0 (north) ; blurNS = 0 if hlat < -R1EW (south) ; So this correction adds to two opposite quadrants of sky blurNS = ( (1-2*(EW GT 0.0d0))*(1-2*(pa GT 0.0d0))*(hlat/R1EW) +1.0d0 ) > 0.0d0 < 1.0d0 ; Andy's "feathering" of the NS transition at hlng = 0 (through Sun) ; and hlng=180 (through anti-solar direction). ; blurNS LT 1.0d0 identifies the two opposite quadrants with no slab ; It adds an exponential dropoff from 1 (at the NS direction) to zero ; moving into the quadrants with no slab. blurEW = (blurNS LT 1.0d0)*(exp(-abs_hlng/(R1/3.0d0))+exp(-(pi-abs_hlng)/(R1/2.0d0))) ;blurEW = 0.0d0 zslab += REW*abs(EW)*abs(cos_elo)*(blurNS+blurEW) ; Increase slab brightness towards Sun. ; Drop off the slab brightness in antisolar hemisphere zslab *= sunward_hemisphere/sin_elo_power_r + antisolar_hemisphere*sin_elo zld += zslab ENDIF ; If elo=90 and hlat=90 (ecliptic poles), this +DC offset below should ; agree with AQ p329 (22.4), or p330 (24) ; Nov 2 value is 6 {from the 75 constant above top line for zod, ; times 10^(-1/(.007*130))}; but see just below ;zld += 19.0d0 ; DC offset, a case can be made that this should be "zld+=REW" ; @@@@ zld += 12.0d0 ; DC offset, a case can be made that this should be "zld+=REW" ; @@@@ zld /= rsun_ecl[2]^power_r END model_4: BEGIN ; This uses the first of Steve Price's two suggested functional forms. ; See Lamy and Perrin (Astr. Ap. 163, 269 (1986)), with a VSF from Ney ; and Merrill (Science 194, 1051 (1976)), which uses a sum of ; cosines up through power 3, all divided by sin(elo)^2.3. cos_elo = cos(elo) sin_elo = sin(elo) sin_elo_power_r = sin_elo^power_r ; Change angular scale zld = 65.0d0+120.0d0*cos_elo sunward_hemisphere = cos_elo GT 0.0d0 antisolar_hemisphere = 1B-sunward_hemisphere sunward = where(sunward_hemisphere, complement=antisolar) tmp = cos_elo*cos_elo IF sunward[0] NE -1 THEN BEGIN zld[sunward] += tmp[sunward]*(-185.0d0+65.0d0*cos_elo[sunward]) zld[sunward] /= sin_elo_power_r[sunward] ENDIF IF antisolar[0] NE -1 THEN BEGIN zld[antisolar] += tmp[antisolar]*(154.0d0+88.0d0*cos_elo[antisolar]) ENDIF ; Ecliptic latitude and longitude relative to Sun hlat = reform(rzod[1,*]) abs_hlng = abs(reform(rzod[0,*])) ; Dropoff exponential in sin(rLat), with elo slope ; The cusp is removed by using a hyperbola of rLat ; in the exponential ; This has a flat bottom, pulls the rest of the curve down ; by 1.5 degrees, with the transition at about 1.5 degrees ; This may generate another cusp at +-90 deg., hopefully ; too small to see ... {to fix, multiply by 90/88.5} ; The elongation offset below comes from observed slopes on ; the hlat - elo plot of log10 (Observed/Model) ; An offset range between about 30 to 45 is permissable. width = 1.5d0/dpr ; rzod[1,*] is ecliptic latitude ; Hyperbola: tmp[0] = 0 tmp = hlat/width tmp = width*(sqrt(1.0d0+tmp*tmp)-1.0d0) zld /= 10.0d0^( abs(sin(tmp))/(0.009d0*dpr*(elo+40.0d0/dpr)) ) zld += 8.0d0*(1.0d0-cos(tmp)) ; In the sunward hemisphere, include term propto cos(tmp) zld += sunward_hemisphere*( 30.0d0*(1.0d0/sin_elo_power_r-1.0d0)*cos(tmp) ) ; This gaussian contribution is used with camera 2 when making ; C2{eq,np,sp} 3 grid files to be used in allskymap.for. It is ; empirically found to clear out a broad ecliptic band in the ; eq sidereal map, so the sidereal residue has "no memory" of ; the ecliptic. Here, "no memory" means down below a few ADUs, ; as best we can see... ; This arbitrarily defines the "sidereal component" to be subtracted ; from all maps, and "if it looks good, it's good". ; It doesn't strictly constrain the zodiacal-model parameters, ; although we should probably minimize the subsequent changes. ; Gaussian contribution sigma = 16 deg., 19 May 2008, works ; best for camera 1 tmp = hlat/(sqrt(512.0d0)/dpr) zld += 6.0d0*exp(-tmp*tmp) ; Subtract the double-exponential Gegenschein terms, as in our ; Icarus manuscript ; The weight formula here is an empirical match to the observed ; Gegenschein distribution near the antisolar direction. ; It uses 2% of the sine^2 of twice the angle times the antisolar ; elongation (deg.) to diminish the geometric weighting of the ; x,y components. ; This is ugly, but reproduces the lozenge shape that the ; Gegenschein has. ; This formula gets somewhat crazy for aEloP > 90 degrees, so we ; reduce Gegenschein to < 0.1 ADU at aEloP > 90 deg., using a ; Gaussian rolloff tmp = pi-elo ; Elongation from anti-solar direction exp_lat = 5.5d0*exp(-tmp/(4.0d0/dpr))+38.5d0*exp(-tmp/(19.5d0/dpr)) ; Latitude exp_lng = 5.5d0*exp(-tmp/(4.0d0/dpr))+38.5d0*exp(-tmp/(27.0d0/dpr)) ; Longitude arg_lng = (pi-abs_hlng)/tmp arg_lat = hlat/tmp arg_lng = arg_lng*arg_lng arg_lat = arg_lat*arg_lat exp_elo = arg_lng*exp_lng+arg_lat*exp_lat weight = (1.0d0-0.02d0*arg_lng*arg_lat*tmp*dpr) > 0.0d0 ; Roll it down fast for elo <= 120 degrees i = where(elo LE 120.0d0/dpr) IF i[0] NE -1 THEN BEGIN tmp = (120.0d0/dpr-elo[i])/(sqrt(300.0d0)/dpr) weight[i] *= exp(-tmp*tmp) ENDIF CASE keep_gegenschein OF 0: zld += weight*exp_elo 1: IF silent LE 1 THEN message, /info, 'retaining Gegenschein' ENDCASE ; Double exponential, with a 3 degree hyperbola ; derived from the deviation of PA from 90 degrees. tmp = (abs(pa)-halfpi)/(3.0d0/dpr) zld += (8800.0d0*exp(-(sqrt(1.0d0+tmp*tmp)-1.0d0)/10.0d0)-1200.0d0)* $ exp(-elo/(10.0d0/dpr)) IF NOT noslab THEN BEGIN ; NS is max when Sun is 90 deg away from node ; 78.25 is longitude ascending node of dust plane on ecliptic NS = sin(rsun_ecl[0]-78.25d0/dpr) above_plane = double(NS GT 0.0d0) below_plane = 1.0d0-above_plane ; Something is added to either north (NS < 0) or south (NS > 0) ; hemisphere. The north/south transition is blurred between ; ecl latitudes -hlat0=-5 and +hlat0=+5 by a linear transition. ; Earth above plane of dust: NS > 0 ; in antisolar hemisphere the dust brightness is biased to the south ; zslab = 1 if hlat < -hlat0 (southern hemisphere) ; = 0 if hlat > +hlat0 (northern hemisphere) ; with linear transition between -hlat0 and +hlat0 ; Earth below the plane of dust: NS < 0 ; in antisolar hemisphere the dust brightness is biased to the north ; zslab = 0 if hlat < -hlat0 (southern hemisphere) ; = 1 if hlat > +hlat0 (northern hemisphere) ; with linear transition between -hlat0 and +hlat0 ; Slope of transition region between decl -hlat0 and +hlat0 ; Slope=-1: above plane, Slope=+1: below plane tmp = 1.0d0-2.0d0*above_plane ; Linear transition in hlat across [-hlat0,+hlat0] ; tmp(hlat=-hlat0) = 1.0 (Earth above plane) or 0.0 (Earth below plane) ; tmp(hlat= 0) = 0.5 ; tmp(hlat=+hlat0) = 0.0 (Earth above plane) or 1.0 (Earth below plane) tmp = (tmp*(hlat/(5.0d0/dpr))+1.0d0)/2.0d0 ; Limit to range [0,1] below hlat=-hlat0 and above hlat=+hlat0 tmp >= 0.0d0 tmp <= 1.0d0 ; Increase slab brightness towards Sun. ; Drop off in antisolar hemisphere tmp *= sunward_hemisphere/sin_elo_power_r+antisolar_hemisphere*sin_elo tmp *= 1.00d0-0.25d0*below_plane+2.00d0*cos_elo*sunward_hemisphere zld += 6.00d0*abs(NS)*tmp ENDIF InitVar, camera, 1 zld += ([7.0d0,7.0d0,-3.0d0])[camera-1] zld /= rsun_ecl[2]^power_r END ENDCASE IF silent LE 0 THEN message, /info, $ 'zodiacal light distribution, model'+string(version,format='(F5.2)') CASE n_elements(rrsz) OF 0 : 1 : zld = zld[0] ELSE: zld = reform(zld,rrsz[1:*],/overwrite) ENDCASE RETURN, zld & END