;+
; NAME:
;  vty16_sol_terms_diurnal
; PURPOSE: (one line)
;  return the diurnal solar terms in sinusoidal expansion
; DESCRIPTION:
;  return the diurnal solar terms such that
;    s[i,lon] = sol_terms[i,0] + Re( sol_terms[i,1]*exp(_i*  p) )$
;                          + Re( sol_terms[i,2]*exp(_i*2*p) ) + 
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;   sol_terms = vty16_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol, n_terms)
; INPUTS:
;  dist_sol_au : sub-solar distance, in AU, float
;  albedo : hemispheric albedo, unitless, float or float[n_loc]
;  lat : latitude, radian, float or float[n_loc]
;  h_phase0 : hour angle at zero phase
;  lat_sol : subsolar latitude, radian, float
;  n_terms : number of terms in the expansion
; OPTIONAL INPUTS:
;  sol_norm_1au - Normal incident solar flux at 1 AU (scalar, erg/cm^2/s)
;    Default 1367.6e3
; OUTPUTS:
;  sol_terms : the complex diurnal solar terms in a sinusoidal expansion, 
;    ergcm^2/s, complex[n_loc, n_term+1] or complex[n_term+1]
; OPTIONAL INPUTS:
;  h_max - Hour angle at sunset, radian, float or float[n_loc]
; PROCEDURE:
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;  Resubmitted to Icarus.
;  Eq 3.2-4a,  3.2-4b,  3.2-4c
; MODIFICATION HISTORY:
;  Written 2012 Jan 8, by Leslie Young, SwRI
;  2012 Apr 19 LAY, change offset to h_phase0
;  2016 Mar 22 LAY. Modified for inclusion in vty16 library.
;-

function vty16_sol_terms_diurnal,dist_sol_au, albedo, lat, h_phase0, lat_sol, n_term, $
                           sol_norm_1au = sol_norm_1au, h_max=h_max
    
    if n_elements(sol_norm_1au) eq 0 then sol_norm_1au = 1367.6e3
    _i = complex(0.,1.)

    coslat = cos(lat)
    sinlat = sin(lat)

    sinlat0 = sin(lat_sol)
    coslat0 = cos(lat_sol)

    sol_ss = sol_norm_1au * (1-albedo) / dist_sol_au^2 ; erg/cm^2/s

    n_loc = n_elements(albedo * lat)
    s_2d = (size(albedo*lat, /n_dim) gt 0)
    n_time = n_elements(lat0)
    
    sol_terms = complexarr(n_loc, n_term+1)
    h_max = acos((-1 > (-tan(lat)*tan(lat_sol)) < 1))
    cosh = cos(h_max)
    sinh = sin(h_max)

    indx = where(h_max gt 0., nindx)
    if nindx eq 0 then return, sol_terms ; all zero

    sol_terms[*,0] = ( 0. > ( (sinlat*sinlat0)*h_max + (coslat*coslat0)*sin(h_max)) / !pi )  * sol_ss
    if n_term ge 1 then begin
       sol_terms[*,1] = 2*sinlat*sinlat0*sinh/(!pi) + $
                    (coslat*coslat0)*(h_max + cosh*sinh)/(!pi)
       sol_terms[*,1] *= exp( _i * h_phase0)  * sol_ss
    endif
    for m = 2, n_term do begin
       cosmh = cos(m*h_max)
       sinmh = sin(m*h_max)
       sol_terms[*,m] = sinlat*sinlat0*sinmh*2/(m*!pi) + $
                    (coslat*coslat0)*(m*cosh*sinmh-sinh*cosmh)*2/(!pi*(m^2-1))
       sol_terms[*,m] *= exp( _i * m * h_phase0)  * sol_ss
    endfor

    if s_2d then return, sol_terms else return, reform(sol_terms)

end

pro vty16_sol_terms_diurnal_test

    deg = !pi/180.
    _i = complex(0,1.)

    dist_sol_au = 1.
    albedo = 0.
    n_terms = 7
    p = findgen(360.)*deg

    h_phase0 = 90*deg

    lat = 0.*deg & lat_sol = 00.*deg
    ss = ( sin(lat)*sin(lat_sol) + cos(lat)*cos(lat_sol)*cos(p + h_phase0) > 0.)
    sol_terms = vty16_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol, n_terms, $
                                        sol_norm_1au=1., h_max=h_max)
    plot, p/deg, ss, xr=[0,360], /xs, yr=[-.5,1.5], th=2
    se = cos(0*p) * sol_terms[0]
    oplot, p/deg, se, co=x11rgb('green')
    for m=1, n_terms do    se += re(exp(_i*m*p) * sol_terms[m])
    oplot, p/deg, se, co=x11rgb('blue')
    se2 = vty16_thermwave(sol_terms, p, 0., 1.,0.)
    oplot, p/deg, se2, co=x11rgb('red'), li=2
    print, lat/deg, lat_sol/deg, h_max/deg, re(sol_terms), $
           for='(2F6.2,F7.2,8F8.3)'


end