pro ulf_make_mud, station_id, start_date, file_list

;Created:
;20161013 Mark Chutter
;
;Args:
;station_id 3 character station ID
;start_date first day of month in the form yymmdd
;file_list  list of .dat from which data will be extracted
;
;Description:
;I think MUD stands for Monthly ULF Database
;ulf_make_mud takes a list of input .dat files and makes a monthly MUD
;file whose first data point is the first point after midnight on the
;first day of the month - no data is padded with 0

;it appears the MUD format includes a 32 byte header header:
;file_flag = bytarr(4) 3 letter station ID and NULL
;num_channels = 0
;sample_rate = 0
;sample_size = 0
;start_year = 0
;start_month = 0
;start_day = 0
;station_name = bytarr(4) 3 letter station ID and NULL
;num_seconds = long(0)
;reserved1 = long(0)
;reserved2 = long(0)

;example:
;f_0 = findfile('/mirl/ULF/incoming/database/LYR/053122*.dat', count=count)
;f_1 = findfile('/mirl/ULF/incoming/database/LYR/06??22*.dat', count=count)
;ulf_make_mud, 'LYR', '20220601', [f_0, f_1]

mud_dir = '/mirl/ULF/Database/'

;60 days may seem like overkill, but when daily .dat files contain
;many days, trouble can arise
NUM_TIMES = 5184000L ;86400 seconds * 60 days
NUM_POINTS = 51840000L ;86400 seconds * 10 per second * 60 days
mud_time = lonarr(NUM_TIMES)
mud_data = intarr(2, NUM_POINTS)

file_flag_bytes = bytarr(4) 
file_flag_bytes[0:2] = byte('MUD') 
file_flag_bytes[3] = 0

num_channels = 2
sample_rate = 10 ;in Hz
sample_size = 16 ;in bits

year = strmid(start_date, 0, 4)
year_short = strmid(start_date, 2, 2)
month = strmid(start_date, 4, 2)
day = strmid(start_date, 6, 2)
year_short_long = long(year_short);luv the name
month_long = long(month)
day_long = long(day)

if month eq '01' or $
   month eq '03' or $
   month eq '05' or $
   month eq '07' or $
   month eq '08' or $
   month eq '10' or $
   month eq '12' then num_days = 31

if month eq '04' or $
   month eq '06' or $
   month eq '09' or $
   month eq '11' then num_days = 30

if month eq '02' then begin

   if year eq '2000' or $
      year eq '2004' or $
      year eq '2008' or $
      year eq '2012' or $
      year eq '2016' or $
      year eq '2020' or $
      year eq '2024' then num_days = 29 $
   else num_days = 28

endif

station_name_bytes = bytarr(4) 
station_name_bytes[0:2] = byte(station_id) 
station_name_bytes[3] = 0

year_fix = fix(year)
month_fix = fix(month)
day_fix = fix(day)
num_seconds = num_days * 86400L
reserved1 = 0L
reserved2 = 0L

mud_name = mud_dir + station_id + '_' + year + '_' + month + '.MUD' 
print, mud_name
openw, mud_lun, mud_name, /get_lun
writeu, mud_lun, file_flag_bytes
writeu, mud_lun, num_channels
writeu, mud_lun, sample_rate
writeu, mud_lun, sample_size
writeu, mud_lun, year_fix
writeu, mud_lun, month_fix
writeu, mud_lun, day_fix
writeu, mud_lun, station_name_bytes
writeu, mud_lun, num_seconds
writeu, mud_lun, reserved1
writeu, mud_lun, reserved2

for i=0, n_elements(file_list)-1 do begin
print, file_list[i]
filepath = file_list[i]
  ulf_read_dat, $
     filepath, $
     date_str, $
     time_str, $
     lat_str, $
     latd_str, $
     lon_str, $
     lond_str, $
     data_count, $
     data_arr

  ;get GPS day, month, year
  gps_day = long(strmid(date_str, 0, 2))
  gps_month = long(strmid(date_str, 2, 2))
  gps_year = long(strmid(date_str, 4, 2))
  ;get GPS hours, minutes, seconds
  gps_hour = long(strmid(time_str, 0, 2))
  gps_min = long(strmid(time_str, 2, 2))
  gps_sec = long(strmid(time_str, 4, 2))

  if year_short_long gt 01 then begin

     if (month_long gt gps_month) or $
        (month_long eq 1L and gps_month eq 12L) then begin

        time0 = 3600L * gps_hour + $
                60L * gps_min + $
                gps_sec

        early_secs = 86400L - time0
        mud_secs = (data_count / 10L) - early_secs
        mud_samps = data_count - (early_secs * 10L)
        data_offset = early_secs * 10L

        gps_sync = ishft((data_arr[0, data_offset] AND '8000'x), -15)
        mud_time[0:mud_secs-1L] = $
           lindgen(mud_secs)

        mud_data[0, 0:mud_samps-1L] = $
           fix(data_arr[0, data_offset:data_count-1L] AND '0fff'x) - 2048
        mud_data[1, 0:mud_samps-1L] = $
           fix(data_arr[1, data_offset:data_count-1L] AND '0fff'x) - 2048

     endif else if (month_long eq gps_month) then begin

        time_index0 = 86400L * (gps_day-1L) + $
                      3600L * (gps_hour) + $
                      60L * gps_min + $
                      gps_sec

        data_index0 = time_index0 * 10L
        mud_time[time_index0:time_index0+(data_count / 10L)-1L] = $
           time_index0 + lindgen(data_count / 10L)
        mud_data[0, data_index0:data_index0+data_count-1L] = $
           fix(data_arr[0, *] AND '0fff'x) - 2048
        mud_data[1, data_index0:data_index0+data_count-1L] = $
           fix(data_arr[1, *] AND '0fff'x) - 2048

     endif

  endif

endfor

for i = 0, num_seconds-1 do begin
   writeu, mud_lun, mud_time[i]
   out_data0 = mud_data[0, i * 10L: i * 10L + 9]
   writeu, mud_lun, out_data0
   out_data1 = mud_data[1, i * 10L: i * 10L + 9]
   writeu, mud_lun, out_data1
endfor

close, mud_lun
free_lun, mud_lun

end