;Color bar added. 4/24/2008  Hyomin Kim
;Modified to incorporate UNH web service. 4/24/2008  Hyomin Kim
;A code to generate survey plots from the UNH ULF data which are saved in MUD format. 4/13/2007  Hyojin Kim
;
;Modified:
;May 30, 2018 MWC now use ulf_tick_labels to make labels for X axis
;May  8, 2019 MWC fixed time range errors
;Oct 16, 2019 MWC added 'new' Halley, with 20 Hz sampling

;example: ulf_make_plot_cdf, '20171104', '000000', '240000', plot_dir, 'ps', 'NAL', 1, 2048, 200, 0.2, 0, 1
pro ulf_make_plot_cdf, $
   date, $
   start_time, $
   end_time, $
   plot_dir, $
   file_format, $
   station_in, $
   maxfreq, $
   nfft, $
   step, $
   max_B_field, $
   f_low, $
   f_high, $
   upper, $
   lower


;;====== The followings will be located in batchprocess.pro
plot_dir = '/mirl/ULF/data_plot/' + station_in
;file_format   = 'ps'
;maxfreq       = 1
;nfft	      = 2^11
;step          = 300   ;sec
;max_B_field   = 0.2   ;nT/sec
f_low         = 0             ;Hz
f_high        = maxfreq   ;Hz


;====== Specify start/end date and time (this section could be in 'batch_process')
;date        = '20141101'       ;yyyymmdd format
;start_time  = '000000'         ;hhmmss format
;end_time    = '240000'         ;hhmmss format


;====== Create station names for file name
if station_in eq 'SPA' then begin
   station_code_one_letter = 'A'
endif
if station_in eq 'MCM' then begin
   station_code_one_letter = 'M'
endif
if station_in eq 'NAL' then begin
   station_code_one_letter = 'N'
endif
if station_in eq 'LYR' then begin
   station_code_one_letter = 'L'
endif
if station_in eq 'HOR' then begin
   station_code_one_letter = 'H'
endif
if station_in eq 'SDY' then begin
   station_code_one_letter = 'S'
endif
if station_in eq 'POK' then begin
   station_code_one_letter = 'P'
endif
if station_in eq 'IQA' then begin
   station_code_one_letter = 'I'
endif
if station_in eq 'ISR' then begin
   station_code_one_letter = 'R'
endif
if station_in eq 'GAK' then begin
   station_code_one_letter = 'G'
endif
if station_in eq 'SNK' then begin
   station_code_one_letter = 'Y'
endif
if station_in eq 'PG2' then begin
   station_code_one_letter = '2'
endif
if station_in eq 'PG3' then begin
   station_code_one_letter = '3'
endif
if station_in eq 'PG4' then begin
   station_code_one_letter = '4'
endif
if station_in eq 'PG5' then begin
   station_code_one_letter = '5'
endif
if station_in eq 'JBS' then begin
   station_code_one_letter = 'J'
endif
if station_in eq 'KSS' then begin
   station_code_one_letter = 'K'
endif
if station_in eq 'VNA' then begin
   station_code_one_letter = 'U'
endif
if station_in eq 'HAL' then begin
   station_code_one_letter = 'B'
endif



;====== Extracting the date and time
yyyy         = strmid(date, 0, 4)
mm        	 = strmid(date, 4, 2)
dd           = strmid(date, 6, 2)
start_hour   = strmid(start_time, 0, 2)
start_min    = strmid(start_time, 2, 2)
start_second = strmid(start_time, 4, 2)

end_time_cdf = end_time
if end_time eq '240000' then end_time_cdf = '235959'
end_hour_cdf     = strmid(end_time, 0, 2)
end_min_cdf      = strmid(end_time, 2, 2)
end_second_cdf   = strmid(end_time, 4, 2)

end_hour     = strmid(end_time, 0, 2)
end_min      = strmid(end_time, 2, 2)
end_second   = strmid(end_time, 4, 2)
spm          = 60.0
mph          = 60.0
hpd          = 24.0
sph          = 3600.0
start_time_sec = start_hour*sph+start_min*mph+start_second  ;start time in seconds
end_time_sec   = end_hour*sph+end_min*mph+end_second       ;end time in seconds
duration_sec   = end_time_sec - start_time_sec
duration       = duration_sec/sph   ;in decimal hour.

month_array  = ["Jan","Feb","Mar","Apr","May","Jun",$
                "Jul","Aug","Sep","Oct","Nov","Dec"]
month = month_array[fix(mm)-1]
start_epoch = cdf_parse_tt2000(dd + '-' + month + '-' + yyyy + ' ' + $
                               start_hour + ':' + $
                               start_min + ':' + $
                               start_second + '.000')
end_epoch = cdf_parse_tt2000(dd + '-' + month + '-' + yyyy + ' ' + $
                             end_hour_cdf + ':' + $
                             end_min_cdf + ':' + $
                             end_second_cdf + '.000')
if end_time eq '240000' then begin
end_epoch = cdf_parse_tt2000(dd + '-' + month + '-' + yyyy + ' ' + $
                             end_hour_cdf + ':' + $
                             end_min_cdf + ':' + $
                             end_second_cdf + '.999')
endif

;====== Set up X axis (time scale)
ulf_tick_labels, double(start_time_sec), double(end_time_sec), $
                 xtick_count, tick_values, xtick_labels

;====== Initialize data variables
sample_rate     = 10.;Hz
if station_in eq 'KSS' or $
   station_in eq 'VNA' then begin
   sample_rate     = 20.;Hz
endif

period          = 1./sample_rate ;sec
datacount       = long(duration * sph * sample_rate);double(float(duration * sph * sample_rate))
;times           = findgen(datacount) * period/sph

;====== Set up parameters for FFT process
nfft         = nfft             ;number of points per fft: 1024 -> 0.47 mHz
step         = step  ;step size between fft's (secs) (YOU CAN SPECIFY THE STEP SIZE FOR RESOLUTION, HYOJIN)
missing      = double(9999999.0)  ;missing data value
yminor_count = 5
;system_gain  = 4.43               ;system gain = 4.43 V/(nT*Hz)--> applies only to UNH ULF system
max_B_field  = max_B_field               ;nT*Hz (nT/s)
maxfreq      = maxfreq                 ;highest frequency plotted, in Hz
;pstart       = 0 / sph         ;---> fix this (Hyomin 4/1/08)
pstart       = start_time_sec / sph
nsum_value   = duration * 5       ;for reduction of data size of time series plot
f_nyquist    = sample_rate/2     ;Hz
f_low        = f_low             ;Hz
f_high       = f_high              ;Hz
coeff       = digital_filter(f_low/f_nyquist, f_high/f_nyquist, 100, 50)
;Determine spectral array size and number of arrays for spectral images
skip          = round(step / period)
ncols         = round((duration * spm * mph) / period)
if ((ncols mod skip) ne 0) then print,'Error in step size!!!!!!!'
ncols         = ncols/skip
hzperbin      = 1/ (nfft * period)
nrows         = round(maxfreq / hzperbin)
maxfreq_calc  = hzperbin * nrows
spectra       = dblarr(ncols,nrows)

;====== Read in data CDF file
cdf_file = '/mirl/ULF/cdf/' + station_in + '/' + $
           yyyy + '/' + mm + '/mica_ulf_' + $
           strlowcase(station_in) + '_' + date + '_v00.cdf'

file_found = file_test(cdf_file)
if file_found eq 0 then return

cdf_id = cdf_open(cdf_file)

cdf_control, cdf_id, var = 'dbdt_x', get_var_info = rsurvey
nrec_in = rsurvey.maxrec+1

cdf_varget, cdf_id, 'Epoch', epoch, $
  rec_start = 0L, rec_count = nrec_in
epoch = reform(epoch)
cdf_varget, cdf_id, 'dbdt_x', file_array_bx, $
  rec_start = 0L, rec_count = nrec_in
file_array_bx = reform(file_array_bx)
cdf_varget, cdf_id, 'dbdt_y', file_array_by, $
  rec_start = 0L, rec_count = nrec_in
file_array_by = reform(file_array_by)

w_range = where(epoch ge start_epoch and $
                epoch lt end_epoch, datacount)
if datacount lt 1 then begin
   print, 'No data found?!?!?!'
   stop
endif

times = findgen(datacount) * period/sph
plot_array_bx = convol(file_array_bx[w_range], coeff) ;digital filtering
plot_array_by = convol(file_array_by[w_range], coeff) ;digital filtering

;====== Specify graphic format
if file_format eq 'x' then begin
    graphic_type = '.jpg'
	set_plot,'X'
endif else begin
    graphic_type = '.ps'
	set_plot,'PS'
endelse

;====== Specify directory and filename where plots will be saved
file_subdir_plot = file_subdir_plot
file_name_plot   = date + '_' + start_time + '_' + end_time +'_' + $
                   station_code_one_letter + graphic_type
file_path_plot   = plot_dir + '/' + file_name_plot
file_path_base   = plot_dir + '/' + $
                   date + '_' + start_time + '_' + end_time +'_' + $
                   station_code_one_letter

;====== Setup graphic device and plot
loadct, 39
!p.background=255
!p.color=0

if !d.name eq 'PS' then begin
   device, filename=file_path_plot, bits=8, /color, /inches, $
           xoffset=0., yoffset=0., xsize=8.5, ysize=11;was .75, .5, 7 and 9.5 MWC 20170615
   print, 'Device:', !d.name
endif
if !d.name eq 'X' then begin
   device, true = 24, decomposed = 0, retain = 2
   window, 0, retain=2, xsize=640, ysize=828, title=title
   print, 'Device:', !d.name
endif


;====== Set up the plotting window dimensions
leftedge  = 0.05 ;bugger
leftedge = 0.13
rightedge = 0.87
topedge   = 1.0 ;bugger
topedge   = 0.9
dataysize  = 0.16                    ;vertical size of the line plots
gapsize    = 0.04                   ;vertical size of the gap between some panels
xwinsize   = rightedge - leftedge    ;horizontal panel size
x0  = leftedge
x2  = leftedge
x4  = leftedge
x6  = leftedge
;x8  = leftedge
;x10 = leftedge
x1  = rightedge
x3  = rightedge
x5  = rightedge
x7  = rightedge
;x9  = rightedge
;x11 = rightedge
y0  = topedge - dataysize
y2  = y0 - dataysize - gapsize
y4  = y2 - dataysize - gapsize
y6  = y4 - dataysize - gapsize
;y8  = y6 - dataysize - gapsize
;y10 = y8 - dataysize - gapsize
y1  = y0 + dataysize
y3  = y1 - dataysize - gapsize
y5  = y3 - dataysize - gapsize
y7  = y5 - dataysize - gapsize
;y9  = y7 - dataysize - gapsize
;y11 = y9 - dataysize - gapsize


;====== Set up plot positions
position1 = [x0,y0,x1,y1]
position2 = [x2,y2,x3,y3]
position3 = [x4,y4,x5,y5]
position4 = [x6,y6,x7,y7]




doy   = 100 ; fixed (temp)
winsize = convert_coord(xwinsize, dataysize, /normal, /to_device) ;Convert coordinates for the window size from normal to device coords



;print,'Step size = ',(duration/ncols) * spm * mph,' sec
;print,'nrows = ', nrows, ' ,  ncols = ', ncols
;print,'f_max = ', maxfreq_calc,' Hz'


;====== Setup Y axis labels for FFT
if fix(maxfreq) gt 1 and fix(maxfreq) le 5 then begin
    freq_ytickcount = maxfreq
    freq_yticklabel = strarr(freq_ytickcount + 1)
    freq_yticksize = maxfreq / freq_ytickcount
    freq_yticklabel = strarr(freq_ytickcount + 1)
    for y = 0, freq_ytickcount do begin
        freq_yticklabel(y) = strtrim(string(format='(F7.1)', $
                                            freq_yticksize * y), 1)
    endfor
endif

if maxfreq eq 1.000 then begin
    freq_ytickcount = 5.0
    freq_yticklabel = strarr(freq_ytickcount + 1)
    freq_yticksize = maxfreq / freq_ytickcount
    freq_yticklabel = strarr(freq_ytickcount + 1)
    for y = 0, freq_ytickcount do begin
        freq_yticklabel(y) = strtrim(string(format='(F7.1)', $
                                            freq_yticksize * y), 1)
    endfor
endif

if maxfreq gt 0.1 and maxfreq lt 1.0 then begin
    freq_ytickcount = maxfreq*10
    freq_yticklabel = strarr(freq_ytickcount + 1)
    freq_yticksize = maxfreq / freq_ytickcount
    freq_yticklabel = strarr(freq_ytickcount + 1)
    for y = 0, freq_ytickcount do begin
        freq_yticklabel(y) = strtrim(string(format='(F7.1)', $
                                            freq_yticksize * y), 1)
    endfor
endif

if maxfreq eq 0.1 then begin
    freq_ytickcount = 5.0
    freq_yticklabel = strarr(freq_ytickcount + 1)
    freq_yticksize = maxfreq / freq_ytickcount
    freq_yticklabel = strarr(freq_ytickcount + 1)
    for y = 0, freq_ytickcount do begin
        freq_yticklabel(y) = strtrim(string(format='(F7.2)', $
                                            freq_yticksize * y), 1)
    endfor
endif


;====== Set up plot options
charsize1 = 0.8
charsize2 = 1


;====== Plot graphs (time-series and spectrogram; one stations (4 panels) per page)
page_index = 1   ;determine number of page to be plotted in PS format. Prints data from up to 3 stations in one page.
spec_array = dblarr(2, datacount)

;====== Plot time-series (bx)
plot, times, plot_array_bx, nsum = nsum_value, pos=position1, $
    xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
    xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
    ystyle=1, yticklen=-0.01, ycharsize=charsize1, $
    yrange=[-max_B_field, max_B_field], $
    ytitle = station_in + ' bx dB/dt (nT/sec)'

;====== Plot FFT spectrogram (bx)
plot, [0,1], [0,1], /noerase, /nodata, pos=position2,$
    xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
    xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
    yticks=freq_ytickcount, yminor= yminor_count,  ytickname=freq_yticklabel, $
    ystyle=1, yticklen=-0.01, ycharsize=charsize1,$
    ytitle= station_in + ' bx Freq (Hz)'
spec_array[0, *] = times
spec_array[1, *] = plot_array_bx
calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
greyplot, spectra(*,*), position2, winsize, doy, times, pstart, duration, missing, upper, lower

;====== Plot time-series (by)
plot, times, plot_array_by, nsum = nsum_value, pos=position3, /noerase, $
    xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
    xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
    ystyle=1, yticklen=-0.01, ycharsize=charsize1, $
    yrange=[-max_B_field, max_B_field], $
    ytitle = station_in + ' by dB/dt (nT/sec)'


;====== Plot FFT spectrogram (by)
plot, [0,1], [0,1], /noerase, /nodata, pos=position4, $
    xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
    xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
    yticks=freq_ytickcount, yminor= yminor_count,  ytickname=freq_yticklabel, $
    ystyle=1, yticklen=-0.01, ycharsize=charsize1,$
    ytitle= station_in + ' by Freq (Hz)'
spec_array[0, *] = times
spec_array[1, *] = plot_array_by
calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
greyplot, spectra(*,*), position4, winsize, doy, times, pstart, duration, missing, upper, lower

;======Print main titles
main_title1 = yyyy +'/'+mm+'/'+dd+'  '+start_hour+':'+start_min+'-'+end_hour+':'+end_min+' UT'
xyouts, 0.45, .92, main_title1, /normal, alignment=0.5, charsize=charsize2
main_title2 = 'Page ' + strtrim(string(1),2) + '/' + strtrim(string(page_index),2)
xyouts, 0.82, .92, main_title2, /normal, alignment=0.5, charsize=charsize2

x_title = 'Universal Time (hh:mm)'
xyouts, 0.45, y6-0.05, x_title, /normal, alignment=0.5, charsize=charsize2



;==========================================================================
; Reset graphic device
print, file_path_base

if !d.name eq 'PS' then begin
    device, /close
    cmd_str = '/usr/bin/convert ' + file_path_base + '.ps '+ file_path_base + '.jpg'
    spawn, cmd_str
endif
;if !d.name eq 'WIN' then begin
;	WRITE_JPEG, '/tmp/test1.jpg'; [, /ORDER] [, /PROGRESSIVE] [, QUALITY=value{0 to 100}] [, TRUE={1 | 2 | 3}], tvrd()
;    window, /close
;endif

end