;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

;example: unh_cdf, '20170321', '000000', '240000', plot_dir, 'NAL', 2048, 200, 0.2
pro unh_cdf, $
   date, $
   start_time, $
   end_time, $
   plot_dir, $
   station_in, $
   nfft, $
   step, $
   max_B_field

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)


;;====== The followings will be located in batchprocess.pro
database_dir = '/media/mirl-ulf/4TB_HDD1/mirl/ULF/database'
;plot_dir = '/mirl/ULF/data_plot/' + station_in
plot_dir = '/home/mirl-ulf/ULF/Data_Files/' + station_in + '_plots/'
maxfreq       = 1
nfft	      = 2^11
step          = 60   ;step size between fft's in seconds
max_B_field   = 0.2   ;nT/sec
has_z = 0

;====== 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

system_gain = 1.

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 'SDY' then begin
   latitude = '66.5916 N'
   longitude = '5.0568 W'
   station_name = 'SDY>Sondrestrom'
   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
if station_in eq 'MCM' then begin
   latitude = '77.85 S'
   longitude = '166.67 W'
   station_name = 'MCM>McMurdo'
   if year_int lt 2015 then begin
      system_type = 'AUG'
   endif else begin
      system_type = 'NJIT'
   endelse
endif
if station_in eq 'SPA' then begin
   latitude = '90.0000 S'
   longitude = '00. W'
   station_name = 'SPA>South Pole'
   if year_int lt 2015 then begin
      system_type = 'AUG'
   endif else begin
      system_type = 'NJIT'
   endelse
   has_z = 1
endif
if station_in eq 'JBS' then begin
   latitude = '76.5 S'
   longitude = '164.2 E'
   station_name = 'JBS>Jang Bogo Station'
   system_type = 'NJIT'
endif
if station_in eq 'KSS' then begin
   latitude = '62.2 S'
   longitude = '57.8 W'
   station_name = 'KSS>King Sejong Station'
   system_type = 'NJIT'
endif
if station_in eq 'VNA' then begin
   latitude = '70.67 S'
   longitude = '8.25 W'
   station_name = 'VNA>Neumayer'
   system_type = 'NJIT'
endif
if station_in eq 'HAL' then begin
   latitude = '75.5 S'
   longitude = '26.6 W'
   station_name = 'HAL>Halley Station'
   system_type = 'NJIT'
   has_z = 1
endif
if station_in eq 'PG2' then begin
   latitude = '84.42 S'
   longitude = '57.96'
   station_name = 'PG2'
   system_type = 'AAL-PIP'
endif
if station_in eq 'PG3' then begin
   latitude = '84.81 S'
   longitude = '37.63'
   station_name = 'PG3'
   system_type = 'AAL-PIP'
endif
if station_in eq 'PG4' then begin
   latitude = '83.34 S'
   longitude = '12.25'
   station_name = 'PG4'
   system_type = 'AAL-PIP'
endif
if station_in eq 'PG5' then begin
   latitude = '81.96 S'
   longitude = '5.71'
   station_name = 'PG5'
   system_type = 'AAL-PIP'
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_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.

;====== Specify MUD sources
;database_path = database_dir + '/' + station_in + '_' + $
;                strtrim(string(yyyy), 1) + '_' + $
;                strtrim(string(mm), 1) + '.MUD'
database_path = '/home/mirl-ulf/ULF/Data_Files/MUD/' + station_in + '_' + $
		strtrim(string(yyyy), 1) + '_' + $
		strtrim(string(mm), 1) + '.MUD'

;====== Initialize data variables
sample_rate     = 10.    ;Hz
period          = 1./sample_rate   ;sec
datacount       = long(duration * sph * sample_rate);double(float(duration * sph * sample_rate))
db_times        = dindgen(datacount) * period/sph

if station_in eq 'VNA' or $
   station_in eq 'KSS' or $
   (station_in eq 'HAL' and year_int gt 2018) then begin

   nfft = 2^12
   sample_rate = 20. ;Hz
   period =  1. / sample_rate  ;sec
   datacount = long(duration * sph * sample_rate)
   db_times = dindgen(datacount) * period/sph

endif

;====== Set up parameters for FFT process
nfft         = nfft             ;number of points per fft: 1024 -> 0.47 mHz
missing      = double(9999999.0)  ;missing data value
yminor_count = 5
;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
;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)
spectra       = dblarr(ncols,nrows)


;
;====== Read in data (stored in various formats)
array_bx = dblarr(datacount)
array_by = dblarr(datacount)
array_bz = dblarr(datacount)
fromSecond = (long(dd) - 1) * 86400 + $
             (long(start_hour)*3600) + $
             (long(start_min)*60) + $
             start_second
toSecond = long(fromSecond) + long(duration_sec)

if system_type eq 'UNH' then begin

   file_found = file_test(database_path)
   if file_found eq 0 then return
   
   system_gain = 4.43         ;system gain = 4.43 V/(nT*Hz)
   file_array = readmud(database_path, fromSecond, toSecond, sample_rate)
   ;check for no data in day
   if abs(min(file_array[1,*])) lt 1.d-6 and $
      abs(max(file_array[1,*])) lt 1.d-6 then return
   array_bx = file_array[0, *] / system_gain ;x signal in terms of nT
   array_by = file_array[1, *] / system_gain ;y signal in terms of nT
   ;->Change data in V (original data format) to in nT by dividing by the system gain ->applies only to UNH ULF system.

endif

;10 Hz MCM and SPA data from before November 25, 2014
new_harddrive = '/media/mirl-ulf/4TB_HDD1/'
if system_type eq 'AUG' then begin
   
   doy = ulf_year_month_day_to_doy(year_str, month, day_str)

   if station_in eq 'MCM' then $
      database_path = new_harddrive + 'mirl/ULF/incoming/database/' + $
                      station_in + '_converted/mc' + $
                      year_str_2 + $
                      strtrim(string(doy), 1) + $
                      '_01.dat'

   if station_in eq 'SPA' then $
      database_path = new_harddrive + '/mirl/ULF/incoming/database/' + $
                      station_in + '_converted/sp' + $
                      year_str_2 + $
                      strtrim(string(doy), 1) + $
                      '_01.dat'

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

   if station_in eq 'MCM' then $
      ulf_read_sp, database_path, data_count, times, file_array, 2
   if station_in eq 'SPA' then $
      ulf_read_sp, database_path, data_count, times, file_array, 3
   array_bx = file_array[0, *]
   array_by = file_array[1, *]
   if has_z eq 1 then array_bz = file_array[2, *]

endif

if system_type eq 'NJIT' then begin

   database_path = new_harddrive + '/mirl/ULF/incoming/database/' + $
                   station_in + '_converted/' + $
                   strtrim(string(yyyy), 1) + '_' + $
                   strtrim(string(mm), 1) + '_' + $
                   strtrim(string(dd), 1) + '_' + $
                   station_in + '_dbdt_v2.txt'

   if station_in eq 'HAL' then $
      database_path = new_harddrive + '/mirl/ULF/incoming/database/' + $
                      station_in + '_converted/' + $
                      station_in + '_' + $
                      strtrim(string(yyyy), 1) + '_' + $
                      strtrim(string(mm), 1) + '_' + $
                      strtrim(string(dd), 1) + $
                      '.txt'

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

   if station_in eq 'MCM' then $
      ulf_read_sp, database_path, data_count, times, file_array, 2
   if station_in eq 'JBS' then $
      ulf_read_sp, database_path, data_count, times, file_array, 2
   if station_in eq 'KSS' then $
      ulf_read_neu, database_path, data_count, times, file_array, 2
   if station_in eq 'VNA' then $
      ulf_read_neu, database_path, data_count, times, file_array, 2
   if station_in eq 'SPA' then $
      ulf_read_sp, database_path, data_count, times, file_array, 3
   if station_in eq 'HAL' then begin
      if year_int gt 2018 then begin
         ;data in 2019 sampled at 20 Hz?
         ulf_read_neu, database_path, data_count, times, file_array, 3
      endif else begin
         ulf_read_sp, database_path, data_count, times, file_array, 3
      endelse
   endif
   array_bx = file_array[0, *]
   array_by = file_array[1, *]
   if has_z eq 1 then array_bz = file_array[2, *]

endif

if system_type eq 'AAL-PIP' then begin

   pg_restore, date, station_in, array_bx, array_by

   ;check for no data
   if abs(min(array_by)) lt 1.d-6 and $
      abs(max(array_by)) lt 1.d-6 then return

endif

spec_array = dblarr(2, datacount)
;set freq in skt file
maxfreq_f = float(maxfreq)
;for i = 0, nrows-1 do print, '    [',strcompress(string(i+1), /remove_all),'] = ', maxfreq_f * i / nrows + 1. / nrows / 2.

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

skeleton_base = ' -delete ./'
cdf_file = new_harddrive + '/mirl/ULF/cdf/' + station_in + $
           '/' + year_str + '/' + month + '/mica_ulf_' + $
           strlowcase(station_in) + '_' + date + '_v00'
if station_in eq 'VNA' then $
cdf_file = new_harddrive + '/mirl/ULF/cdf/' + station_in + $
           '/' + year_str + '/' + month + '/mica_ulf_' + $
           strlowcase(station_in) + '_' + date + '_v00'

if station_in eq 'VNA' then begin
   skt_name = 'mica_ulf_20hz'
endif else begin
   if has_z eq 0 then skt_name = 'mica_ulf'
   if has_z eq 1 then skt_name = 'mica_ulf_3axis'
endelse
if station_in eq 'HAL' and $
   year_int gt 2018 then begin
;   skt_name = 'mica_ulf_20hz'
   skt_name = 'mica_ulf_3axis'
endif

if station_in eq 'KSS' then begin
   skt_name = 'mica_ulf_20hz'
endif

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
if station_in eq 'JBS' then begin

   cdf_attput, cdf_id, 'PI_affiliation', 0, 'Kyung Hee University'
   cdf_attput, cdf_id, 'PI_name', 0, 'Khan-Hyuk Kim'

endif

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
if has_z eq 1 then cdf_varput, cdf_id, 'dbdt_z', array_bz, /zvariable

spec_times = findgen(ncols) * (duration * spm * mph) / ncols
cdf_varput, cdf_id, 'Epoch_spectra', $
            start_time_long64 + long64(spec_times) * 1000000000LL
spec_array[0, *] = db_times
spec_array[1, *] = array_bx
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

if has_z eq 1 then begin

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

endif

;make 5 Hz spectra
maxfreq       = 5
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

if has_z eq 1 then begin

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

endif

cdf_close, cdf_id

end