;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; This IDL routine makes plots similar to those that are distributed 
;;; with the SWAP data sets.
;;;
;;; Usage:    swapplots, rtn
;;;
;;; RTN:  this is a structure returned by the routine getswappplotdata.pro
;;;       Here are the expected members of the structure:
;;;
;;;  .coin_Hz    event rate, Hz, two dimensional:  time x emergy
;;;              - this is the data that will be behind the pseudo-color
;;;                image plot at the top of the plot set.  The problem
;;;                is that the .coin_Hz array, though rectangular, is
;;;                evenly sampled in neither of its dimensions
;;;              - The key to this routine is that it resamples the
;;;                coin Hz data over time and along a log scale in energy
;;;
;;;  .energies   energies - this array is used to autoscale the plot area
;;;
;;;  .doys, startDoys  Mid Times?A - one of these will be used to set horiz 
;;;                                  scale and control resampling
;;;
;;;  .thetas, .phis, .rjs:     geometries for line plots
;;; 
;;;  Refer to the IDL documentation to understand these plotting commands.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
pro swapplots, rtn, data1k=data1k

  nSize=1000L                            ;;;  Time resampling dimension

  coinDims = size(rtn.coin_hz,/dim)      ;;; The size of .coin_Hz
  coinWid = coinDims[0]
  hgt = coinDims[1]

  mn=min(rtn.doys,max=mx)                ;;; Range of DOYs
  xRng=[mn,mx]
  dxRng = xRng[1] - xRng[0]
  gen1k = lindgen(nSize)                 ;;; Linear integers along time axis

  ;;; Scale linear integers into DOY values - assigns time to each
  ;;;   point along resampled time dimension

  xDoyVals = xRng[0] + dxrng * (gen1k + 0.5d0) / double(nSize)

  ; DESCRIPTION of VALUE_LOCATE
  ;
  ;   VALUE_LOCATE locates the positions of given values within a
  ;   reference array.  The reference array need not be regularly
  ;   spaced.  This is useful for various searching, sorting and
  ;   interpolation algorithms.

  ;;; horizontal array mapping from resampled array into time array:

  vlocDoys = value_locate(rtn.doys,xDoyVals)     ;;; long[nSize]

  ;;; Make an vector of all ones - the height of the image

  onehgt = make_array(hgt,val=1L)           ;;; LONG[hgt], all ones

  ;;; make vLocDays 2-dimensional i.e. stretch it by hgt

  vlocDoysArr = vlocDoys # onehgt  ;;; LONG[nSize,hgt]
  vlocDoysArr2 = vlocDoysArr + (make_array(nSize,val=coinWid) # lindgen(hgt))

  ;;; DOUBLE[nSize,hgt]: Ratio of delta COIN, Hz to delta VALUE_LOCATEd DOYs


  valToDoysArr = (rtn.coin_hz[vlocDoysArr2+1]-rtn.coin_hz[vlocDoysArr2]) $
                 / (rtn.doys[vlocDoysArr+1]-rtn.doys[vlocDoysArr])

  dDoysArr = (xDoyVals # onehgt) - rtn.doys[vlocDoysArr]

  ;;; Here's the beef, the rest is prologue:
  ;;;   Interpolate i.e. reasample

  data1k = rtn.coin_hz[vlocDoysArr2] ;+ (valToDoysArr * dDoysArr)

  help,valToDoysArr,dDoysArr,vlocDoysArr2

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;;; Data have been resampled, the rist of this routine is about plotting
  ;;;   and annotation.  The TVIMAGE.PRO routine is crucial for putting the
  ;;;   2D color image inside the plot axes.

  !x.range = xRng
  !x.style=1
  !x.charsize=0.001
  !x.margin=[25,3]

  !y.style=17
  !y.charsize=2

  !y.margin=[2,2]

  !p.multi=[0,1,3]   ;;; configure plot area for three plots arranged vertically

  loadct,0,/silent

  en2=[min(rtn.energies,max=mx),mx]
  plot,xRng,en2,/nodata,ytitle='Energy/Q [eV/q]',/ylog   ;;; Plot axes, no data
                                                         ;;; - to set up TVIMAGE

  loadct,33,/silent         ;;; start playing with color tables so pseudocolor
                            ;;;    does not overwrite plotting colors
  tvlct,red,grn,blu,/get    ;;; get color lookup table
  red[0]=0                  ;;; - protect black
  grn[0]=0
  blu[0]=0
  red[255]=255              ;;; - protect white
  grn[255]=255
  blu[255]=255
  tvlct,red,grn,blu         ;;; - reload lookup table with blk & white protected

  ;;; display image using loaded pseudocolor lookup tables

  tvimage,bytscl(alog10(data1k>.1 <150)),/overplot

  loadct,0,/silent               ;;; reload gray lookup table

  !p.multi=[3,1,3]
  plot,xRng,en2,/nodata,ytitle='Energy/Q [eV/q]',/ylog,/noerase ;;;overplot axes

  ;;; change to 6 plots stacked vertical to put 4 plots in remaining 2/3
  ;;; - normal XY line plots after this

  !p.multi=[4,1,6]         
  xs = rtn.doys
  plot,xs,total(rtn.coin_hz,2),/ylog,ytitle='Coin Rate!CEnergy Sum!c[cnts/samp]'

  sxs = rtn.startdoys
  plot,sxs,rtn.thetas,ytitle='Sun Theta!CAngle [deg]'

  plot,sxs,rtn.phis,ytitle='Sun Phi!CAngle [deg]'

  !x.charsize=2
  !y.margin=[4,2]
  plot,sxs,rtn.rjs,ytitle='Jupiter-S/C!CRange [Rj]'

  rjmin = min(rtn.rjs,max=rjmax)
  ypos = 1.15 * rjmin - 0.15 * rjmax
  xyouts,xRng[0],ypos,/data,'DOY=>',align=1.0,charsiz=1.01,/noclip
  
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; This is test code to drive a test case
;;; - jupswap.lst file contains SWAP FITS data

strTyp=size({a:0b},/type)
rtnTyp=size(rtn,/type)

;;; call the GETSWAPPLOTDATA.PRO routine to read the list of 10 files,
;;; covering 10 days, in flat ASCII file jupswap.lst ....

if rtnTyp ne strTyp then rtn=getswapplotdata('jupswap.lst')

;;; ... and plot the results

swapplots, rtn, data1k=data1k

end