;+ ; NAME: ; smei_zld_weekly ; PURPOSE: ; CATEGORY: ; camera/idl/zld ; CALLING SEQUENCE: FUNCTION smei_zld_weekly, t, r , $ elo_angle = elo_angle , $ pa_angle = pa_angle , $ degrees = degrees , $ silent = silent , $ zld_component = zld_component ; INPUTS: ; t array[1]; type: time structure ; r array[2,n]; type: float ; Sun-centered ecliptic longitude and latitude ; OPTIONAL INPUT PARAMETERS: ; elo_angle=elo_angle ; array[n]; type: float ; solar elongation of locations r ; (calculated from r if not specified) ; pa_angle=pa_angle ; array[n]; type: float ; position angle of locations r: angle from ; ecliptic north measured counter-clockwise ; (calculated from r if not specified) ; /degrees if set then angles are in degrees (default: radians) ; zld_component = zld_component ; array; type: integer; default: none ; By default, both components (analytic and residue ; map) are returned. ; The values ZLD_MDL__FNC_WEEKLY and ZLD_MDL__RES_WEEKLY ; can be used in zld_component to extract either one ; of both components. The constants are defined in ; include file smei_zld_parts.pro ; OUTPUTS: ; R array; type: double ; zodiacal brightness at locations r ; INCLUDE: @compile_opt.pro ; On error, return to caller @smei_zld_parts.pro ; CALLS: ; InitVar, ToDegrees, IsTime, TimeSet, TimeGet, TimeUnit, TimeOp ; SyncDims, IsType, big_eph, sgp4_orbit_axis, sphere_distance ; AngleRange, cvsmei ; PROCEDURE: ; Correction consists of an analytic part, plus ; a residual map based on a sequence of 52 weekly residues. ; ; The analytic correction consists of a strictly time-depenent ; part that matches cameras 2 and 3, and corrects for discontinuities ; in camera 3 introduced by different "bad-pixel" masks, and a second ; part that depends on the orientation of the Coriolis orbital plane ; relative to the Sun. ; ; The residue maps are stored in $SMEIDB/c3zldresidue.fts.gz. The ; residual correction is calculated by lineary interpolation between ; two neigbouring weekly maps. ; MODIFICATION HISTORY: ; APR-2009, Paul Hick (UCSD/CASS) ; DEC-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Removed /add_residue, and instead added keyword zld_component. ;- InitVar, silent , 0 get_components = IsType(zld_component, /defined) IF silent LE 0 THEN $ message, /info, 'time-dependent analytic zld component' dpm = ToDegrees(degrees=degrees) ; Degrees Per Mystery unit dpr = 180.0d0/!dpi ; Degrees Per Radian uday = TimeUnit(/day) CASE 1 OF IsTime(t) : tt = t IsTime(t,/orbit): tt = TimeSet(smei=t) ELSE : tt = TimeSet(t) ENDCASE ;doy = TimeGet(tt,TimeUnit(/day),/scalar) ; Integer doy truncated doy = TimeGet(tt,/doy,/scalar) ; Doy with fraction for time of day week = (1+doy/7) < 52 ; ================================== ; Empirical correction to match camera 2 and 3, and ; correction for offsets between cam 3 masks. ; This correction depends only on time, not on location in the sky. zlda = 7.8*cos(doy/55.0+2.4+0.00002*doy*doy)-6.06-0.03*doy ; Matches cams 2 & 3 ; > 972: 2006_056; <= 996: 2006_082 ; > 1136: 2006_264; <= 1250: 2007_021 ; > 1250: 2007_021; <= 1325: 2007_099 ; > 1434: 2007_256; < 1581: 2008_043 ; > 1591: 2008_054; <= 1655: 2008_138 ; >= 1750: 2008_247; <= 1901: 2009_044 ; >= 1915: 2009_059 ts = TimeGet(/_ydoy,TimeGet(tt,/bot,uday),upto=uday,/scalar) ; this is a separate component in Andy's model and isn't taken away until after the zld (arg3 ) is generated. CASE 1 OF '2006_056' LT ts AND ts LT '2006_082': zlda -= 0.50*TimeOp(/subtract,TimeSet(ts),TimeSet('2006_049'),uday) ;(i- 966) '2006_264' LT ts AND ts LT '2007_021': zlda += 0.22*TimeOp(/subtract,TimeSet(ts),TimeSet('2006_263'),uday) ;(i-1136) '2007_021' LT ts AND ts LT '2007_099': zlda += 28.0-0.04*TimeOp(/subtract,TimeSet(ts),TimeSet('2007_020'),uday) ;(i-1250) '2007_256' LT ts AND ts LT '2008_043': zlda += 52.0-0.15*TimeOp(/subtract,TimeSet(ts),TimeSet('2007_256'),uday) ;(i-1435) '2008_054' LT ts AND ts LT '2008_138': zlda += 52.0-0.15*TimeOp(/subtract,TimeSet(ts),TimeSet('2007_256'),uday) ;(i-1435) '2008_240' LT ts AND ts LT '2008_245': zlda -= 15.0 '2008_247' LT ts AND ts LT '2008_261': zlda += 40.0 '2008_261' LT ts AND ts LT '2009_044': zlda += 70.0 '2009_059' LT ts : zlda += 70.0 ELSE: ENDCASE ;=================================== ; Analytic time-dependent correction ; Get geocentric equatorial J2000 coordinates of Sun rsun = big_eph(tt, $ body = 'sun' , $ center = 'earth' , $ /to_sphere , $ /onebody , $ /to_equatorial , $ /degrees ) ; Get geocentric equatorial J2000 coordinates of normal ; to Coriolis orbital plane rnorm = sgp4_orbit_axis(tt,/degrees,/to_sphere,/precess) ; Get elongation (angle from Sun) of normal to orbital plane rnorm = sphere_distance(rsun,rnorm,/degrees) arg = (doy-1)*(360.0/!earth.julianyr) ; Andy originally assumed that the normal tracked the ; mean Sun in longitude, but at an 8.8 degree higher ; latitude (=90-orbit inclination) ;rsun[0] -= 280.0d0+arg ;rsun[1] += 8.8d0 ; Inclination Coriolis orbital axis to equatorial plane ; Elongation of normal to orbital axis (this is Andy's dSMEI) ;rsun = sphere_distance(rsun,[0,0,1],/degrees) arg = (arg < 360.0)-78.25+90.0 arg = AngleRange(arg,/pi,/degrees) ; Sun-centered ecliptic longitude and latitude in degrees szr = size(r) rr = reform(r,szr[1],szr[szr[0]+2]/szr[1]) rr = rr[0:1,*]*dpm ; Lng, lat in degrees rr[0,*] = AngleRange(rr[0,*], /pi, /degrees) ; The elongation and position angle (in degrees) are ; calculated from rr if not provided as input CASE IsType(elo_angle,/defined) OF 0: elo_sun = sphere_distance(rr,rsun[0:1], eps=pa_sun, /degrees) 1: BEGIN elo_sun = elo_angle*dpm pa_sun = pa_angle *dpm END ENDCASE elo_sun = reform(elo_sun,/overwrite); Loose leading dummy dimension pa_sun = reform(pa_sun ,/overwrite) eps = (elo_sun-20.0) > 0.0 arg = abs(arg) arg_small = arg LT 45.0 arg_large = arg GT 135.0 zld = 0.0 zld += (arg_small OR arg_large)*(1.0+0.5*arg_large)*cos(2.0*arg/dpr)*200.0*exp(-(eps+cos(rr[1,*]/dpr))/10.0) zld -= 250.0*exp(-eps/10.0)*(-0.3+rnorm/15.0)*cos(pa_sun/dpr)^2 zld -= 60.0*exp(-eps/15.0-abs((week-41)/4.0)) zld -= 40.0*exp(-eps/ 5.0-abs((week-41)/4.0)) zld *= ((5.0-rr[1,*])/10.0 > 0.0 < 1.0) ; Linear 1 -> 0 from -5 -> +5 degrees ; I had 0.05 instead of 0.0 < 1.0 here, I think I know why but I don't know how ; it got there... zld -= 20.0*exp(-eps/35.0-abs((week- 9)/4.0)) zld -= 20.0*exp(-eps/ 5.0-abs((week- 9)/4.0)) zld -= 5.0*exp(-eps/40.0-abs((week-41)/4.0)) zld += 40.0*exp(-eps/ 5.0-abs((week-41)/4.0)) zld -= 2.0 zld += zlda IF get_components THEN BEGIN CASE (where(zld_component EQ ZLD_MDL__FNC_WEEKLY))[0] NE -1 OF 0: zld *= 0 1: whatis, zld, title='Analytic time-dependent correction' ENDCASE ENDIF ; Add a correction from the residual zld map residue_file = filepath(root=getenv('SMEIDB'),'c3zldresidue.fts.gz') week1 = readfits(residue_file, ext= fix(week) ) week2 = readfits(residue_file, ext=(fix(week) mod 52)+1) CASE IsType(week1,/array) OF 0: IF silent LE 0 THEN message, /info, 'no weekly residue added' 1: BEGIN ; Extract two weekly residue maps bracketing the input time ; residue is [720,360,2] array. residue = [ [[week1]], [[week2]] ] ; Convert helioecliptic lng/lat to array indices rr = cvsmei(from_helioecliptic=rr, /degrees, /to_map, mode='eclsun') ; Interpolate on residue map and add to zld. week = replicate(week, size(rr[0,*],/dim)) ;zld_residue = interpolate(residue, reform(rr[0,*]), reform(rr[1,*]), week, missing=bad) ;bad is undefined, it's the same to leave it out. should it be undefined? zld_residue = interpolate(residue, rr[0,*], rr[1,*], week) IF get_components THEN BEGIN CASE (where(zld_component EQ ZLD_MDL__RES_WEEKLY))[0] NE -1 OF 0: zld_residue *= 0 1: IF silent LE 0 THEN whatis, zld_residue, title='Residual weekly correction' ENDCASE ENDIF ELSE zld_residue *= 0 zld += zld_residue END ENDCASE residueavg = readfits(residue_file,ext=0) CASE IsType(residueavg,/array) OF 0: IF silent LE 0 THEN message, /info, 'no static residue map added' 1: BEGIN ; Convert helioecliptic lng/lat to array indices ;rr = cvsmei(from_helioecliptic=rr,/degrees, /to_map, mode='eclsun') ; Interpolate on residue map and add to zld. zld_residueavg = interpolate(residueavg, rr[0,*], rr[1,*]) IF get_components THEN BEGIN CASE (where(zld_component EQ ZLD_MDL__RES_DUMBBELL))[0] NE -1 OF 0: zld_residueavg *= 0 1: IF silent LE 0 THEN whatis, zld_residueavg, title='Residual dumbbell component' ENDCASE ENDIF ELSE zld_residueavg *= 0 zld += zld_residueavg END ENDCASE CASE size(r,/n_dim) GT 1 OF 0: zld = zld[0] 1: SyncDims, zld, size=[szr[0]-1,szr[2:szr[0]],szr[0]+1,szr[szr[0]+2]/szr[1]] ENDCASE RETURN, zld & END