function readmud, FilePath, FromSecond, ToSecond, sample_rate


;FilePath='c:\LYR_2007_01.MUD'
;FromSecond = long(long(3600) * 24)
;ToSecond = long(long(3600) * 24 * 2)

FILE_NOT_FOUND = 1
NOT_ENOUGH_DATA = 2

; open the MUD file:
file_found = file_test(filepath)
if file_found eq 0 then return, FILE_NOT_FOUND

openr, dat_fh, filepath, /get_lun

; read header
FileFlag = bytarr(4)
NumChannels = 0
SampleRate = 0
SampleSize = 0
StartYear = 0
StartMonth = 0
StartDay = 0
StationName = bytarr(4)
NumSeconds = long(0)
Reserved1 = long(0)
Reserved2 = long(0)

readu, dat_fh, FileFlag
readu, dat_fh, NumChannels
readu, dat_fh, SampleRate
readu, dat_fh, SampleSize
readu, dat_fh, StartYear
readu, dat_fh, StartMonth
readu, dat_fh, StartDay
readu, dat_fh, StationName
readu, dat_fh, NumSeconds
readu, dat_fh, Reserved1
readu, dat_fh, Reserved2

; set sample buffer size
SampleBufferSize = SampleRate * SampleSize / 8

; allocate data array
Data_Array  = fltarr(2, (long(ToSecond)-long(FromSecond)) * sample_rate)

; allocate temp buffer
Buffer_Array = intarr(10, 2)

; move to the 1st data position
Offset = long(32)
Offset = Offset + FromSecond * (4 + SampleBufferSize * NumChannels)
POINT_LUN, dat_fh, Offset

; read samples
Second = long(0)
FOR n = long(0), long(ToSecond)-long(FromSecond)-1 DO BEGIN
    readu, dat_fh, Second

    ; read 10-sample chunk of 2 channels
    readu, dat_fh, Buffer_Array

    ; assign the samples to the data array
    FOR m = 0, 9 DO BEGIN
        Data_Array(0, n * 10 + m) = (Buffer_Array(m, 0) + 2048.0) / 4096.0 * 20.0 - 10.0;
        Data_Array(1, n * 10 + m) = (Buffer_Array(m, 1) + 2048.0) / 4096.0 * 20.0 - 10.0;
    ENDFOR
ENDFOR

; close the file
close, dat_fh
free_lun, dat_fh

return, Data_Array
end



function read_sp, sp_station, filepath, FromSecond, ToSecond, sample_rate

FILE_NOT_FOUND = 1
NOT_ENOUGH_DATA = 2

if sp_station eq 'SPA' then begin

   num_chan = 3

endif

if sp_station eq 'MCM' then begin

   num_chan = 2

endif

file_found = file_test(filepath)
if file_found eq 0 then return, FILE_NOT_FOUND

openr, dat_fh, filepath, /get_lun

; allocate data array
Data_Array  = fltarr(num_chan, (long(ToSecond)-long(FromSecond)) * sample_rate)
text_line = ''
;read header
readf, dat_fh, text_line
readf, dat_fh, text_line
readf, dat_fh, text_line

data_count = 0
done = -1
while not eof(dat_fh) and done eq -1 do begin

   readf, dat_fh, text_line
   line_parts = strsplit(strcompress(text_line), ' ', /extract)
   time = long(float(line_parts[0]))

   if time ge FromSecond and time lt ToSecond then begin

      data_array[0, data_count] = float(line_parts[1])
      data_array[1, data_count] = float(line_parts[2])
      if num_chan eq 3 then begin
         data_array[2, data_count] = float(line_parts[3])
      endif

      data_count = data_count + 1L

   endif

   if time ge ToSecond then done = 1

endwhile

; close the file
close, dat_fh
free_lun, dat_fh

return, Data_Array
end


function read_cdf, filepath, FromSecond, ToSecond, sample_rate, yyyy, mm, dd

file_found = file_test(filepath)
if file_found eq 0 then return, FILE_NOT_FOUND

num_chan = 2
;find the start of the day in tt2000
day_start = cdf_parse_tt2000(dd + '-' + 'Jan' + '-' + yyyy + ' 00:00:00.000')

cdf_id = cdf_open(filepath)

cdf_control, cdf_id, var = 'dbdt_x', get_var_info = rsurvey
nrec_in = rsurvey.maxrec+1
datacount = nrec_in
cdf_varget, cdf_id, 'Epoch', epoch, $
  rec_start = 0L, rec_count = nrec_in
epoch = reform(epoch)
cdf_varget, cdf_id, 'dbdt_x', dbdt_x, $
  rec_start = 0L, rec_count = nrec_in
dbdt_x = reform(dbdt_x)
cdf_varget, cdf_id, 'dbdt_y', dbdt_y, $
  rec_start = 0L, rec_count = nrec_in
dbdt_y = reform(dbdt_y)

w_range = where(epoch ge day_start + FromSecond * 1000000000LL and $
                epoch lt day_start + ToSecond * 1000000000LL, count)

Data_Array  = fltarr(num_chan, count)
Data_Array[0, *] = dbdt_x(w_range)
Data_Array[1, *] = dbdt_y(w_range)
cdf_close, cdf_id

return, data_array
end



;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: web_plot, '20161031', '000000', '240000',
;plot_dir, 'ps', '100000000', 1, 2048, 200, 0.2, 0, 1
pro web_plot_x, $
   date, $
   start_time, $
   end_time, $
   plot_dir, $
   file_format, $
   stations, $
   maxfreq, $
   nfft, $
   step, $
   max_B_field, $
   f_low, $
   f_high

FILE_NOT_FOUND = 1
NOT_ENOUGH_DATA = 2

sp_dir = '/mirl/ULF/incoming/database/'

database_dir = '/mirl/ULF/Database/'

;determine which stations were requested from stations string
NAL = 0
LYR = 0
HOR = 0
SDY = 0
SNK = 0
IQA = 0
ISR = 0
SPA = 0
MCM = 0
PG2 = 0
PG3 = 0
PG4 = 0
PG5 = 0
JBS = 0
VNA = 0

if strmid(stations, 0, 1) eq '1' then NAL = 1
if strmid(stations, 1, 1) eq '1' then LYR = 1
if strmid(stations, 2, 1) eq '1' then HOR = 1
if strmid(stations, 3, 1) eq '1' then SDY = 1
if strmid(stations, 4, 1) eq '1' then SNK = 1
if strmid(stations, 5, 1) eq '1' then IQA = 1
if strmid(stations, 6, 1) eq '1' then ISR = 1
if strmid(stations, 7, 1) eq '1' then SPA = 1
if strmid(stations, 8, 1) eq '1' then MCM = 1
if strmid(stations, 9, 1) eq '1' then PG2 = 1
if strmid(stations, 10, 1) eq '1' then PG3 = 1
if strmid(stations, 11, 1) eq '1' then PG4 = 1
if strmid(stations, 12, 1) eq '1' then PG5 = 1
if strmid(stations, 13, 1) eq '1' then JBS = 1
if strmid(stations, 14, 1) eq '1' then VNA = 1

f_low         = 0             ;Hz
f_high        = maxfreq   ;Hz


;====== Create station names for file name
station_names = strarr(1)
if NAL eq 1 then begin
	station_names = 'N'
endif
if LYR eq 1 then begin
	station_names = station_names + 'L'
endif
if HOR eq 1 then begin
	station_names = station_names + 'H'
endif
if SDY eq 1 then begin
	station_names = station_names + 'S'
endif
if SNK eq 1 then begin
	station_names = station_names + 'P'
endif
if IQA eq 1 then begin
	station_names = station_names + 'I'
endif
if ISR eq 1 then begin
	station_names = station_names + 'R'
endif
if SPA eq 1 then begin
	station_names = station_names + 'A'
endif
if MCM eq 1 then begin
	station_names = station_names + 'M'
endif
if PG2 eq 1 then begin
	station_names = station_names + '2'
endif
if PG3 eq 1 then begin
	station_names = station_names + '3'
endif
if PG4 eq 1 then begin
	station_names = station_names + '4'
endif
if PG5 eq 1 then begin
	station_names = station_names + '5'
endif
if JBS eq 1 then begin
	station_names = station_names + 'J'
endif
if VNA eq 1 then begin
	station_names = station_names + 'N'
endif
print, station_names


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

ulf_tick_labels, double(start_time_sec), double(end_time_sec), xtick_count, tick_values, xtick_labels

;====== Specify MUD sources
database_path_NAL = database_dir + '/NAL_' + strtrim(string(yyyy), 1) + '_' + strtrim(string(mm), 1) + '.MUD'
database_path_LYR = database_dir + '/LYR_' + strtrim(string(yyyy), 1) + '_' + strtrim(string(mm), 1) + '.MUD'
database_path_HOR = database_dir + '/HOR_' + strtrim(string(yyyy), 1) + '_' + strtrim(string(mm), 1) + '.MUD'
database_path_SDY = database_dir + '/SDY_' + strtrim(string(yyyy), 1) + '_' + strtrim(string(mm), 1) + '.MUD'
database_path_SNK = database_dir + '/SNK_' + strtrim(string(yyyy), 1) + '_' + strtrim(string(mm), 1) + '.MUD'
database_path_IQA = database_dir + '/IQA_' + strtrim(string(yyyy), 1) + '_' + strtrim(string(mm), 1) + '.MUD'
database_path_ISR = database_dir + '/ISR_' + 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)
times           = findgen(datacount) * period/sph
data_array_NAL  = dblarr(3, datacount)
data_array_LYR  = dblarr(3, datacount)
data_array_HOR  = dblarr(3, datacount)
data_array_SDY  = dblarr(3, datacount)
data_array_SNK  = dblarr(3, datacount)
data_array_IQA  = dblarr(3, datacount)
data_array_ISR  = dblarr(3, datacount)
data_array_JBS  = dblarr(3, datacount)
data_array_VNA  = dblarr(3, datacount*2)

;====== 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
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!!!!!!!'
print, 'skip, period, ncols ', skip, period, ncols
ncols         = ncols/skip
hzperbin      = 1/ (nfft * period)
nrows         = round(maxfreq / hzperbin)
maxfreq_calc  = hzperbin * nrows
spectra       = dblarr(ncols,nrows)
help, spectra

;====== Read in data (stored in database) from the stations selected to be plotted
no_stations   = NAL + LYR + HOR + SDY + SNK + IQA + ISR + SPA + MCM + PG2 + PG3 + PG4 + PG5 + JBS + VNA
plot_array_Bx = dblarr(no_stations*2, datacount)
plot_array_By = dblarr(no_stations*2, datacount)
status_array = intarr(no_stations*2)

station_code  = strarr(no_stations)
fromSecond    = (long(dd) - 1) * 86400 + (long(start_hour)*3600) + (long(start_min)*60) + start_second
toSecond     = long(fromSecond) + long(duration_sec)

fromSecond_sp = (long(start_hour)*3600) + (long(start_min)*60) + start_second
toSecond_sp   = long(fromSecond_sp) + long(duration_sec)

system_gain   = 4.43  ;system gain = 4.43 V/(nT*Hz) --> applies only to UNH ULF system
station_index = 0

if NAL eq 1 then begin

  temp_array = readmud(database_path_NAL, fromSecond, toSecond, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     reformed_array = reform(temp_array[0, *])
     temp_array[0, *] = convol(reformed_array, coeff) ;digital filtering
     reformed_array = reform(temp_array[1, *])
     temp_array[1, *] = convol(reformed_array, coeff) ;digital filtering
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] / system_gain ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_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.

  endelse

  station_code(station_index/2)  = 'NAL'
  station_index = station_index + 2

endif

if LYR eq 1 then begin

  temp_array = readmud(database_path_LYR, fromSecond, toSecond, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     reformed_array = reform(temp_array[0, *])
     temp_array[0, *] = convol(reformed_array, coeff) ;digital filtering
     reformed_array = reform(temp_array[1, *])
     temp_array[1, *] = convol(reformed_array, coeff) ;digital filtering
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] / system_gain ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_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.

  endelse

  station_code(station_index/2)  = 'LYR'
  station_index = station_index + 2

endif

if HOR eq 1 then begin

  temp_array = readmud(database_path_HOR, fromSecond, toSecond, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     reformed_array = reform(temp_array[0, *])
     temp_array[0, *] = convol(reformed_array, coeff) ;digital filtering
     reformed_array = reform(temp_array[1, *])
     temp_array[1, *] = convol(reformed_array, coeff) ;digital filtering
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] / system_gain ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_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.

  endelse

  station_code(station_index/2)  = 'HOR'
  station_index = station_index + 2
   
endif

if SDY eq 1 then begin

  temp_array = readmud(database_path_SDY, fromSecond, toSecond, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     reformed_array = reform(temp_array[0, *])
     temp_array[0, *] = convol(reformed_array, coeff) ;digital filtering
     reformed_array = reform(temp_array[1, *])
     temp_array[1, *] = convol(reformed_array, coeff) ;digital filtering
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] / system_gain ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_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.

  endelse

  station_code(station_index/2)  = 'SDY'
  station_index = station_index + 2

endif

if SNK eq 1 then begin

  temp_array = readmud(database_path_SNK, fromSecond, toSecond, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     reformed_array = reform(temp_array[0, *])
     temp_array[0, *] = convol(reformed_array, coeff) ;digital filtering
     reformed_array = reform(temp_array[1, *])
     temp_array[1, *] = convol(reformed_array, coeff) ;digital filtering
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] / system_gain ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_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.

  endelse

  station_code(station_index/2)  = 'SNK'
  station_index = station_index + 2

endif

if IQA eq 1 then begin

  temp_array = readmud(database_path_IQA, fromSecond, toSecond, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     reformed_array = reform(temp_array[0, *])
     temp_array[0, *] = convol(reformed_array, coeff) ;digital filtering
     reformed_array = reform(temp_array[1, *])
     temp_array[1, *] = convol(reformed_array, coeff) ;digital filtering
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] / system_gain ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_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.

  endelse

  station_code(station_index/2)  = 'IQA'
  station_index = station_index + 2

endif

if ISR eq 1 then begin

  temp_array = readmud(database_path_ISR, fromSecond, toSecond, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     reformed_array = reform(temp_array[0, *])
     temp_array[0, *] = convol(reformed_array, coeff) ;digital filtering
     reformed_array = reform(temp_array[1, *])
     temp_array[1, *] = convol(reformed_array, coeff) ;digital filtering
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] / system_gain ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_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.

  endelse

  station_code(station_index/2)  = 'ISR'
  station_index = station_index + 2

endif

if SPA eq 1 then begin

   filepath = sp_dir + 'SPA_converted/' + $
              strtrim(string(yyyy), 1) + '_' + $
              strtrim(string(mm), 1) + '_' + $
              strtrim(string(dd), 1) + '_SPA_dbdt_v2.txt'

  temp_array = read_sp('SPA', filepath, fromSecond_sp, toSecond_sp, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endelse

  station_code(station_index/2)  = 'SPA'
  station_index = station_index + 2

endif


if MCM eq 1 then begin

   filepath = sp_dir + 'MCM_converted/' + $
              strtrim(string(yyyy), 1) + '_' + $
              strtrim(string(mm), 1) + '_' + $
              strtrim(string(dd), 1) + '_MCM_dbdt_v2.txt'

  temp_array = read_sp('MCM', filepath, fromSecond_sp, toSecond_sp, sample_rate)

  if n_elements(temp_array) eq 1L then begin

     if temp_array eq FILE_NOT_FOUND then begin

        status_array[station_index] = 0
        status_array[station_index+1] = 0

     endif

  endif else begin

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endelse

  station_code(station_index/2)  = 'MCM'
  station_index = station_index + 2

endif

if PG2 eq 1 then begin

   filepath = '/mirl/ULF/cdf/PG2/' + $
              strtrim(string(yyyy), 1) + '/' + $
              strtrim(string(mm), 1) + '/' + $
              'mica_ulf_pg2_' + $
              strtrim(string(date), 1) + $
              '_v00.cdf'

   file_found = file_test(filepath)
   if file_found eq 1 then begin

      temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate, yyyy, mm, dd)

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endif else begin

     status_array[station_index] = 0
     status_array[station_index+1] = 0

  endelse

  station_code(station_index/2)  = 'PG2'
  station_index = station_index + 2

endif


if PG3 eq 1 then begin

   filepath = '/mirl/ULF/cdf/PG3/' + $
              strtrim(string(yyyy), 1) + '/' + $
              strtrim(string(mm), 1) + '/' + $
              'mica_ulf_pg3_' + $
              strtrim(string(date), 1) + $
              '_v00.cdf'

   file_found = file_test(filepath)
   if file_found eq 1 then begin

      temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate, yyyy, mm, dd)

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endif else begin

     status_array[station_index] = 0
     status_array[station_index+1] = 0

  endelse

  station_code(station_index/2)  = 'PG3'
  station_index = station_index + 2

endif


if PG4 eq 1 then begin

   filepath = '/mirl/ULF/cdf/PG4/' + $
              strtrim(string(yyyy), 1) + '/' + $
              strtrim(string(mm), 1) + '/' + $
              'mica_ulf_pg4_' + $
              strtrim(string(date), 1) + $
              '_v00.cdf'

   file_found = file_test(filepath)
   if file_found eq 1 then begin

      temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate, yyyy, mm, dd)

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endif else begin

     status_array[station_index] = 0
     status_array[station_index+1] = 0

  endelse

  station_code(station_index/2)  = 'PG4'
  station_index = station_index + 2

endif


if PG5 eq 1 then begin

   filepath = '/mirl/ULF/cdf/PG5/' + $
              strtrim(string(yyyy), 1) + '/' + $
              strtrim(string(mm), 1) + '/' + $
              'mica_ulf_pg5_' + $
              strtrim(string(date), 1) + $
              '_v00.cdf'

   file_found = file_test(filepath)
   if file_found eq 1 then begin

      temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate, yyyy, mm, dd)

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endif else begin

     status_array[station_index] = 0
     status_array[station_index+1] = 0

  endelse

  station_code(station_index/2)  = 'PG5'
  station_index = station_index + 2

endif


if JBS eq 1 then begin

   filepath = '/mirl/ULF/cdf/JBS/' + $
              strtrim(string(yyyy), 1) + '/' + $
              strtrim(string(mm), 1) + '/' + $
              'mica_ulf_jbs_' + $
              strtrim(string(date), 1) + $
              '_v00.cdf'

   file_found = file_test(filepath)
   if file_found eq 1 then begin

      temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate, yyyy, mm, dd)

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     plot_array_Bx[station_index+1, *] = temp_array[0, *] ;x signal in terms of nT
     plot_array_By[station_index, *] = times
     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endif else begin

     status_array[station_index] = 0
     status_array[station_index+1] = 0

  endelse

  station_code(station_index/2)  = 'JBS'
  station_index = station_index + 2

endif


if VNA eq 1 then begin

   filepath = '/mirl/ULF/cdf/VNA/' + $
              strtrim(string(yyyy), 1) + '/' + $
              strtrim(string(mm), 1) + '/' + $
              'mica_ulf_neu_' + $
              strtrim(string(date), 1) + $
              '_v00.cdf'

   file_found = file_test(filepath)
   if file_found eq 1 then begin

      temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate, yyyy, mm, dd)

     status_array[station_index] = 1
     status_array[station_index+1] = 1
     plot_array_Bx[station_index, *] = times
     for i_fudge = 0, 863999 do $
        plot_array_Bx[station_index+1, i_fudge] = $
        (temp_array[0, i_fudge * 2] + temp_array[0, i_fudge * 2 + 1]) /2.
;     plot_array_By[station_index, *] = times
;     plot_array_By[station_index+1, *] = temp_array[1, *] ;y signal in terms of nT

  endif else begin

     status_array[station_index] = 0
     status_array[station_index+1] = 0

  endelse

  station_code(station_index/2)  = 'VNA'
  station_index = station_index + 2

endif



;====== Specify graphic format
;file_format = 'x'
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_names + graphic_type
file_path_plot   = plot_dir + '/' + file_name_plot
file_path_base   = plot_dir + '/' + $
                   date + '_' + start_time + '_' + end_time +'_' + $
                   station_names


;====== 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.75, yoffset=0.5, xsize=7.0, ysize=9.5
;   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
;if !d.name eq 'WIN' then begin
;    window, 0, retain=2, xsize=1024 ,ysize=768, title=title
;    print, 'Device:', !d.name
;endif



;====== Set up the plotting window dimensions
xsize = 650.
ysize = 820.
start_x = 60.
stop_x = 580.
x_len = stop_x - start_x
n_plots = 16
y_gap = 0.2
top_space = 50
bottom_space = 40
rel_size = intarr(no_stations)
rel_size[*] = 1
layout_plots, rel_size, y_gap, ysize, top_space, bottom_space, plot_loc
y_len = plot_loc[1,0] - plot_loc[0,0]

;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

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

;position1 = [x0,y0,x1,y1]
positions = fltarr(no_stations, 4)
positions[*,0] = start_x / xsize
positions[*,1] = plot_loc[0,*] / ysize
positions[*,2] = stop_x / xsize
positions[*,3] = plot_loc[1,*] / ysize

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


spec_array = dblarr(2, datacount)

for plot_index = 0, station_index/2-1 do begin

   ;====== Plot FFT spectrogram (Bx)
   if plot_index eq 0 then begin
      plot, [0,max(times)], [0,1], /noerase, /nodata, $
            pos=positions[plot_index,*],$
            xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
            xticks=xtick_count, xminor=xminor_count, $
            xtickv = tick_values, xtickname=xtick_labels, $
            yticks=freq_ytickcount, yminor= yminor_count,  $
            ytickname=freq_yticklabel, $
            ystyle=1, yticklen=-0.01, ycharsize=charsize1,$
            ytitle= strtrim(station_code(plot_index),2) + ' Bx Freq (Hz)'
   endif else begin
      plot, [0,max(times)], [0,1], /noerase, /nodata, $
            pos=positions[plot_index,*],$
            xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
            xticks=1, xminor=xminor, $
            xtickv = [-1,25], xtickname=[' ',' '], $;xtick_labels_blank, $
            yticks=freq_ytickcount, yminor= yminor_count,  $
            ytickname=freq_yticklabel, $
            ystyle=1, yticklen=-0.01, ycharsize=charsize1,$
            ytitle= strtrim(station_code(plot_index),2) + ' Bx Freq (Hz)'
   endelse
   spec_array(0, *) = plot_array_Bx(plot_index*2, *)
   spec_array(1, *) = plot_array_Bx(plot_index*2+1, *)
   calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing

   if plot_index eq 0 then begin
      greyplot_multi, spectra(*,*), positions[plot_index,*], winsize, $
                      doy, times, pstart, duration, missing, $
                      upper, lower, /colorbar
   endif else begin
      greyplot_multi, spectra(*,*), positions[plot_index,*], winsize, $
                      doy, times, pstart, duration, missing, $
                      upper, lower
   endelse

   title = 'Global ULF ' + date
   xyouts, 0.45, 0.95, title, /normal, alignment=0.5, charsize=charsize2

endfor


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