;+
; NAME:
;  vty16_fig3_1
; PURPOSE: (one line)
;  Plot Young 2016, Fig3.1
; DESCRIPTION:
;  Plot Young 2016, Fig3.1
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
; CATEGORY:
;  VT3D
; CALLING SEQUENCE:
;  vty16_fig3_1, phase, flux_sol, flux_sol_1, flux_sol_7
; INPUTS:
;  none
; OPTIONAL OUTPUTS:
;  phase: phase of the period running from 0 to 2 pi radian
;  flux_sol : 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 files 
;  vty16_fig3_1.png - as in Young 2016
;  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_1, phase, flux_sol, flux_sol_1, flux_sol_7

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

;-------------------------
; set up for correct plotting
;-------------------------
loadct, 0
device,/decompose
!p.background = 'ffffff'xl
!p.color=0
!p.charsize=1.5

;-------------------------
; constants
;-------------------------
physconstants    ; load physical constants from layoung library
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 = '(A-20, ": ", A-20,"=",F15.6," ",A-10)'
forme = '(A-20, "_sol_term_bare: ", A-20,"=",E15.3," ",A-10)'
formi = '(A-20, "_sol_term_bare: ", A-20,"=",I8,"        ",A-10)'
forms = '(A-20, "_sol_term_bare: ", A-48)'

;-------------------------
; Insolation on an bare area
;-------------------------
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 = vt3d_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 = vt3d_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol[0], 1)
flux_sol_t_7 = vt3d_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 = vt3d_solwave(flux_sol_t_1,phase)
flux_sol_7 = vt3d_solwave(flux_sol_t_7,phase)

; Plot it
window, 0, xs=840, ys=514
!p.color = 0
plot, phase/deg, flux_sol, /nodata, $
      xtit='Phase, deg', xstyle=3, xticks=4,xminor=6,$
      ytit='Insolation, erg cm!U-2!N s!U-1!N', yrange=[-100,400]*15.
; Exact as thick, gray
oplot, phase/deg, flux_sol, color='a0a0a0'xl, thick=6
oplot, 270-[5,45], 5000+94+[0,0], color='a0a0a0'xl, thick=6
xyouts, 270, 5000, 'Numerical'
; M=1 as dashed
oplot, phase/deg, flux_sol_1, line=0, thick=2
oplot, 270-[5,45], 4500+94+[0,0], line=0, thick=2
xyouts, 270, 4500, 'Approximation, M=1'
; M=7 as dot-dashed
oplot, phase/deg, flux_sol_7, line=2, thick=2
oplot, 270-[5,45], 4000+94+[0,0], line=2, thick=2
xyouts, 270, 4000, 'Approximation, M=7'
tv2im, funcname+'.png', /corner

forprint, phase/deg, flux_sol, flux_sol_1, flux_sol_7, $
          TEXTOUT = funcname + '.txt', $
          comment = string('phase', 'numeric', 'M=1', 'M=7', for='(A7, 3A8)'), $
          FORMAT = '(F7.2, 3F8.2)'

end