;+
; NAME:
;  vty16_pluto_mssearch_resub_func
; PURPOSE: (one line)
;  Make the data for a Pluto movie for Young 2013 ApJL resubmission
; DESCRIPTION:
;  Make the data for a Pluto movie for Young 2013 ApJL resubmission
; CATEGORY:
;  VT
; CALLING SEQUENCE:
;  vty16_fig5_7_still, run, res, yr_still
; INPUTS:
;  run : name of this run.
;  res: structure with the results of the run, suitable for plotting.
;    See vty16_pluto_mssearch_resub_func
;  yr_still : year for which to make the plot
; SIDE EFFECTS:
;    write a file fig/<run>/mov/<run>_<iii>.png
;    where <iii> is the time index relative to the start of the last period
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2012 October 13, by Leslie Young, SwRI (pluto_dps12_2.pro)
;  Oct 21 2012.  Make function (pluto_mssearch_func.pro)
;  Dec 13 2012.  Add start phase as options
;  2016 Apr 06 LAY. Modified for inclusion in vty16 library.
;-

pro vty16_fig5_7_still, run, res, yr_still, saveplot=saveplot, openplot=openplot, showrise=showrise
    white = 'ffffff'xl
    red = '0000ff'xl
    orange = '00a5ff'xl
    green = '00ff81'xl
    gold = total([255, 215, 00] * [1L, 256L, 65536L]) ; 255 215       0             gold1
    blue = total([0, 191, 255] * [1L, 256L, 65536L])
    purple = total([255,48,155] * [1L, 256L, 65536L]) ; 155  48 255         purple1

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

    ;-------------------------
    ; constants
    ;-------------------------
    physconstants
    deg = !pi/180.
    day = 24. * 3600.
    year = 365.25 * day
    km = 1e5                    ; cm per km
    au = 149597870.691e5        ; cm per AU
    sol_norm_1au = 1367.6e3     ; solar constant at 1 AU, erg/cm^2/s
    time_day =  6.3865509d * day ; %%%
    time_year = 247.92065d * year ; %%%
    gravacc = 80.                 ; cm/s^2
    radius = 1163.e5              ; Sicardy et al 2011

    xs = 10 * 120
    ys = 6.25 * xs/10
    xm = round(0.3 * xs)
    ym = ys - xm

    if !d.window ne 2 or !d.x_size ne 1280 or !d.y_size ne 700 then $
       window, 2, xs=1280, ys=700 ; 1280 x 700 is allowed
    !p.color = 0 & !p.background = 'ffffff'xl
    charsizeold = !p.charsize
    !p.charsize = 2

    lat = res.lat
    nlat = -lat

    ; get the indices for solstice and equinox
    n_time_per_period = res.n_time_per_period
    n_time = n_time_per_period  ; **************************
    yr = res.yr
    indx_yr = lindgen(n_elements(yr))
    l0 = res.lat0
    nl0 = -l0
    r = res.r
    ang = res.ang

    l0_south = min(l0, i) & yr_ssolstice = yr[i]
    l0_north = max(l0, i) & yr_nsolstice = yr[i]
    indx1 = where(l0[1:n_time_per_period-1]-l0[0:n_time_per_period-2] gt 0, compl=indx2)
    l0_eq = min(l0[indx1]^2, i1) & indx_eq1 = indx1[i1]
    l0_eq = min(l0[indx2]^2, i2) & indx_eq2 = indx2[i2]
    if r[indx_eq1] le r[indx_eq2] then begin
       indx_pequinox = indx_eq1
       indx_aequinox = indx_eq2
    endif else begin
       indx_pequinox = indx_eq2
       indx_aequinox = indx_eq1
    endelse
    yr_pequinox = yr[indx_pequinox]
    yr_aequinox = yr[indx_aequinox]
    indx_still = round(interpol(float(indx_yr), yr, yr_still))
    wset,2

    i_time = n_time - n_time_per_period + indx_still
                                ; orbital plot
   i_time0 = n_time - n_time_per_period
   i_time2 = i_time - i_time0  ; time relative to last period

   erase

   ;========================================================
   ; PLOT 1 - the ORBIT
   ;========================================================
   ; general
   dlat = lat[1]-lat[0]
   x = -r * cos(ang)
   y = -r * sin(ang)
   plot, r,ang+!pi, pos=[40./float(xs),10./float(ys),xm/float(xs),ym/float(ys)], /iso, /polar, $
         xs=7, xr=minmax(x)+[-7,7], ys=7, yr=minmax(y)+[-7,7]
   n_time_per_section = n_time_per_period / 12.
   for section = 0, 11,1 do begin
      if section mod 2 then co='808080'xl else co='e0e0e0'xl
      r0 = section*n_time_per_section
      r1 = r0+n_time_per_section-1
      if section eq 11 then r2 = 0 else r2 = r1+1
      polyfill, [0,x[r0:r1],x[r2],0], [0,y[r0:r1],y[r2],0], co=co
   endfor
   oplot, r, ang+!pi, thi=2, /pol
   c = cos(indgen(361)*deg) & s = sin(indgen(361)*deg)
   ; get the angle for the axis
   x0 = x[indx_pequinox]
   y0 = y[indx_pequinox]
   dx_axis = [-6,6] * y0/r[indx_pequinox]
   dy_axis = [6,-6] * x0/r[indx_pequinox]
   for section = 0, 11,1 do begin
      x0 = x[section*n_time_per_section]
      y0 = y[section*n_time_per_section]
      xyouts, (r[section*n_time_per_section]+6+4*ch_size(/x))*cos(ang[section*n_time_per_section]+!pi), $
              (r[section*n_time_per_section]+6+2*ch_size(/x))*sin(ang[section*n_time_per_section]+!pi)-ch_size(/y)/2, $
              string(round(yr[section*n_time_per_section + i_time0]), for='(I4)'), $
              al=0.5
      oplot, x0+dx_axis, y0+dy_axis, th=4
      is_volatile = res.isv[*,section*n_time_per_section + i_time0]
      n_loc = n_elements(is_volatile)
      for i_loc = 0, n_loc-1 do begin
         if is_volatile[i_loc] then color = 'CCCCCC'xl else color='333333'xl
         polyfill, x0+5*cos(nlat[i_loc]+[-1,-1,1,1]*0.5*dlat)*[-1,1,1,-1], $
                   y0 + 5*sin(nlat[i_loc]+[-1,-1,1,1]*0.5*dlat), co=color
      endfor
      oplot, x0+c*5, y0+s*5
   endfor
   ; ORBIT - specific to this time
   x0 = x[i_time2]
   y0 = y[i_time2]
   oplot, [0,x[i_time2]], [0, y[i_time2]], th=5, co=red
   oplot, x0+dx_axis, y0+dy_axis, th=4, co=red
   is_volatile = res.isv[*,i_time2 + i_time0]
   for i_loc = 0, n_loc-1 do begin
      if is_volatile[i_loc] then color = 'CCCCCC'xl else color='333333'xl
      polyfill, x0+5*cos(nlat[i_loc]+[-1,-1,1,1]*0.5*dlat)*[-1,1,1,-1], $
                y0 + 5*sin(nlat[i_loc]+[-1,-1,1,1]*0.5*dlat), co=color
   endfor
   oplot, x0+c*5, y0+s*5, co=red

   ;========================================================
   ; PLOT 2 - the PROJECTION
   ;========================================================
   nlat0 = -res.lat0[i_time]
   is_volatile = res.isv[*,i_time]
   indx_volatile = where(is_volatile, nindx_volatile, compl=indx_not_volatile, ncompl=nindx_not_volatile)
   indxcomp = 2-is_volatile ; is_volatile=1 -> indx=1, is_volatile=0 -> indx=2
   albedo = res.albedo_comp[indxcomp]
   emis = res.emis_comp[indxcomp]
   name_loc = 'latband_60'
   name_loc2map = name_loc+'->map_60_30'
   albedo_map = loc2map(albedo, name_loc2map=name_loc2map)
   ;albedo_map = (rebin(albedo_map,360,180, /sample)*255)
   albedo_map = albedo_map*255
   scale = min(res.r)/res.r[i_time]

   pos_30 = [40./float(xs),ym/float(ys),xm/float(xs),1.-80./ys]
   str = "A!dV!n="+string(res.albedo_comp[1],fo='(F4.2)')
   str += ", !Me!dV!n="+string(res.emis_comp[1],fo='(F4.2)')
   str += ", A!dS!n="+string(res.albedo_comp[2],fo='(F4.2)')
   str += ", !Me!dS!n="+string(res.emis_comp[2],fo='(F4.2)')
   str += ", !4C!X="+string(res.therminertia[0],fo='(E8.2)')+' erg/cm!u2!n s!u1/2!n K'
   str += ", N!d2!n="+string(res.mvtot,fo='(F5.1)')+' g/cm!u2!n'
   xyouts, 0.5, pos_30[3]+40./ys, str, al=0.5, chars=2.5, /norm
   pos_x0 = (pos_30[0] + pos_30[2])/2
   pos_y0 = (pos_30[1] + pos_30[3])/2
   pos_dx = (pos_30[2] - pos_30[0])/2
   pos_dy = (pos_30[3] - pos_30[1])/2
   pos = [pos_x0-scale*pos_dx, pos_y0-scale*pos_dy,pos_x0+scale*pos_dx, pos_y0+scale*pos_dy]

   map_set, nlat0 / deg, 180 * deg, 0., /GRID, /ORTHOGRAPHIC, /ISOTROPIC,$
            pos=pos, /noer, /nobor           
   map_res = map_image(rotate(albedo_map,2), xstart, ystart, missing=white)
   tv, map_res, xstart, ystart
   MAP_GRID, latdel=30, LABEL=0, /HORIZON, co=red, glinethick=2

   ;========================================================
   ; PLOT 3 - the ALBEDO and PRESSURE
   ;========================================================
   plotsym, 0, fill=0
   a = res.a[n_time - n_time_per_period:n_time-1]
   tv = res.tv[n_time - n_time_per_period:n_time-1]
   p = n2vp(tv)
   pos = [xm/float(xs)+170./float(xs), 40./float(ys), 1-40./float(xs),ym/float(ys)-80./ys]
   posa = [pos[0], pos[1], pos[2], (pos[1]+pos[3])/2]
   posb = [pos[0], (pos[1]+pos[3])/2, pos[2], pos[3]]
   yearr = [1865,2113]
   plot, yr, a, pos=posa, /noe, $
         xtitl = 'Year', xtickint = 50, xr=yearr, xs=3,$
         ytitl = 'Albedo', ytickint = 0.2, yr=[0.0, 1.0], ytickname=[' ','','','','',' ']
   oplot, yr_pequinox*[1,1], !y.crange,li=2
   oplot, yr_aequinox*[1,1], !y.crange,li=2
   oplot, yr_nsolstice*[1,1], !y.crange,li=2
   oplot, yr_ssolstice*[1,1], !y.crange,li=2
   xyouts, yr_pequinox,!y.crange[1], 'Equinox', ori=90, al=1
   xyouts, yr_aequinox,!y.crange[1], 'Equinox', ori=90, al=1
   xyouts, yr_nsolstice,!y.crange[1], 'Solstice', ori=90, al=1
   xyouts, yr_ssolstice,!y.crange[1], 'Solstice', ori=90, al=1
   oplot, yr[[i_time2]], a[[i_time2]], ps=8, co=red,syms=2

   yrange = 10.^[floor(min(alog10(p))/3.)*3, ceil(max(alog10(p))/3.)*3]
   ytickv = 10.^(arrgen(alog10(yrange[0]), alog10(yrange[1]), 3) )
   yticks = n_elements(ytickv)-1
   plot, yr, p, pos=posb, /noe, $
         xtickint = 50, xtickna=[' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '], xr=yearr,xs=3,$
         ytitl = 'Pressure (!4l!Xbar)', yticks=yticks, ytickv= ytickv,/ylo
   oplot, yr_pequinox*[1,1],yrange,li=2
   oplot, yr_aequinox*[1,1],yrange,li=2
   oplot, yr_nsolstice*[1,1],yrange,li=2
   oplot, yr_ssolstice*[1,1],yrange,li=2
   oplot, yr[[i_time2]], p[[i_time2]], ps=8, co=red,syms=2
   if keyword_set(showrise) then begin
      oplot, [1988.5, 2006.5], interpol(p, yr, [1988.5, 2006.5]), ps=-8, co=orange,syms=1, thick=2
   endif
   xyouts, yr[i_time2], yrange[1], strtrim(string(p[i_time2], for='(E9.1)'),2) + ' !4l!Xbar', al=(yr[i_time2]-yr[0])/(max(yr)-yr[0])

   ;========================================================
   ; PLOT 4 - the TEMPERATURE AND MASS OF SLAB
   ;========================================================
   ts = res.ts[*,i_time]
   tv = res.tv[i_time]
   ms = res.ms[*,i_time]
   mach = res.v[*,i_time] / res.vs[*,i_time]
   indx_volatile = where(res.isv[*,i_time], nindx_volatile, compl=indx_not_volatile, ncompl=nindx_not_volatile)
   pos = [xm/float(xs)+170./float(xs), (ym+40)/float(ys), 1-40./float(xs), 1-90./float(ys)]
   xtitl = ''
   plot, nlat/deg, ts, pos=pos, /noer, ps=10, $
         title='Temperature and Volatile Mass for '+string(yr[i_time2],for='(F6.1)'), $
         xr=[-90,90], /xs,xticks=6,xtitl=xtitl, $
         yr=[0,60],yticks=6, ytitl='Surface Temperature'
   ms_plot_scale = (min(tv)-10)/max(ms)
   lat_delta = nlat[1]-nlat[0]
   for i_indx = 0, indexcount(indx_volatile)-1 do begin
      indx = indx_volatile[i_indx]
      polyfill, nlat[indx]/deg + [-.5,.5,.5,-.5]*lat_delta/deg, ms[indx]*[0,0,-1,-1]* ms_plot_scale + tv, co=blue
   endfor
   oplot, nlat/deg, ts, ps=10
   ms_max = max(ms, imax)
   if nlat[imax] lt 0 then al=0 else al=1.
   xyouts, nlat[imax]/deg, tv + 0.5*ch_size(/y), string(tv,for='(F5.1)')+' K', al=al
   xyouts, nlat[imax]/deg, tv - ms_max*ms_plot_scale-2*ch_size(/y), string(ms_max,for='(F5.1)')+' g/cm!u2!n max', al=al
   AXIS, yr=[0,60], yax=0, yticks=6
   nlat0 = -res.lat0[i_time]
   nlat0str = string(round(nlat0/deg), for='(I3)')
   austr = string(res.r[i_time],for='(F4.1)')
   oplot, nlat0/deg+[0,0], !y.crange,li=2
   xyouts, nlat0/deg, 60-2.5*ch_size(/y), '  !4k!X'+sunsymbol()+'='+nlat0str+'!uo!n, R='+austr, al = 1- (nlat0 lt 0)

   ; plot the velocity
   vsign = sgn(mach)
   width = linscl(alog10(1e-4 > abs(mach)  < 100), -2.,1.,0.1,1.,/clip)
   indx_south = where(vsign le 0, nsouth)
   indx_north = where(vsign ge 0, nnorth)
   if nnorth eq res.n_lat_edge then nsouth = 0
   if nsouth eq res.n_lat_edge then nnorth = 0
   if nsouth gt 0 then begin
      for ii=0,nsouth-1 do begin
         i = indx_south[ii]
         x0 = nlat[i]/deg
         x1 = nlat[i+1]/deg
         y0 = 60-4.5*ch_size(/y)-width[i]*ch_size(/y)
         y1 = 60-4.5*ch_size(/y)+width[i]*ch_size(/y)
         if mach[i] gt 1. then color='ff'xl else color=byte(255*(1-width[i]))*'000101'xl
         polyfill, [x0,x1,x1,x0,x0], [y0,y0,y1,y1,y0], co=color
         plots, [x1,x0], [y0,y0]
         if ii eq nsouth-1 then plots, [x1,x1], [y0-ch_size(/y),y1+ch_size(/y)]
         if ii eq 0 then plots, [x0,x0, x0+ch_size(/x), x0,x0], [y0,y0-ch_size(/y), (y0+y1)/2, y1+ch_size(/y),y1]
         plots, [x0,x1], [y1,y1]
      endfor
   endif
   if nnorth gt 0 then begin
      for ii=0,nnorth-1 do begin
         i = indx_north[ii]
         x0 = nlat[i]/deg
         x1 = nlat[i+1]/deg
         y0 = 60-4.5*ch_size(/y)-width[i]*ch_size(/y)
         y1 = 60-4.5*ch_size(/y)+width[i]*ch_size(/y)
         if mach[i] gt 1. then color='ff'xl else color=byte(255*(1-width[i]))*'010001'xl
         polyfill, [x0,x1,x1,x0,x0], [y0,y0,y1,y1,y0], co=color
         ;polyfill, [x0,x1,x1,x0,x0], [y0,y0,y1,y1,y0], co=vcolor[i]
         plots, [x1,x0], [y0,y0]
         if ii eq 0 then plots, [x0,x0], [y0-ch_size(/y),y1+ch_size(/y)]
         if ii eq nnorth-1 then plots, [x1,x1, x1-ch_size(/x), x1,x1], [y0,y0-ch_size(/y), (y0+y1)/2, y1+ch_size(/y),y1]
         plots, [x0,x1], [y1,y1]
      endfor
   endif

if keyword_set(saveplot) then begin
   tv2im, run+'_'+string(i_time2,fo='(I03)')+'.png', /co, /tr ; ##########
   if keyword_set(openplot) then begin
      spawn, 'open ' + run+'_'+string(i_time2,fo='(I03)')+'.png'
   endif
endif




end