;+
; NAME:
;  vty16_fig4_2
; PURPOSE: (one line)
;  Plot Young 2016, Fig4-2
; DESCRIPTION:
;  Plot Young 2016, Fig4-2
; CATEGORY:
;  VT3D
; CALLING SEQUENCE:
;  vty16_fig4_2
; INPUTS:
;  none
; OPTIONAL OUTPUTS:
;  
; SIDE EFFECTS:
;  creates vty16_fig_14.png, the plot in the Young 2016.
;  vty16_fig3_14_temp_surf.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_fig4_2_plot, dist_sol_au, tod, theta_sva, temp_term_0, temp, p, $
                       name_therminertia, n_per_period

funcname = 'vty16_fig4_2'
run = 'vty16_fig4_2'

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

n_au = n_elements(dist_sol_au)
xval = dist_sol_au & xr=[30,80] & xtit='Heliocentric Distance (AU)'

; control over how things are placed
ypos_theta = [.85, 0.98]
ypos_day   = [.70, 0.80]
ypos_temp  = [.40, 0.68]
ypos_p     = [.05, 0.38]

color_species = ['00B000'xl,'4040B0'xl]
line_species = [2, 4]
ori_species = [30, 90]
thick_species = [2,1]
ori_species = [30, 90.]
;line_species = [4,2]

;-----------------------------------------------------
;  plot the theta's

; the yvaltxt are for output to a text file
yvaltxt = fltarr(n_au, 5)
theta_str = '!4H!X'
i_species = 0
i_therminertia = 0
i_is_volatile = 0
yval = theta_sva[i_species, *, i_therminertia, i_is_volatile,0]
plot, xval, yval, pos=[.12,ypos_theta[0],.98,ypos_theta[1]], /nodata, $
      xr=xr, xstyle=1, xtit='Distance, AU', $;xtickn=REPLICATE(' ', 6), $
      ylo=1, yr=[.1,1e4], ys=1, ytit='Thermal parameter, '+theta_str
; theta_substrate
for i_therminertia = 0, 1 do begin
   yval = theta_sva[0, *, i_therminertia, 0, 0]
   yvaltxt[*,i_therminertia] = yval
   oplot, xval, yval, thick=2
   xyouts, xval[40], yval[40]/4, theta_str+'!DS!N'+' '+name_therminertia[i_therminertia]
endfor
; theta_volatile slab
yval = theta_sva[0, *, 0, 2, 1]
yvaltxt[*,2] = yval
oplot, xval, yval, thick=2
xyouts, xval[2], yval[2]*1.5, theta_str+'!DV!N'
; theta_atmosphere
xyoutstr = ['N!D2!N', 'CH!D4!N']
for i_species = 1,0,-1 do begin
   yval = reform(theta_sva[i_species, *, 0, 0, 2])
   yvaltxt[*,3+i_species] = yval
   color = color_species[i_species]
   line = line_species[i_species]
   oplot, xval, yval, thick=2, co=color, line=line
   xyouts, interpol(xval,yval,.3)+1, .3, theta_str+'!DA!N'+' '+xyoutstr[i_species], co = color
endfor
forprint, xval, yvaltxt[*,0], yvaltxt[*,1], yvaltxt[*,2], yvaltxt[*,3], yvaltxt[*,4], $
          TEXTOUT = funcname + '_table1.txt', $
          comment = string('AU', 'Theta_S_5tiu', 'Theta_S_2100tiu', 'Theta_V_1g', 'Theta_A_N2', 'Theta_A_CH4',for='(A5, 5A16)'), $
          FORMAT = '(F5.1, 5E16.4)'


;-----------------------------------------------------
; Plot the temperatures
i_species = 0
i_therminertia = 0
i_is_volatile = 0
xval = dist_sol_au
yvaltxt = fltarr(n_au, 3*4)
yval = temp_term_0[i_species, *, i_therminertia, i_is_volatile]
yval_min = fltarr(n_au)
yval_max = fltarr(n_au)
plot, xval, yval, pos=[.12,ypos_temp[0],.98,ypos_temp[1]], /nodata, /noerase, $
      xr=xr, xstyle=1, xtickn=REPLICATE(' ', 6), $
      ylo=0, yr=[23,51], ystyle=1, ytit='Temperature, K'
; average equatorial
oplot, xval, yval
; just substrate
i_therminertia = 0
yval_dawn = reform(temp[0, *, i_therminertia, 0, n_per_period/4])
yval_dusk = reform(temp[0, *, i_therminertia, 0, n_per_period*3/4])
for i_au=0,n_au-1 do yval_min[i_au] =min(temp[0, i_au, i_therminertia, 0, *])
for i_au=0,n_au-1 do yval_max[i_au] =max(temp[0, i_au, i_therminertia, 0, *])
yvaltxt[*, 0:3 ] = [ [yval_min],  [yval_dawn],  [yval_dusk],  [yval_max] ]
polyfill, [xval, reverse(xval)], [yval_min,reverse(yval_max)], co='c0c0c0'xl
oplot, xval, yval_min, thick=2
oplot, xval, yval_max, thick=2
polyfill, [xval, reverse(xval)], [yval_dawn,reverse(yval_dusk)], co='a0a0a0'xl
oplot, xval, yval_dawn
oplot, xval, yval_dusk
; plot one more time to get the ticks over the solid polyfill
plot, xval, yval, pos=[.12,ypos_temp[0],.98,ypos_temp[1]], /nodata, /noerase, $
      xr=xr, xstyle=1, xtickn=REPLICATE(' ', 6), $
      ylo=0, yr=[23,51], ystyle=1, ytit='Temperature, K'
; plot the legend
xleg = 70.
dxleg = 5.
yleg = 45.
dylegmin = 0.
dylegmax = 5.
dylegdawn = dylegmin + (dylegmax-dylegmin)/4
dylegdusk = dylegmax - (dylegmax-dylegmin)/4
dylegchar = -.5
xyouts, xleg+0, yleg+(dylegmin+dylegmax)/2., 'No Volatiles ', al=1.0
polyfill, xleg+[ [0,dxleg], reverse([0,dxleg])], yleg+[[dylegmin,dylegmin],[dylegmax,dylegmax]], co='c0c0c0'xl
oplot, xleg+[0,dxleg], yleg+[0,0], thick=2
oplot, xleg+[0,dxleg], yleg+[dylegmax,dylegmax], thick=2
xyouts, xleg+dxleg, yleg + dylegmin + dylegchar, ' Min'
xyouts, xleg+dxleg, yleg + dylegmax + dylegchar, ' Max'
polyfill, xleg+[ [0,dxleg], reverse([0,dxleg])], yleg+[[dylegdawn,dylegdawn],[dylegdusk,dylegdusk]], co='a0a0a0'xl
oplot, xleg+[0,dxleg], yleg+[dylegdawn,dylegdawn]
oplot, xleg+[0,dxleg], yleg+[dylegdusk,dylegdusk]
xyouts, xleg+dxleg, yleg + dylegdawn + dylegchar, ' Dawn'
xyouts, xleg+dxleg, yleg + dylegdusk + dylegchar, ' Dusk'

; CH4 and N2 atm
for i_species = 1, 0, -1 do begin
   color = color_species[i_species]
   line = line_species[i_species]
   ori = ori_species[i_species]
   thick = thick_species[i_species]
   yval_dawn = reform(temp[i_species, *, 0, 1, n_per_period/4])
   yval_dusk = reform(temp[i_species, *, 0, 1, n_per_period*3/4])
   for i_au=0,n_au-1 do yval_min[i_au] =min(temp[i_species, i_au, 0, 1, *])
   for i_au=0,n_au-1 do yval_max[i_au] =max(temp[i_species, i_au, 0, 1, *])
   yvaltxt[*, (4+4*i_species): (7+4*i_species) ] = [ [yval_min],  [yval_dawn],  [yval_dusk],  [yval_max] ]
   polyfill, [xval, reverse(xval)], [yval_min,reverse(yval_max)], thick=thick, orientation=ori, color=color, spac=0.2
   oplot, xval, yval_min, thick=2, line=line, color=color
   oplot, xval, yval_max, thick=2, line=line, color=color
   polyfill, [xval, reverse(xval)], [yval_dawn,reverse(yval_dusk)], thick=thick, orientation=ori, color=color,  spac=0.1
   oplot, xval, yval_dawn, line=line, color=color
   oplot, xval, yval_dusk, line=line, color=color
; plot the legend
   yleg = 38 - 7 * i_species
   xyouts, xleg+0, yleg+(dylegmin+dylegmax)/2., xyoutstr[i_species]+' ', al=1.0
   polyfill, xleg+[ [0,dxleg], reverse([0,dxleg])], yleg+[[dylegmin,dylegmin],[dylegmax,dylegmax]], $
             thick=thick, orientation=ori, color=color,  spac=0.2
   oplot, xleg+[0,dxleg], yleg+[0,0], thick=2, line=line, color=color
   oplot, xleg+[0,dxleg], yleg+[dylegmax,dylegmax], thick=2, line=line, color=color
   xyouts, xleg+dxleg, yleg + dylegmin + dylegchar, ' Min'
   xyouts, xleg+dxleg, yleg + dylegmax + dylegchar, ' Max'
   polyfill, xleg+[ [0,dxleg], reverse([0,dxleg])], yleg+[[dylegdawn,dylegdawn],[dylegdusk,dylegdusk]], $
             thick=thick, orientation=ori, color=color, spac=0.1
   oplot, xleg+[0,dxleg], yleg+[dylegdawn,dylegdawn], line=line, color=color
   oplot, xleg+[0,dxleg], yleg+[dylegdusk,dylegdusk], line=line, color=color
   xyouts, xleg+dxleg, yleg + dylegdawn + dylegchar, ' Dawn'
   xyouts, xleg+dxleg, yleg + dylegdusk + dylegchar, ' Dusk'
endfor 


forprint, xval, $
          yvaltxt[*, 0], yvaltxt[*, 1], yvaltxt[*, 2], yvaltxt[*, 3], $
          yvaltxt[*, 4], yvaltxt[*, 5], yvaltxt[*, 6], yvaltxt[*, 7], $
          yvaltxt[*, 8], yvaltxt[*, 9], yvaltxt[*,10], yvaltxt[*,11], $
          TEXTOUT = funcname + '_table3.txt', $
          comment = string('AU', $
                           'Bare_Min', 'Bare_Dawn', 'Bare_Dusk', 'Bare_Max', $
                           'N2_Min',   'N2_Dawn',   'N2_Dusk',   'N2_Max', $
                           'CH4_Min',  'CH4_Dawn',  'CH4_Dusk',  'CH4_Max', $
                           for='(A5, 12A10)'), $
          FORMAT = '(F5.1, 12F10.2)'

;-----------------------------------------------------
; plot the temperature over a day for selected AU
i_species = 0
i_therminertia = 0
i_is_volatile = 0
n_tod = n_elements(tod)
xval = tod
yvaltxt = fltarr(n_tod, 5*3)
for i_au = 0, 40, 10 do begin
   d = dist_sol_au[i_au]
   d_str = string(d,form='(I2)')+" AU"
   if i_au eq 0 then ytickn='' else ytickn=REPLICATE(' ', 6)
   if i_au eq 0 then ytit='Temperature, K' else ytit = ''
   xpos = linscl( d+[0,10], 30,80,.12,.98)
   plot, xval, xval, pos=[xpos[0],ypos_day[0],xpos[1],ypos_day[1]], /nodata, /noerase, $
         xr=[0,1], xstyle=1, xtickn=REPLICATE(' ', 6), xticks=4, $
         ylo=0, yr=[23,51], ystyle=1, ytickn=ytickn, ytit=ytit
   yval = reform(temp[0, i_au, i_therminertia, 0, *])
   yvaltxt[*,i_au*3/10] = yval
   oplot, xval, yval, th=2
   for i_species = 0, 1 do begin
      color = color_species[i_species]
      line = line_species[i_species]
      yval = reform(temp[i_species, i_au, i_therminertia, 1, *])
      yvaltxt[*,i_au*3/10 + 1] = yval
      oplot, xval, yval, th=2, line=line, color=color
   endfor
   for i=1,3 do oplot, [.25,.25]*i,!y.crange, li=1
   xyouts, 0.57, !y.crange[0]+1, d_str, al=0.5
endfor
forprint, xval, $
          yvaltxt[*, 0], yvaltxt[*, 1], yvaltxt[*, 2], $
          yvaltxt[*, 3], yvaltxt[*, 4], yvaltxt[*, 5], $
          yvaltxt[*, 6], yvaltxt[*, 7], yvaltxt[*, 8], $
          yvaltxt[*, 9], yvaltxt[*,10], yvaltxt[*,11], $
          yvaltxt[*,12], yvaltxt[*,13], yvaltxt[*,14], $
          TEXTOUT = funcname + '_table2.txt', $
          width = 200, $
          comment = string('  TOD', $
                           ' 30AU_Bare', '   30AU_N2', '  30AU_CH4', $
                           ' 40AU_Bare', '   40AU_N2', '  40AU_CH4', $
                           ' 50AU_Bare', '   50AU_N2', '  50AU_CH4', $
                           ' 60AU_Bare', '   60AU_N2', '  60AU_CH4', $
                           ' 70AU_Bare', '   70AU_N2', '  70AU_CH4', for='(A5, 15A10)'), $
          FORMAT = '(F5.1, 15F10.2)'


;-----------------------------------------------------
; plot the pressures
i_species = 0
i_therminertia = 0
i_is_volatile = 0
xval = dist_sol_au
yvaltxt = fltarr(n_au, 2*4)
yval = p[i_species, *, i_therminertia, i_is_volatile]
plot, xval, yval, pos=[.12,ypos_p[0],.98,ypos_p[1]], /nodata, /noerase, $
      xr=xr, xstyle=1, xtit='Heliocentric Distance, AU',$
;      ylo=1, yr=[1e-7,1e2], ystyle=1, ytit='Pressure, !4l!Xbar', ytickv=[1e-6,1e-3,1e0], yticks=2
      ylo=1, yr=[1e-7,1e2], ystyle=1, ytit='Pressure, !4l!Xbar', yticks=9, ytickname=[' ','',' ',' ','',' ',' ','',' ',' ','']
i_therminertia = 0
for i_species = 0, 1 do begin
  ; atm
   color = color_species[i_species]
   line = line_species[i_species]
   ori = ori_species[i_species]
   thick = thick_species[i_species]
   yval_dawn = reform(p[i_species, *, i_therminertia, 1, n_per_period/4])
   yval_dusk = reform(p[i_species, *, i_therminertia, 1, n_per_period*3/4])
   for i_au=0,n_au-1 do yval_min[i_au] =min(p[i_species, i_au, 0, 1, *])
   for i_au=0,n_au-1 do yval_max[i_au] =max(p[i_species, i_au, 0, 1, *])
   polyfill, [xval, reverse(xval)], 1e2< [yval_dawn,reverse(yval_dusk)] > 1e-7, $
             thick=thick, orientation=ori, color=color, spac=0.1
   oplot, xval, yval_dawn, line=line, color=color
   oplot, xval, yval_dusk, line=line, color=color
   polyfill, [xval, reverse(xval)], 1e2< [yval_min,reverse(yval_max)] > 1e-7, $
             thick=thick, orientation=ori, color=color, spac=0.2
   oplot, xval, yval_min, thick=2, line=line, color=color
   oplot, xval, yval_max, thick=2, line=line, color=color
   yvaltxt[*, (0+4*i_species): (3+4*i_species) ] = [ [yval_min],  [yval_dawn],  [yval_dusk],  [yval_max] ]
end

forprint, xval, $
          yvaltxt[*, 0], yvaltxt[*, 1], yvaltxt[*, 2], yvaltxt[*, 3], $
          yvaltxt[*, 4], yvaltxt[*, 5], yvaltxt[*, 6], yvaltxt[*, 7], $
          TEXTOUT = funcname + '_table4.txt', $
          width = 200, $
          comment = string('AU', $
                           'N2_Min',   'N2_Dawn',   'N2_Dusk',   'N2_Max', $
                           'CH4_Min',  'CH4_Dawn',  'CH4_Dusk',  'CH4_Max', $
                           for='(A5, 8A12)'), $
          FORMAT = '(F5.1, 8E12.5)'

tv2im, funcname + '.png', /corner, /true

end