;+ ; NAME: ; vt3d_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 = vt3d_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 OUTPUTS: ; 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 vt3d library. ;- function vt3d_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 vt3d_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 = vt3d_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 = vt3d_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