;+
; NAME:
;  vty16_fig3_3
; PURPOSE: (one line)
;  Plot Young 2016, Fig3.3
; DESCRIPTION:
;  Plot Young 2016, Fig3.3
; CATEGORY:
;  VT3D
; CALLING SEQUENCE:
;  vty16_fig3_3, tod_vt3d, tsurf_vt3d, phase, temp, temp_1, temp_7, temp
; INPUTS:
;  none
; OPTIONAL OUTPUTS:
;  tod_vt3d: phase of the period running from 0 to 2 pi radian
;  tsurf_vt3 : absorbed solar flux, erg/cm^2/s
;  flux_sol_1 : approximate absorbed solar flux, M=1, erg/cm^2/s
;  flux_sol_7 : approximate absorbed solar flux, M=7, erg/cm^2/s
; SIDE EFFECTS:
;  creates vty16_fig3_1.tif and vty16_fig3_1.txt
; MODIFICATION HISTORY:
;  Written 2012 Sep 29, by Leslie Young, SwRI
;  2016 Mar 22 LAY. Modified for inclusion in vty16 library.
;-

pro vty16_fig3_3, tod_vt3d, tsurf_vt3d, phase, temp, temp_1_zeta0, temp_7_zeta0, temp_7_iter_zeta0

;-------------------------
; this function name
;-------------------------
funcname = 'vty16_fig3_3'

;-------------------------
; constants
;-------------------------
physconstants
deg = !pi/180.
sol_norm_1au = 1367.6e3      ; solar constant at 1 AU, erg/cm^2/s
s_per_day = 24. * 3600.
s_per_year = 365.25 * s_per_day
cm_per_au = 149597870.691e5    ; cm per AU
km_per_au = 149597870.691      ; km per AU
_i = complex(0,1)

;-------------------------
; formats for printing
;-------------------------
formf = '("' + funcname + ': ", A-20,"=",F15.6," ",A-10)'

;-------------------------
; Insolation on an asteroid
;-------------------------
period = .9424218 * s_per_day
n_phase = 720 ; number per period
phase = dindgen(n_phase) * 2.d * !dpi / n_phase
lat = 30. * deg ; latitude in radians
lon = 0. ; east longitude in radians
lat_sol = replicate(2.24*deg, n_phase) ; subsolar latitude in radians as a function of time
; subsolar east longitude in radians; decreases with time for
; prograde objects
lon_sol = (!pi/2 - phase)

; Calculate the exact insolation
; mu0 is the cosine of the solar illumination angle
mu0 = vty16_solar_mu(lat, lon, lat_sol, lon_sol)
albedo = 0.6
dist_sol_au = 9.487299
flux_sol_ss = (1-albedo) * sol_norm_1au / (dist_sol_au^2)
flux_sol = flux_sol_ss * reform(mu0)

; Calculate the sinusoidal expansion
h_phase0 = lon - lon_sol[0]
; flux_sol_t_1 and flux_sol_t_7 are the M=1 and M=7 arrays of coefficient
flux_sol_t_1 = vty16_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol[0], 1)
flux_sol_t_7 = vty16_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol[0], 7)
; flux_sol_1 and flux_sol_7 are the M=1 and M=7 approximations
flux_sol_1 = vty16_solwave(flux_sol_t_1,phase)
flux_sol_7 = vty16_solwave(flux_sol_t_7,phase)

; Calculate the temperature terms
flux_int = 0.
emis = 1.
freq = 2.d * !dpi / period
therminertia = 16e3 ; Mimas region 1
dens = 1.0
specheat = 3.e6
thermcond = therminertia^2/(dens * specheat)
z_skin = vty16_skindepth(dens, specheat, thermcond, freq)

temp_eq = (flux_sol / (!phys.sigma * emis) ) ^ 0.25

temp_terms_1 = vty16_temp_terms_bare(flux_sol_t_1,flux_int,emis,freq,therminertia)
temp_terms_1_iter = vty16_temp_terms_bare_iter(flux_sol_t_1,flux_int,emis,freq,therminertia, thermcond)
temp_terms_7 = vty16_temp_terms_bare(flux_sol_t_7,flux_int,emis,freq,therminertia)
temp_terms_7_iter = vty16_temp_terms_bare_iter(flux_sol_t_7,flux_int,emis,freq,therminertia, thermcond)
dtemp_dzeta = 0.
temp_1_zeta0 = vty16_thermwave(temp_terms_1, phase, dtemp_dzeta, z_skin, 0.)
temp_7_zeta0 = vty16_thermwave(temp_terms_7, phase, dtemp_dzeta, z_skin, 0.)
temp_7_iter_zeta0 = vty16_thermwave(temp_terms_7_iter, phase, dtemp_dzeta, z_skin, 0.)


temp_0 = vty16_temp_term0_bare(float(flux_sol_t_1[0]), flux_int, emis)
phi_e = vty16_dfluxdtemp_emit(emis, temp_0)
phi_s = vty16_dfluxdtemp_substrate(freq, therminertia)
help, temp_0, phi_e, phi_s
print, 'Theta', 4 * phi_s / phi_e, form=formf

vt3d_thermpro,tsurf_vt3d,tod_vt3d,tdeep_vt3d,teq_vt3d,balance_vt3d,tdarr,zarr0, $
            rhel=dist_sol_au,         $          ; heliocentric distance, AU
            alb=albedo,           $              ; bolometric albedo
            ti=therminertia,             $       ; thermal inertia, erg-cgs. 
            rot=period/s_per_day,           $    ; Rotation period, days
; Optional physical keyword parameters:
            emvty=emis,       $                ; emissivity                   (default 1.0)
            rho=dens,           $              ; density, g cm-3.  Can be an array, giving the density of each slab. (default 1.0)
            cp=specheat,               $       ; specific heat, erg g-1 K-1 (default H20)
            plat=lat,         $         ; latitude (radians)           (default 0.0)
            sunlat=lat_sol[0],     $              ; Subsolar latitude (radians)  (default 0.0)
            heatflow=flux_int,   $                  ; $                ; Endogenic heat flow, erg cm-2 s-1 (default 0.0)
; Optional simulation keyword parameters:
            ntinc=1000, $           ; Time increments per day                      (default 5000)
            nday=20                    ; Number of days per run                       (default 4)

window, 0, xs=840, ys=514
!p.color = 0
plot, phase/deg, temp_eq, /nodata, $
      xtit='Phase, deg', /xstyle, xticks=4,xminor=6,$
      ytit='Temperature, K', yrange=[60,92], /ystyle
; thick gray = numeric
reorder_by_phase, tod_vt3d-h_phase0+!pi, tsurf_vt3d, x, y
oplot, x/deg, y, co='a0a0a0'xl, th=4
; dashed M=1
oplot, phase/deg, temp_1_zeta0, li=2
; dot-dashed M=7
oplot, phase/deg, temp_7_zeta0, li=3
; dot-dot-dot-dash M=7, iter
oplot, phase/deg, temp_7_iter_zeta0, li=4
tv2im, funcname + '.tif', /corner


end