; Code written by Chrystal Moser-Gauthier Sept 25, 2023
; will make cdfs of the new DAQs data (grant MICA-2022)

;example: ulf_make_cdfs, '20170321', '000000', '240000', plot_dir, 'NAL', 2048, 3
pro ulf_make_cdfs, $
   date, $
   start_time, $
   end_time, $
   plot_dir, $
   station_in, $
   nfft, $
   axis_count

;==========================================================================
; Station info
if station_in eq 'UNH' then begin
   latitude = '43.1463 N'
   longitude = '70.9447 W'
   station_name = 'UNH>University of New Hampshire'
   system_type = 'UNH'
endif
if station_in eq 'NAL' then begin
   latitude = '78.5540 N'
   longitude = '11.5580 E'
   station_name = 'NAL>Ny-Alesund'
   system_type = 'UNH'
endif
if station_in eq 'LYR' then begin
   latitude = '78.0886 N'
   longitude = '16.0242 E'
   station_name = 'LYR>Longyearbyen'
   system_type = 'UNH'
endif
if station_in eq 'HOR' then begin
   latitude = '77.0035 N'
   longitude = '15.3277 E'
   station_name = 'HOR>Hornsund'
   system_type = 'UNH'
endif
if station_in eq 'ISR' then begin
   latitude = '78.0366 N'
   longitude = '13.3823 E'
   station_name = 'ISR>Isfjord Radio'
   system_type = 'UNH'
endif
if station_in eq 'NAQ' then begin
   latitude = '61.1567 N'
   longitude = '45.4254 W'
   station_name = 'NAQ>Narsarsuaq'
   system_type = 'UNH'
endif
if station_in eq 'IQA' then begin
   latitude = '63.4537 N'
   longitude = '68.3065 W'
   station_name = 'IQA>Iqaluit'
   system_type = 'UNH'
endif
if station_in eq 'SNK' then begin
   latitude = '56.3218 N'
   longitude = '79.1387 W'
   station_name = 'SNK>Sanikiluaq'
   system_type = 'UNH'
endif

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

filename     = 'ULF-' + station_in + '-' + date + '_0000'
filerootdir  = '/home/mirl-ulf/ULF/Data_Files/' + station_in + '_incoming/'


filename_dat = filename + '.dat'
filepath_dat = filepath(filename_dat, root_dir=filerootdir)

if axis_count eq 2 then begin
    arr_ind = 0
    sample_rate = 10.0
endif
if axis_count eq 3 then begin
    arr_ind = 2
    sample_rate = 6.6666
    if station_in eq "NAL" then begin
        sample_rate = 10.0
    endif
endif

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

datacount   = ulong(sample_rate * 3600 * 24)
system_gain = float(150.0 * 121.0 * 31.0 /1000000.0) ;150uV/(nT Hz) * preamp_gain * bufferboard gain
npans       = 2

;==========================================================================
; Read dat file and Re-arrange arrays

file_array      = readdat(filepath_dat, datestr, timestr, datacount_actual, axis_count)
index           = long(0)
data_array      = fltarr(npans + 1, datacount)
plot_array      = fltarr(npans, datacount)
for index = 0, n_elements(file_array(0,*)) - 1 do begin
  data_array(0, index) = index/sample_rate                       ;time stamps
  data_array(1, index) = file_array(arr_ind, index)/system_gain  ;x signal
  data_array(2, index) = file_array(1, index)/system_gain        ;y signal
;  data_array(3, index) = file_array(2, index)/system_gain       ;z signal
endfor

times    = data_array(0,*)
array_bx = data_array(1,*)
array_by = data_array(2,*)

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

minfrom         = 0.0
minto           = 14400.0

period          = 1./ sample_rate     ;Sec
duration        = (minto - minfrom) / 60.0      ;duration of each panel, in hours
step            = (minto - minfrom) / 12.0  ;step size between fft's (secs) (YOU CAN SPECIFY THE STEP SIZE FOR RESOLUTION, HYOJIN)
maxvolt         = mean(data_array(1,*))+0.2;voltage range (10: -10V ~ +10V) 
maxfreq         = 1.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
nsum_value      = duration*5

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)


year_str        = strmid(date, 0, 4)
year_str_2      = strmid(date, 2, 2)
year_int        = fix(year_str)
month_array     = ["Jan","Feb","Mar","Apr","May","Jun",$
                  "Jul","Aug","Sep","Oct","Nov","Dec"]
month           = strmid(date, 4, 2)
month_str       = month_array[fix(strmid(date, 4, 2))-1]
day_str         = strmid(date, 6, 2)

start_epoch     = day_str + '-' + month_str + '-' + year_str + ' 00:00:00.000'
start_time_long64 = cdf_parse_tt2000(start_epoch)

db_times        = times/sph
spec_array      = dblarr(2, datacount)

;==========================================================================
; Create CDF file and assign attributes


skeleton_base = ' -delete ./'
new_harddrive = '/media/mirl-ulf/4TB_HDD1'
cdf_file      = new_harddrive + '/mirl/ULF/cdf/' + station_in + $
                '/' + year_str + '/' + month + '/mica_ulf_' + $
                strlowcase(station_in) + '_' + date + '_v00'
skt_name      = 'mica_ulf'

cmd_base      = '/usr/local/cdf/bin/skeletoncdf -cdf ' + cdf_file
cmd           = cmd_base + skeleton_base + skt_name
spawn, cmd
spawn, 'chmod 644 ' + cdf_file + '.cdf'
cdf_id        = cdf_open(cdf_file)

cdf_attput, cdf_id, 'Longitude', 0, longitude
cdf_attput, cdf_id, 'Latitude', 0, latitude
cdf_attput, cdf_id, 'Station_name', 0, station_name

cdf_varput, cdf_id, 'Epoch', $
            start_time_long64 + long64(dindgen(datacount) / sample_rate * 1.d9)
cdf_varput, cdf_id, 'dbdt_x', array_bx, /zvariable
cdf_varput, cdf_id, 'dbdt_y', array_by, /zvariable

spec_times = findgen(ncols) * (duration * spm * mph) / ncols
cdf_varput, cdf_id, 'Epoch_spectra', $
            start_time_long64 + long64(spec_times) * 1000000000LL

; Calculate Spectra
spec_array[0, *] = db_times
spec_array[1, *] = array_bx

;var_idx = cdf_varcreate(cdf_id, 'spectra_x_1Hz', /zvariable)
;var_idy = cdf_varcreate(cdf_id, 'spectra_y_1Hz', /zvariable)

calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
cdf_varput, cdf_id, 'spectra_x_1Hz', transpose(spectra), /zvariable

spec_array[1, *] = array_by
calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
cdf_varput, cdf_id, 'spectra_y_1Hz', transpose(spectra), /zvariable

print, spectra

;make 5 Hz spectra
maxfreq       = 5.0
nrows         = round(maxfreq / hzperbin)
spectra       = dblarr(ncols,nrows)

spec_array[1, *] = array_bx
calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
cdf_varput, cdf_id, 'spectra_x_5Hz', transpose(spectra), /zvariable

spec_array[1, *] = array_by
calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
cdf_varput, cdf_id, 'spectra_y_5Hz', transpose(spectra), /zvariable

cdf_close, cdf_id

end