;+
; NAME:
;	smei_zld_dumbbell
; PURPOSE:
;	Calculates analytic dumbbell-shaped correction to the zodiacal light
;	distribution seen by SMEI
; CATEGORY:
;	ucsd/camera/idl/zld
; CALLING SEQUENCE:
	FUNCTION smei_zld_dumbbell, r, $
		degrees 		= degrees		, $
		silent			= silent		, $
		zld_component	= zld_component
; INPUTS:
;	r			array[2] or array[2,n]; type: double
;					Sun-centered ecliptic coordinates (longitude and
;					latitude for points where correction is needed).
; OPTIONAL INPUT PARAMETERS:
;	/degrees	if set, all 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_DUMBBELL and ZLD_MDL__RES_DUMBBELL
;					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:
;	Result		scalar or array[n]
;					brightness of zodiacal light correction in
;					ADUs at nominal gain.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
	@smei_zld_parts.pro
; CALLS:
;	ToRadians, ToDegrees, IsType, EulerRotate, SyncDims, cvsmei
; PROCEDURE:
;	There is a discontinuity of < 0.6 ADUs on the ecliptic.
;
;	The residual map is stored in $SMEIDB/c3zldresidue.fts.gz
; MODIFICATION HISTORY:
;	APR-2009, Paul Hick (UCSD/CASS)
;	APR-2009, Paul Hick (UCSD/CASS)
;		Added /add_residue to add interpolated value from residue
;		zld map to analytic zld correction.
;	DEC-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Couple of bugfixes. Fairly substantial rewrite.
;		Removed /add_residue, and instead added keyword zld_component
;	MAR-2010, John Clover (UCSD/CASS; jclover@ucsd.edu)
;		Bugfix with Sun centered ecliptic domain. Coords were 0
;		and positive in array to rightmost column, and 0 to negative
;		to leftmost column, causing some issue, latitudes are
;		shifted to be -180 to 180 degrees.
;		Added exclusion zone for 0-15 degrees elongation 
;		from Sun and added dropoff starting at 75 degrees
;		elongation and added exclusion outside of 90 degrees
;		elongation from the Sun
;-

InitVar, silent		, 0

get_components = IsType(zld_component,/defined)

IF silent LE 0 THEN	$
	message, /info, 'analytic static dumbbell-shaped zld component'

rpm = ToRadians(degrees=degrees)	; Radians Per Mystery unit
dpm = ToDegrees(degrees=degrees)	; Degrees Per Mystery unit
pi  = !dpi
dpr = 180.0d0/pi					; Degrees Per Radian

szr = size(r)

; Sun-centered ecliptic longitude and latitude in radians

rr = reform(r,szr[1],szr[szr[0]+2]/szr[1])
rr = rr[0:1,*]*rpm

rr[0,*] = AngleRange(rr[0,*], /pi)
ss = rr

; save the elongation before rotation
elocen = acos(cos(ss[0,*])*cos(ss[1,*]))
elocen = reform(elocen,/overwrite)

; North and south hemisphere treated slightly different,
; so split all points up into north and south

south = where( rr[1,*] LE 0.0d0, complement=north )

in_south = south[0] NE -1			; Locations in south present ?
in_north = north[0] NE -1			; Locations in north present ?

; Southern hemisphere

; For sun-centered ecliptic coordinates the x-axis points to
; the Sun; the z-axis to ecliptic north.

; The following rotation puts the x-axis at -15 degrees latitude
; south of the Sun, effectively shifting the southern part of the
; dumbbell 15 degrees south.

IF in_south THEN	$
	ss[*,south] = EulerRotate( [0.0d0,15.0d0,0.0d0]/dpr, rr[*,south] )

; Northern hemisphere

; The following rotation puts the x-axis at +21 degrees latitude
; north of the Sun, effectively shiftinbg the northern part of the
; dumbbell 21 degrees north.

IF in_north THEN	$
	ss[*,north] = EulerRotate( [0.0d0,-21.0d0,0.0d0]/dpr, rr[*,north] )


; Calculate "elongation" and "position angle" in the rotated
; coordinate frames.

; Elongation. Range is [0,+pi]

elo = acos(cos(ss[0,*])*cos(ss[1,*]))
elo = reform(elo,/overwrite)		; Loose leading dummy dimension

; Position angle from ecliptic north (counterclockwise)
; Range is [-pi,pi]

pa = atan(sin(ss[0,*]), tan(ss[1,*]))
pa = reform(pa,/overwrite)			; Loose leading dummy dimension

elo = 60.0*exp(-elo/(16.0/dpr))

pa = abs(pa)

IF in_north THEN pa[north] = pi-pa[north]

xs = abs(rr[0,*])/(6.5d0/dpr)-abs(rr[1,*])/(1.0d0/dpr)+15.0d0
IF in_north THEN xs[north] += 5.0d0

elo *= 0.15d0+0.85d0*exp(-6.0d0*cos(0.6d0*pa)^4.0d0)
pa   = 25.0d0*exp(-acos(cos(2.0d0*ss[0,*])*cos(0.7d0*ss[1,*]))/(10.0/dpr))

near_pole = xs LE 0.0d0
zld = (elo+pa)*near_pole+(elo*exp(-xs/10.0d0)+pa*exp(-xs/8.0d0))*(1-near_pole)

idx = where(elocen gt 75.0/dpr)
if idx[0] ne -1 then zld[idx] *= exp(-(elocen[idx]-(75.0/dpr))*(elocen[idx]-(75.0/dpr))/(120.0/(dpr*dpr)))


;IF in_north THEN zld[north] += 1.0d0
IF in_south THEN zld[south] *= 1.5d0

idx = where(elocen lt 15.0/dpr)
if idx[0] ne -1 then zld[idx] = 0.0
idx = where(elocen gt 90.0/dpr)
if idx[0] ne -1 then zld[idx] = 0.0

IF get_components THEN BEGIN
	CASE (where(zld_component EQ ZLD_MDL__FNC_DUMBBELL))[0] NE -1 OF
	0: zld *= 0
	1: IF silent LE 0 THEN whatis, zld, title='Analytic dumbbell component'
	ENDCASE
ENDIF ELSE zld *= 0

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