;+
; NAME:
;  vty16_fig3_14
; PURPOSE: (one line)
;  Plot Young 2016, Fig3-14
; DESCRIPTION:
;  Plot Young 2016, Fig3-14
; CATEGORY:
;  VT3D
; CALLING SEQUENCE:
;  vty16_fig3_14
; INPUTS:
;  none
; OPTIONAL OUTPUTS:
;  temp_0_map float[n_lon,n_lat,n_t]
;   surface temperature as a function of longitude, latitude and 
;   time (aka sub-solar longitude)
;   latitudes at cell center are [-88, -84, ... 84, 88]
;   longitudes at cell center are [2, 6, ... 354, 358]
;   sub-solar longitudes are [0, -1, -2, ... -1079, -1080]
;     (mod 360, this is [0, -1, -2, ... -359, 0]
;   
; SIDE EFFECTS:
;  creates the plot in the Young 2016.
;   vty16_fig3_14_043.png, vty16_fig3_14_087.png, vty16_fig3_14_0167.png, vty16_fig3_14_cb.png
;    vty16_fig3_14.fits - output of vt3s_thermpro
; MODIFICATION HISTORY:
;  Written 2013 Jun 5, by Leslie Young, SwRI
;  Based on ms_plutoseason/figs/__working_copy/_mimas.pro
;  2016 Mar 24 LAY. Modified for inclusion in vty16 library.
;-

pro vty16_fig3_14, temp_0_map

;-------------------------
; this function name
;-------------------------
funcname = 'vty16_fig3_14'
run = 'vty16_fig3_14'

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

;-------------------------
; Options for this run
;-------------------------
n_term = 7 ; number of terms for solar decomposition
n_lat = 45
n_per_period = 360 ; counter

;-------------------------
; constants
;-------------------------
physconstants                ; define standard physical constants in !phys
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
km_per_cm = 1e5                ; km per cm
_i = complex(0,1)              ; sqrt(-1)

;-------------------------
; surface grid
;-------------------------
n_lon = n_lat * 2
locgrid_rect, n_lon, n_lat, n_loc, lon, lat, angarea_delta

;-------------------------
; Defining if locations are in Region 1 or Region 2
;-------------------------
rlon = 96.87 * deg
rlat = 40.9 * deg
clon = -90.8 * deg ; east longitude of center
clat = 0. * deg    ; latitude of center
dlat = lat-clat
dlon = acos(cos(lon - clon))
region2 = ( ( (dlat/rlat)^2 + (dlon/rlon)^2 ) lt 1.)
region1 = 1-region2
is_region2 = where(region2)

;-------------------------
; surface quantities
;-------------------------
albedo = 0.60 * region1 + 0.59 * region2
emis   = 1.00
therminertia = (9. * region1 + 66. * region2) *1e3
dens = 0.934 
specheat = 0.8 * 1e7; 0.8 J/K/g
thermcond = therminertia^2 / (dens * specheat)
radius = 198. * km_per_cm ; radius in cm

;-------------------------
; Initialize insolation on an bare location
;-------------------------

; Date__(UT)__HR:MN     Ob-lon Ob-lat Sl-lon Sl-lat               r        rdot
;******************************************************************************
; 2010-Feb-13 19:10     127.42 -10.23 147.59   2.21  9.487299015475  -7.6211720
;   Ob-lon Ob-lat =
;    Apparent planetographic ("geodetic") longitude and latitude (IAU2006 model)
; of the center of the target seen by the OBSERVER at print-time. This is NOT
; exactly the same as the "sub-observer" (nearest) point for a non-spherical
; target shape, but is generally very close if not an irregular body shape. Light
; travel-time from target to observer is taken into account. Latitude is the
; angle between the equatorial plane and the line perpendicular to the reference
; ellipsoid of the body. The reference ellipsoid is an oblate spheroid with a
; single flatness coefficient in which the y-axis body radius is taken to be the
; same value as the x-axis radius. Positive longitude is to the WEST.
; Units: DEGREES
; 
;  
;  Sl-lon Sl-lat =
;    Apparent planetographic ("geodetic") longitude and latitude of the Sun
; (IAU2006) as seen by the observer at print-time.  This is NOT exactly the same
; as the "sub-solar" (nearest) point for a non-spherical target shape, but is
; generally very close if not an irregular body shape. Light travel-time from Sun
; to target and from target to observer is taken into account.  Latitude is the
; angle between the equatorial plane and the line perpendicular to the reference
; ellipsoid of the body. The reference ellipsoid is an oblate spheroid with a
; single flatness coefficient in which the y-axis body radius is taken to be the
; same value as the x-axis radius. Positive longitude is to the WEST.
; Units: DEGREES

period = .9424218 * s_per_day
freq = 2*!dpi/period
lat_sol = 2.21*deg ; subsolar latitude in radians as a function of time
dist_sol_au = 9.487299
; absorbed solar energy decomposed into constant and 
; first Fourier term, such that 
;   sol[i_loc, phase] = sol_t[i_loc, 0] + 
;                       re( sol_t[i_loc, 1] * exp( _i * phase ) )
; sol_terms_diurnal,dist_sol_au, albedo, lat, lon_sol_ref, lat_sol, n_term)
;h_phase0 = lon - (360-147)*deg  ; 
;lon_sol_phase0 = -148*deg ; sub-solar EAST longitude at zero phase
lon_sol_phase0 = 0.*deg ; sub-solar EAST longitude at zero phase
lat_obs =  -10.23*deg ; sub-observer latitude

h_phase0 = lon - lon_sol_phase0 ; [n_loc]
flux_sol_t = vt3d_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol, 7)
       ; erg cm^-2 s^-1 : complex [n_loc, n_term+1]
sol_ss = sol_norm_1au * (1-albedo) / dist_sol_au^2 ; erg/cm^2/s

;-------------------------
; Calculate the temperature terms
;-------------------------
flux_int = 0.
dtemp_dzeta = 0.
;temp_t = vt3d_temp_terms_bare(flux_sol_t,flux_int,emis,freq,therminertia)
temp_t = vt3d_temp_terms_bare_iter(flux_sol_t,flux_int,emis,freq,therminertia, thermcond)
dtemp_dzeta = 0.
del_typ = 0.25
zeta = -findgen(10/del_typ+1) * del_typ ; [n_z]
n_z = n_elements(zeta)
z_skin = vt3d_skindepth(dens, specheat, thermcond, freq)
temp_phase0 = vt3d_thermwave(temp_t, 0., dtemp_dzeta, -1, zeta)

;============================= 
; Begin the iteration
;============================= 

;-------------------------
; things that don't change with time
; the interior alpha and beta matrix elements
;  alpha_i, beta_i, eta_i
;-------------------------
del_a = replicate(del_typ, n_z)     ;normalized layer width Above, float[n_z]
del_b = replicate(del_typ, n_z)     ;normalized layer width Below, float[n_z]
del = replicate(del_typ, n_z)       ; normalized layer width, float[n_z]
del[0] = del[0]/2.

phi_s = vt3d_dfluxdtemp_substrate(freq, therminertia)                  ; erg cm^-2 s^-1 K^-1 : float, scalar
gamma_J = 0. ; lower boundary, scalar

; phase, defined to 0 at start 
n_per_period = 360 ; counter
n_t = 3* n_per_period + 1 ; counter
delta_t = period/n_per_period ; change in time per time step, float, scalar
tau = delta_t * freq       ; change in phase per time step, radian, scalar
phase = findgen(n_t) * tau ; phase at start of time step
phase_tavg = phase + 0.5 * tau ; phase at middle of time step
lon_sol_arr = lon_sol_phase0 - phase ; solar lon at start of time step

alpha_i = vt3d_alpha_interior(tau, del, del_a) ;& alpha = [0., alpha_i] ; float[n_z-1]
beta_i =  vt3d_beta_interior(tau, del, del_b) ;& beta = [0.0,beta_i]      ; float[n_z-1]
eta_i = 1 - alpha_i - beta_i ;& eta = 1 - alpha - beta

;-------------------------------------------
;initialize surface and interior temperatures
;-------------------------------------------
temp_0 = temp_phase0[*,0] ; [n_loc]
temp_i = temp_phase0[*,1:n_z-1] ; [n_loc, n_z-1]

;-------------------------------------------
; assign arrays to save temperatures
;-------------------------------------------
temp_mat = fltarr(n_loc, n_t, n_z) ; temp at time_n
temp_mat[*,0,0] = temp_0
temp_mat[*,0,1:n_z-1] = temp_i
flux_sol_tavg_arr = fltarr(n_t,n_loc) ; insolation from time_n to time_n+delta_t

spp_tridc = vt3d_cn_tridc(alpha_i, beta_i, spp_indx)

for i_t = 0L, n_t-2 do begin
   ; Each loop is from time time_n to time_n + delta_t
   ; The temperature is defined at the start of the loop at time_n
   ; The step finds the new temperature at time_n + delta_t
   phi_e = vt3d_dfluxdtemp_emit(emis, temp_0) ; erg cm^-1 s^-1 K^-1 : float[n_loc]

   ; total derivative of flux w/r.to temperature at the surface
   ;  in this timestep 
   phi_tot = del[0]*phi_s/tau + phi_e/2. ; erg cm^-1 s^-1 K^-1 : float[n_loc]

   ; the surface beta matrix elements
   beta_0 = (phi_s/del_b[0])/phi_tot ; float[n_loc]

   ; the surface gamma matrix element
   phase_mid = phase_tavg[i_t] ; (i_t + 0.5)*tau
   sol_tavg = sol_ss * vt3d_solar_mu(lat, lon, lat_sol, lon_sol_phase0 - phase_mid)
   flux_sol_tavg_arr[i_t] = sol_tavg
   ; Solar forcing array element, K
   gamma_0 = (sol_tavg - emis*!phys.sigma*( temp_0)^4.)/ phi_tot 

   ; Replace temp_0, temp_i at timestep i_t with values at time i_t+1.
   ; The figure in Young 2015 was made with vt3d_step_expl_nloc.
   ; To use other time-step schemes, uncomment one of the sets of
   ; lines below
   ; vt3d_step_expl_nloc calls the EXPLICIT time step
   ; This is what is in Young 2016
   vt3d_step_expl_nloc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, $
                          temp_0, temp_i
   ; vt3d_step_cn_nloc calls the Crank-Nicholson time step
   ; without passing the tri-diagonal LU decomposition
   ; Uncomment the next two lines
   ; vt3d_step_cn_nloc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, $
   ;                    temp_0, temp_i
   ; vt3d_step_cn_nloc calls the Crank-Nicholson time step
   ; with passing the tri-diagonal LU decomposition
   ; Uncomment the next two lines
   ; vt3d_step_cn_trisol_nloc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, $
   ;     spp_tridc, spp_indx, temp_0, temp_i

   temp_mat[*,i_t+1,0] = temp_0
   temp_mat[*,i_t+1,1:n_z-1] = temp_i

endfor

;----------------------------
; PLOT
;----------------------------
; load the color table
loadct, 13, rgb=rgb
col_r = rgb[*,0]
col_g = rgb[*,1]
col_b = rgb[*,2]
col_true = col_b*256L*256L + col_g*256L + col_r

; temperature range to plot
temp_min = 65.
temp_max = 95.

wlon_sol_plot_deg = [167, 87, 43.] ; Howett et al. 2011, Table 1
wlon_obs_plot_deg = [147., 180, 83] ; Howett et al. 2011, Fig 1
nplot = n_elements(wlon_sol_plot_deg)
temp_0_map = reform(temp_mat[*,*,0],n_lon,n_lat,n_t)
writefits, run+'.fits', temp_0_map

for iplot = 0, nplot - 1 do begin
   sol_plot = -wlon_sol_plot_deg[iplot] * deg
   off_t = ( (lon_sol_phase0 - sol_plot) / tau ) 
   i_t = n_t - n_per_period + off_t
   temp_0_b = round(linscl(temp_mat[*,i_t,0],temp_min,temp_max,0., 255., /clip))
   wlon_sol_deg = (-lon_sol_arr[i_t]/deg) mod 360
   mtitle = string(wlon_sol_deg, for='(I03, " deg")')
   mapdisp, col_true[temp_0_b], /wlo,mtitle=mtitle, /tru
   ALTAZ2HADEC, 30., findgen(361), lat_obs/deg, ha, dec
   oplot, -(wlon_obs_plot_deg[iplot] + (  (ha+180) mod 360 - 180 )), dec, co=255, line=2
   tv2im, run+'_'+string(wlon_sol_deg, for='(I03)')+'.png', /co, /tr
wait,3
end

; plot the scale bar
window, 0, xs=100, ys=339
device, dec=0
loadct,0
erase,'ffffff'xl
loadct,13
n_scale = 300
pos = [10,10,10,10]+[0,0,30,n_scale]
temp_scale =  replicate(1B, 30) # arrgen(temp_min, temp_max, (temp_max - temp_min)/(n_scale-1))
b = round(linscl(temp_scale,temp_min,temp_max,0., 255., /clip)) ; byte
tv,b,pos[0], pos[1]
plot, b, /nodata, pos=pos,/noerase,/dev, /xticks, /yticks, /xmin, /ymin, xtickname=REPLICATE(' ', 2), ytickname=REPLICATE(' ', 2)
axis,/yax, yr=[temp_min,temp_max], ys=1, ytit='Temperature (K)'
tv2im, run+"_cb.png", /co, /tr

end