pro praspec
; Creates spectrograms from Voyager PRA IDL save sets. 
;
; This script reads savesets.txt, a list of IDL save sets. 
; This script uses the CDAWlib package of IDL routines.
; These routines are available from:
; http://spdf.gsfc.nasa.gov/CDAWlib.html
;
; IDL save sets were created using the routine convertvoyHIband.pro
; from the Voyager 1 and 2 PRA Highband data residing at The University of Iowa 
;
; Leonard Garcia 2012 February 10 with minor updates 2016 January 5

; Choose: left channel (0) or right channel(1)?
;
lrchan=1

; Need to create a frequency array and a time array to serve as axes
; for the spectrograms.
; 
;
;From Wang and Carr, 1994
;
;The high-band receiver consisted of 128 channels of 200 kHz bandwidth
;each, with center frequencies spaced at 307.2 kHz intervals from 1.2 MHz to 40.4 MHz.
;
flo=1.2
fhi=40.4
fstep=0.3072
farray=fltarr(128)
for i=0,127 do begin
   farray[i]=flo+fstep*(i)
endfor
;print,farray

;read IDL save set and restore
;filename="v2790716.idl"
;filename=DIALOG_PICKFILE(filter="*.idl")
;read, 'Input filename? ',filename

filename=" " & spacecraft = " " & date= " " & scname=" "
openr,lun,"savesets.txt",/get_lun
while not eof(lun) do begin
   readf,lun,filename
   restore,filename
   shortname=strmid(filename,0,8)
   spacecraft = strmid(shortname,0,2)
   case spacecraft of
      'v1': scname="Voyager 1"
      'v2': scname="Voyager 2"
   endcase
   case lrchan of
      0: channel="L"
      1: channel="R"
   endcase
   date=strmid(shortname,2,8)
; The structure looks similar to this.
;IDL> help,pra,/struc
;** Structure PRA_STRUC, 5 tags, length=528, data length=528:
;   DUM1            LONG               520
;   YMD             LONG            790716
;   MSEC            LONG             46695
;   VALUES          INT       Array[128, 2]
;   DUM2            LONG               520
; corresponding to:
;file separation variable
;year, month, day information
;millisecond decimal value of the day
;Integer array (128,2) for 128 left and right channels (NOTE 128 channels for Hi-band; 70 channels for Lo-band)
;file separation variable
;
; This structure is repeated up to 2300 times, that is:
; IDL> help,pra.values
;<Expression>    INT       = Array[128, 2, 2300]
   prasize=size(pra)
   ntimes=prasize[1]
   spec3=pra.values
; Which channel? Left or Right
   spec=reform(spec3[*,lrchan,*])
   
;time array 48 second resolution starting from pra.msec. There are 1800 48 sec
;intervals in one 24 hour period
   tarray=pra.msec
   ymd=pra.ymd
   ymd1=ymd[0]
; Some of the times are zero. Estimate what it should have been
   for i=0,ntimes-1 do begin
      ttest=tarray[i]
      if (ttest eq 0 AND i ne 0) then begin
         tbefore=tarray[i-1]
         ttest=tbefore+48.0*1000.0
         tarray[i]=ttest
      endif
   endfor
;print,tarray
   tarray=tarray/1000.0
   tarray=tarray/60.0
   tarray=tarray/60.0
;if(ntimes gt 1800) then begin
;   tarray=tarray[0:1800]
;endif
; Rotate the array
   spec2=rotate(spec,3)
   minspec = min(spec2)
   maxspec = max(spec2)
   print,"Min value= ",minspec,"     Max value= ",maxspec
;help,spec,tarray,farray,pra.ymd[0]
   plottitle=scname+" PRA Highband "+shortname+" "+channel+"   date(yymmdd)="+date
;
; Output device (Local=0, PS color=1, PS gray=2, PS BW=3, Tek=4, Zbuf=5, GIF=6, PNG=7, ION=9)? 0
; Create GIFS
   deviceopen,6
; Default is rainbow 13, for black and white plots uncomment the following
;   loadct,0
   spectrogram,spec2,tarray,farray,xrange=[0,24],yrange=[1,41],xstyle=1,ystyle=1,title=plottitle,xtitle="HOURS",ytitle="FREQ (MHz)",cscale=[minspec,maxspec],ctitle="millibels"
   deviceclose
   outfile=shortname+channel+".gif"
   cmd="mv idl.gif "+outfile
   print,cmd
   spawn,cmd
endwhile
free_lun,lun
end