;pro mainprocess, filename, filerootdir, filesubdir
pro mainprocess_combine

;==========================================================================
; data filename, filepath - for one file only

openr, unit, "file_info.txt", /GET_LUN
filesubdir = ''
filename = ''
readf, unit, filesubdir
readf, unit, filename
free_lun, unit

;filename = 'ULF-UNH-20220810_1600'
filerootdir = '/home/mirl-ulf/ULF/Data_Files/'
;filesubdir = ''

filename_dat = filename + '.dat'
filepath_dat = filepath(filename_dat, root_dir=filerootdir, subdir=filesubdir)
filename_ps = filename + '_specgram.ps'
filepath_ps  = filepath(filename_ps,  root_dir=filerootdir, subdir=filesubdir)


;==========================================================================
; Check if complete file, or was split during the day

filename_base = strsplit(filename, '_', /EXTRACT) 
filename_base = filename_base[0] + '*.dat'
other_files = file_search(filepath(filename_base, root_dir=filerootdir, subdir=filesubdir))
num_files = size(other_files)
num_files = num_files[1]

;==========================================================================
; Initialize data variables

datacount= ulong(10 * 3600 * 24)
system_gain = 0.51;
;data_array = fltarr(npans + 1, datacount)
;plot_array = fltarr(npans, datacount)

;==========================================================================
; Read dat file and Re-arrange arrays
npans_actual    = 4                     ;number of panels per page
npans           = npans_actual / 2      ;number of iteration for producing panels
file_array      = readdat(filepath_dat, datestr, timestr, datacount_actual)
index           = long(0)
data_array      = fltarr(npans + 1, n_elements(file_array(0,*)))
openw, unit1, "debug.txt", /GET_LUN
for index = 0, n_elements(file_array(0,*)) - 1 do begin
  data_array(0, index) = file_array(2, index) ;time stamps
  data_array(1, index) = file_array(0, index)/system_gain ;x signal
  data_array(2, index) = file_array(1, index)/system_gain ;y signal
endfor
  printf, unit1, data_array

if (num_files gt 1) then begin

    prev_end_time    = long(data_array(0,-1))
    start_time 	     = long(0)

    for i = 1, num_files-1 do begin
        delvar, file_array
        file_array   = readdat(other_files[i], datestr1, timestr1, datacount_actual1)
        data_array1  = fltarr(npans + 1, n_elements(file_array(0,*)))
        for index = 0, n_elements(file_array(0,*)) - 1 do begin
           data_array1(0, index) = file_array(2, index) ;time stamps
           data_array1(1, index) = file_array(0, index)/system_gain ;x signal
           data_array1(2, index) = file_array(1, index)/system_gain ;y signal
        endfor
	   printf, unit1, data_array1
        start_time 	= data_array1(0,0)
        time_diff  	= start_time - prev_end_time
	print, [n_elements(data_array(0,*)), n_elements(data_array1(0,*))]
	print, [start_time, prev_end_time, time_diff]
        prev_end_time 	= data_array(0,-1) 
     endfor
endif
free_lun, unit1
plot_array      = fltarr(npans, n_elements(data_array(0,*)))
times = data_array(0,*)
stop

;==========================================================================
;Parameter

minfrom        	= round(times(0)/60)
minto          	= round(times(-1)/60)

sampling_rate  	= 10         		;Hz
period         	= 1./ sampling_rate 	;Sec
duration       	= (minto - minfrom) / 60.0	;duration of each panel, in hours
nfft           	= 2^10      		;number of points per fft: 1024 -> 0.47 mHz
step           	= (minto - minfrom) / 12.0  ;step size between fft's (secs) (YOU CAN SPECIFY THE STEP SIZE FOR RESOLUTION, HYOJIN)
;doy           	= TBD         		;days of year -  you can specify the date
npans_actual   	= 4          		;number of panels per page
npans          	= npans_actual / 2 	;number of iteration for producing panels
maxvolt        	= 0.2			;voltage range (10: -10V ~ +10V)
maxfreq        	= 2.0    		;highest frequency plotted, in Hz

missing        	= double(9999999.0) 		;missing data value
spm            	= 60.0
mph            	= 60.0
hpd            	= 24.0
sph            	= 3600.0
pstart         	= minfrom / mph  ;0 / sph

system_gain     = 0.51

;plot titles
name 			= strarr(4)
name(0) 		= '   X Signal (Time-series)'
name(1) 		= '   X Signal (Frequency)  '
name(2) 		= '   Y Signal (Time-series)'
name(3) 		= '   Y Signal (Frequency)  '
title   		= '   UNH'


;==========================================================================
; 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  		= hzperbin * nrows
spectra  		= fltarr(ncols,nrows)
print,'Step size is ',(duration/ncols) * spm * mph,' seconds, ncols = ',ncols
print,'Highest frequency plotted is ',maxfreq,' Hz'


;==========================================================================
; Setup graphic device and plot

!p.background=255
!p.color=0

SET_PLOT, 'PS'

loadct, 39
if !d.name eq 'PS' then begin
    device, filename=filepath_ps, /color, /inches, xoffset=0.75, yoffset=0.5, xsize=7.0, ysize=9.5
    print, 'Device:', !d.name
endif

if !d.name eq 'X' then begin
    window, 0, retain=2, xsize=640, ysize=828, title=title
    print, 'Device:', !d.name
endif

if !d.name eq 'WIN' then begin
    window, 0, retain=2, xsize=1024 ,ysize=768, title=title
    print, 'Device:', !d.name
endif

; Set up for plotting:
contour, [[0,0], [1,1]], /nodata, xstyle=4, ystyle=4

; Set up the plotting window dimensions:
leftedge  = 0.1
rightedge = 0.95
topedge   = 0.96
dataysiz  = 0.15              		;vertical size of the line plots
gapsiz    = 0.08             		;vertical size of the gap between some panels
xwinsiz   = rightedge - leftedge    ;horizontal panel size



;===================================================================
; Calculate spectra

doy = 100 ; fixed

;Convert coordinates for the window size from normal to device coords:
winsize = convert_coord(xwinsiz, dataysiz, /normal, /to_device)
splot = topedge - dataysiz   ;y coord of lower left corner next panel

;Print title (date + time)
yearstr = '20' + strmid(datestr, 4, 2)
monthstr = strmid(datestr, 2, 2)
daystr = strmid(datestr, 0, 2)
hourstr = strmid(timestr, 0, 2)
minstr = strmid(timestr, 2, 2)
secstr = strmid(timestr, 4, 2)
title = title + ' - ' + monthstr + '/' + daystr + '/' + yearstr + ' ' + hourstr + ':' + minstr + ':' + secstr + ' (UTC)'
xyouts, 0.4, 1.1, title, alignment=0.5, charsize=1.6


;Print bottom
datestr = 'Plot generated on ' + systime(/UTC)
xyouts, 0.5, -0.1, datestr, charsize=1


;Setup tick label, etc
xtickcount = 6.0
xticksize = (minto - minfrom) / xtickcount
xticklabel = strarr(xtickcount + 1)
for x = 0, xtickcount do begin
    xticklabel(x) = string(format='(F7.1)', minfrom + xticksize * x)
endfor

freq_ytickcount = 5.0
freq_yticksize = maxfreq / freq_ytickcount
freq_yticklabel = strarr(freq_ytickcount + 1);
for y = 0, freq_ytickcount do begin
    freq_yticklabel(y) = string(format='(F7.1)', freq_yticksize * y)
endfor

timeseries_ytickcount = 4.0
timeseries_yticksize = maxvolt / 2.0
timeseries_yticklabel = strarr(timeseries_ytickcount + 1)
timeseries_yticklabel(2) = '0'
for y = 0, 1 do begin
    timeseries_yticklabel(y) = string(format='(F7.1)', -maxvolt + timeseries_yticksize * y)
    timeseries_yticklabel(y + 3) = string(format='(F7.1)', timeseries_yticksize * (y + 1))
endfor

; Plot spectrogram and time-series
for j=0, npans-1 do begin

    posit=[leftedge, splot, rightedge, dataysiz]

    plot_array(0, *) = data_array(0, *)/sph
    plot_array(1, *) = data_array(j+1, *)

    ; plot frequency (fft)
    calcspec, plot_array, spectra, duration, period, sph, nfft, pstart, missing
    greyplot, spectra(*,*), posit, winsize, doy, times, j, duration, missing

    plot, [0,1], [0,1], /noerase, /nodata, charsize=0.8, $
    pos=[posit(0), posit(1), posit(2), posit(1)+posit(3)],$
    xticks=xtickcount, xminor=10, xstyle=1, xticklen=0.05, xcharsize=0.8, $
    xtitle='Time (minute)',$
    xtickname=xticklabel,$
    yticks=freq_ytickcount, yminor= 5, ystyle=1, yticklen=0.01, ycharsize=0.8,$
    ytitle='Frequency (Hz)',$
    ytickname=freq_yticklabel,$
    title=name(2 * j + 1)

    ; plot time-series
    datafrom = minfrom * spm * sampling_rate
    datato = minto * spm * sampling_rate
    plot, data_array(0,*)/spm, data_array(j+1,*), /noerase, charsize=0.8, $
    pos=[posit(0), posit(1)-posit(3)-gapsiz, posit(2), posit(1)+posit(3)-dataysiz-gapsiz], $
    yrange=[-maxvolt, maxvolt], $
    xrange=[data_array(0,0)/spm, data_array(0,-1)/spm], $
    xticks=xtickcount, xminor=10, xstyle=1, xticklen=0.05, xcharsize=0.8, $
    xtitle='Time (minute)', $
    xtickname=xticklabel, $
    yticks=timeseries_ytickcount, yminor=5, ystyle=1, yticklen=0.01, ycharsize=0.8, $
    ytitle='dB/dt (nT/sec)', $
    ytickname=timeseries_yticklabel, $
    title=name(2 * j)
    
    splot = splot - (dataysiz+gapsiz) * 2

endfor



;==========================================================================
; Reset graphic device

if !d.name eq 'PS' then begin
    device, /close
    print, 'PS device closed'
endif



end