; several necessary routines precede the actual routine itself

;-------------------------------------------------------------
	function ymd2dn,yr,m,d, help=hlp
 
	;----  Days before start of each month (non-leap year)  -----
	idays = [0,31,59,90,120,151,181,212,243,273,304,334,366]
 
	;----  Correct for leap year if month ge 3  -------------
	lpyr = (((yr mod 4) eq 0) and ((yr mod 100) ne 0)) $
            or ((yr mod 400) eq 0) and (m ge 3)
 
	dy = d + idays(m-1) + lpyr
	return, dy
 
	end

;-------------------------------------------------------------

	PRO YDN2MD,YR,DY,M,D, help=hlp
 
	; Days before start of each month.
	YDAYS = [0,31,59,90,120,151,181,212,243,273,304,334,366]
 
	; Correct YDAYS for leap-year.
	IF (((YR MOD 4) EQ 0) AND ((YR MOD 100) NE 0)) $
            OR ((YR MOD 400) EQ 0) THEN YDAYS(2) = YDAYS(2:*) + 1
 
	FOR I = 1, 12 DO BEGIN
	  IF DY LE YDAYS(I) THEN GOTO, NEXT
	ENDFOR
 
	PRINT,' ydn2md: error in day number.'
	M = -1
	D = -1
	RETURN
 
NEXT:	M = I
	D = DY - YDAYS(M-1)
	RETURN
 
	END

;----------------------------------------------------------------------------------

function strreplace,instring,inchar,outchar

len=strlen(inchar)
enter:
   apos=strpos(instring,inchar)
   if (apos lt 0) then goto,break
   instring=strmid(instring,0,apos)+outchar+strmid(instring,apos+len,999)
   goto,enter
break:

return,instring
end

;----------------------------------------------------------------------------------
pro gb_bkg_sub,spec,bkg,specsub

 ; determine background spectrum

 IF (n_params(0) LT 2) THEN BEGIN
   print,'Usage: gb_bkg_sub, spec, bkg [, specsub]'
   print,''
   print,'Determines background BKG for a spectrum SPEC using median of values '
   print,'from low times. '
   print,'Optionally returns background-subtracted spectrum SPECSUB.'
   return
 END

 nx=n_elements(spec[*,0])
 nfr=n_elements(spec[0,*])

 ; create background array
 ; identify low times in light curve: needs to be general, no
 ; assumptions about where frequency is, so use low data
 avespec=total(spec,1) ; sum in time
 ii=sort(avespec)      ; sort to order
 lc=total(spec[*,ii[0:nfr/3L]],2) ; use lowest third of frequency points

 ; Optional background subtraction: median of lowest 60 points
 jj=sort(lc) & kk=jj[0:((nx/10L)<59)] &  bkg=reform(spec[0,*])
 for ij=0,nfr-1 do bkg[ij]=median(spec[kk,ij])
 specsub=spec-rebin(reform(bkg,1,nfr),nx,nfr)

end

;------------------------------------------------------------------------------
function secs_to_str,number,decimal=decimal,wrap=wrap

; returns time string from seconds
; if DECIMAL set, returns DECIMAL places after decimal point, otherwise integer
; if /WRAP set, goes from 24 to 00, default is 24 to 25

; decide whether array or single number
sz = size(number)

if (sz[0] gt 0) then begin ; array
   
   secstring=strarr(sz[1])
   
   for j=0,sz[1]-1 do begin

      starts=strarr(3)
      startv=fix(sixty((number[j]+.01)/3600.0))
      if keyword_set(wrap) then startv[0] = startv[0] mod 24
      for i=0,2 do  $
          starts(i)=strtrim(string(startv(i),format='(i2.2)'),2)
      secstring[j]=starts(0)+':'+starts(1)+':'+starts(2)
   
      if keyword_set(decimal) then begin
         fmt="(f"+string(decimal+2,format='(i1)')+"."+string(decimal,format='(i1)')+")"
         dec=strtrim(string(number[j]-long(number[j])+.00001,format=fmt),2) 
         secstring[j]=secstring[j]+strmid(dec,1,decimal+1)
      endif
   end
 
endif else begin    ;    single number

   starts=strarr(3)
   secstring=' '
   startv=fix(sixty((number+.01)/3600.0))
   if keyword_set(wrap) then startv[0] = startv[0] mod 24
   for i=0,2 do  $
       starts(i)=strtrim(string(startv(i),format='(i2.2)'),2)
   secstring=starts(0)+':'+starts(1)+':'+starts(2)
   
   if keyword_set(decimal) then begin
      fmt="(f"+string(decimal+2,format='(i1)')+"."+string(decimal,format='(i1)')+")"
      dec=strtrim(string(number-long(number)+.00001,format=fmt),2) 
      secstring=secstring+strmid(dec,1,2)
   endif
    
endelse

return,secstring

end
;------------------------------------------------------------------------------
function mpow, arr, power

return, (arr > 0.)^power -(-arr>0.)^power

end
;------------------------------------------------------------------------------
;+
; NAME:
;     PLOT_RSTN_SRS
; PURPOSE:
;     Plots RSTN SRS 25-180 MHz dynamic spectrum SRS files, gets data if needed.
; CATEGORY:
;     RSTN SRS
; CALLING SEQUENCE:
;     plot_rstn_srs, date, tstart, tend [, file=FILE, obs=OBS, /BW]
; INPUTS:
;     DATE      String in the format YYYYMMDD
;     TSTART    Starting time string in the format HHMMSS
;     TEND      Ending time string in the format HHMMSS
; OPTIONAL INPUTS:
;     FILE      Name of file already present: if not supplied, wget is used to get a
;                  file from the NGDC ftp site
;     OBS       Observatory: 4-letter string from PALE (Palehua), LEAR (Learmonth),
;                  SVTO (San Vito), HOLL (Holloman)
;                  Typical time ranges SVTO 5-16, HOLL 13-24, PALE 17-04, 
;                  LEAR 23-10, SGMR 11-22
;     /BW       If set, does black and white, else color table 39
; OUTPUTS:
;     Write and displays Postscript file idl.ps 
; RESTRICTIONS: 
;     Uses the UNIX wget command to get data from ftp site
;     Needs READ_RSTN_SRS, and the JHUAPL (for timeaxis, also located in
;     SSW/radio/ethz) and GSFC Astronomy libraries of routines.
; MODIFICATION HISTORY:
;     Written 14-Sep-2004 by Stephen White
;-

pro plot_rstn_srs, date, tstart, tend, file=file, obs=obs, bw=bw

; date is YYYYMMDD, times HHMMSS, wraps if end before start
; file is optional if file already here, else wgets it
; obs specifies observatory HOLL, SVTO, LEAR, PALE if don't want default
; BW does black and white

IF (n_params(0) LT 3) THEN BEGIN
PRINT,'---------------------------------------------------------------------'
PRINT,'Call: plot_rstn_srs, date, tstart, tend, file=file, obs=obs, /bw'
PRINT,'---------------------------------------------------------------------'
PRINT,'Reads data from SRS file, wgets from ftp site for observatory OBS if'
PRINT,'FILE is not supplied. DATE is UT date YYYYMMDD of event, will try to get right file.'
PRINT,'OBS can be SVTO, HOLL, PALE, LEAR, SGMR'
PRINT,'Typical time ranges SVTO 5-16, HOLL 13-24, PALE 17-04, LEAR 23-10, SGMR 11-22'
PRINT,'/BW prints black-and-white reversed'
PRINT,'Uses UNIX wget command and JHUAPL and GSFC IDL libraries.'
PRINT,'---------------------------------------------------------------------'
return
endif

 ; get rid of semicolons in times if they are there
 tstart=strreplace(tstart,':','') & tend=strreplace(tend,':','')
 ta=ten(fix(strmid(tstart,0,2)),fix(strmid(tstart,2,2)),fix(strmid(tstart,4,2)))*3600.
 tz=ten(fix(strmid(tend,0,2)),fix(strmid(tend,2,2)),fix(strmid(tend,4,2)))*3600.
 if (tz lt ta) then tz=tz+86400.0

 ; get observatory
 if not keyword_set(obs) then begin
    if (ta gt 0.0) then obs='LEAR'
    if (ta gt (8.*3600.0)) then obs='SVTO'
    if (ta gt (15.*3600.0)) then obs='HOLL'
    if (ta gt (21.*3600.0)) then obs='PALE'
    if (ta lt (2.*3600.0)) then obs='PALE'
 endif
 ; parse file prefix
 shobs=strmid(obs,0,2)
 if (obs eq 'LEAR') then shobs='LM'
 if (obs eq 'SGMR') then shobs='K7'

 ; need to be careful with dates: PALE data start with times on day preceding
 ; date of file, early LEAR data start on file 1 day later
 doy=ymd2dn(strmid(date,0,4),strmid(date,4,2),strmid(date,6,2))
 if (obs eq 'PALE') then begin
    ; now need to fix times since ta,tz just raw conversions
    if (ta gt 40000.) then doy=doy+1
    if (ta lt 40000.) then ta=ta+86400.
    if (tz lt 40000.) then tz=tz+86400.
 endif
 if (obs eq 'LEAR') then begin
    if (ta gt 40000.) then doy=doy+1
    ; times in LEAR files go above 24 at day boundary
    if (ta lt 40000.) then ta=ta+86400.
    if (tz lt 40000.) then tz=tz+86400.
 endif
 if (doy lt 1) then begin
    print,'Sorry, I am not programmed to handle end-of-year transitions.'
    print,'Terminating ....'
 endif
 yr=strmid(date,0,4)
 ydn2md,yr,doy,mm,dd
 yymm=strmid(yr,2,2)+string(mm,format='(I2.2)')
 yymmdd=strmid(yr,2,2)+string(mm,format='(I2.2)')+string(fix(dd),format='(I2.2)')

 if not keyword_set(file) then begin
    ; parse file name and location
    file=shobs+yymmdd+'.SRS'
    ; check whether the file is already there
    if (file_test('/tmp/'+file)) then begin
       print,'File /tmp/'+file+' already seems to be there.'
       
    endif else begin
       print,'Copying file '+file+' from NGDC. Will take a while.'
       spawn,'wget ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/SPECTRAL_RSTN/'+obs+yymm+'/'+shobs+yymmdd+'.SRS'
       nf=0 & spawn,'\ls '+file+' | wc -l',nf
       ; try .srs extension if no file found
       if (nf lt 1) then spawn,'wget ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/SPECTRAL_RSTN/'+obs+yymm+'/'+shobs+yymmdd+'.srs'
       if (file_test(shobs+yymmdd+'.srs')) then spawn,'mv '+shobs+yymmdd+'.srs '+shobs+yymmdd+'.SRS'
       nf=0 & spawn,'\ls '+file+' | wc -l',nf
       if (nf lt 1) then begin
          print,'There seems to be a problem: no data file! Check web site'
          print,'ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/SPECTRAL_RSTN/'
          print,'for an observatory with suitable data and specify OBS keyword'
          print,'Terminating ....'
          return
       endif
       spawn,'mv '+file+' /tmp'
    endelse
    file='/tmp/'+file
 endif 

 ; now read the file
 read_rstn_srs, file, ts, fr, sp, fdate
 print,'---------------------------------------------------------------------'
 print,''
 print,'Data run from ',secs_to_str(ts[0]),' to ',secs_to_str(ts[n_elements(ts)-1]),' on ',string(fdate)

 ; extract the desired period
 ix=where(ts ge ta) & i0=ix[0]
 ix=where(ts gt tz) & i1=ix[0]-1
 if (i1-1 le i0) then begin
    print,'Oops, something is wrong with the times: data requested from '+tstart+' to '+tend
    return
 end

 if (date ne fdate) then print,'Warning: dates ',date,' and ',fdate,' appear to differ'
 print,''
 print,'---------------------------------------------------------------------'

 tsecs=ts[i0:i1]
 freq=fr
 spec=sp[i0:i1,*]
 ; dont want spectrum to be byte array or bkgsub doesn't work due to wrapping
 gb_bkg_sub,float(spec),bkg,specsub

 ; set up device
 pltdev=!d.name
 set_plot,'PS'
 !p.font=0
 !p.thick=3.0
 !x.thick=3.0
 !y.thick=3.0
 !p.charthick=3.0
 !p.charsize=1.2
 device,/portr,/inch,xs=6.5,ys=6.0,xoff=1.2,yoff=3.5
 device,/color,bits_per_pixel=8
 !x.margin=[0,0] & !y.margin=[0,0] & !p.multi=[0,1,1,0]
 ; get dimensions
 contour,[[0,0],[1,1]],/nodata, xstyle=4, ystyle = 4
 px = !x.window * !d.x_vsize ;Get size of window in device units
 py = !y.window * !d.y_vsize
 swx = px(1)-px(0)       ;Size in x in device units
 swy = py(1)-py(0)       ;Size in Y
 if keyword_set(bw) then loadct,0 else loadct,39

 nx=n_elements(specsub[*,0])
 nfr=n_elements(freq)
 ifr=where((freq gt 30.0) and (freq lt 81.0)) & smed=median(specsub[*,ifr])
 smx=max(specsub[*,ifr]>(smed))
 smn=stdev(specsub[0:60,ifr])
 print,smn,smx
 smx=min([10.*smn,smx])
 bb=bytscl(mpow(specsub,0.5),min=sqrt(smn)/3.,max=sqrt(smx))
 if keyword_set(bw) then bb=255-bb
 ; bb=bytscl((specsub))
 ; the mapping of index to frequency has a break at 400 = 75 MHz
 ; need to fix this using scaling of TV display
 ; first display lower portion to right area
 tv, bb[*,0:400], px(0), py(0), xsize=swx,$
     ysize=swy*(freq[400]-freq[0])/(freq[nfr-1]-freq[0]), /device
 ; now display upper portion to right area
 tv, bb[*,401:nfr-1], $
     px(0), py(0)+swy*(freq[400]-freq[0])/(freq[nfr-1]-freq[0]), xsize=swx,$
     ysize=swy*(freq[nfr-1]-freq[401])/(freq[nfr-1]-freq[0]), /device
 contour,specsub,tsecs,freq,/noerase, xstyle=5, ystyle = 1,color=0,$
    pos = [px(0),py(0), px(0)+swx,py(0)+swy],/device,$
    xthick=3.0,ythick=3.0,ytitle='Frequency (MHz)',levels=[1.e6]
 timeaxis,title='Time (UT) '+date,size=1.2,color=0

 loadct,0
 xyouts,0.50,1.02,'RSTN/'+obs,/norm,size=1.6,align=0.5
 device,/close & spawn,'gsview idl.ps'
 spawn,'mv idl.ps PS.'+date+'_'+obs

end