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 'HAL' 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_syowa, filepath, FromSecond, ToSecond, sample_rate
FILE_NOT_FOUND = 1
NOT_ENOUGH_DATA = 2
num_chan = 2
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
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]) * 3600. + $
float(line_parts[1]) * 60. + $
float(line_parts[2]))
if time ge FromSecond and time lt ToSecond then begin
data_array[0, data_count] = float(line_parts[3])
data_array[1, data_count] = float(line_parts[4])
if num_chan eq 3 then begin
data_array[2, data_count] = float(line_parts[5])
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_maw, filepath, FromSecond, ToSecond, sample_rate
FILE_NOT_FOUND = 1
NOT_ENOUGH_DATA = 2
file_found = file_test(filepath)
if file_found eq 0 then return, FILE_NOT_FOUND
openr, dat_fh, filepath, /get_lun
; allocate data array
num_chan = 2
Data_Array = fltarr(num_chan, (long(ToSecond)-long(FromSecond)) * sample_rate)
text_line = ''
;read header
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]) * 3600.)
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])
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
month_array = ["Jan","Feb","Mar","Apr","May","Jun",$
"Jul","Aug","Sep","Oct","Nov","Dec"]
month = month_array[fix(mm)-1]
num_chan = 2
;find the start of the day in tt2000
day_start = cdf_parse_tt2000(dd + '-' + month + '-' + 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 daisy_plot_maw, $
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/'
HAL = 1
VNA = 1
SYO = 1
MAW = 1
f_low = 0 ;Hz
f_high = maxfreq ;Hz
;====== 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
sample_rate_20 = 20. ;Hz
period_20 = 1./sample_rate_20 ;sec
datacount_20 = long(duration * sph * sample_rate_20)
times_20 = findgen(datacount_20) * period/sph
sample_rate_2 = 2. ;Hz
period_2 = 1./sample_rate_2 ;sec
datacount_2 = long(duration * sph * sample_rate_2)
times_2 = findgen(datacount_2) * period/sph
data_array_VNA = dblarr(3, datacount_20)
data_array_SYO = dblarr(3, datacount_20)
data_array_MAW = dblarr(2, 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)
missing = double(9999999.0) ;missing data value
yminor_count = 5
max_B_field = max_B_field ;nT*Hz (nT/s)
maxfreq = maxfreq ;highest frequency plotted, in Hz
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
;====== Set up parameters for FFT process
nfft_20 = nfft ;number of points per fft: 1024 -> 0.47 mHz
step_20 = step / 2. ;step size between fft's (secs)
;Determine spectral array size and number of arrays for spectral images
skip_20 = round(step_20 / period_20)
ncols_20 = round((duration * spm * mph) / period_20)
if ((ncols_20 mod skip_20) ne 0) then print,'Error in step size!!!!!!!'
print, 'skip, period, ncols ', skip_20, period_20, ncols_20
ncols_20 = ncols_20 / skip_20
hzperbin_20 = 1/ (nfft_20 * period_20)
nrows_20 = round(maxfreq / hzperbin_20)
maxfreq_calc_20 = hzperbin_20 * nrows_20
spectra_20 = dblarr(ncols_20, nrows_20)
;====== Set up parameters for FFT process - 2 Hz
nfft_2 = nfft ;number of points per fft: 1024 -> 0.47 mHz
step_2 = step ;step size between fft's (secs)
;Determine spectral array size and number of arrays for spectral images
skip_2 = round(step_2 / period_2)
ncols_2 = round((duration * spm * mph) / period_2)
if ((ncols_2 mod skip_2) ne 0) then print,'Error in step size!!!!!!!'
print, 'skip, period, ncols ', skip_2, period_2, ncols_2
ncols_2 = ncols_2 / skip_2
hzperbin_2 = 1/ (nfft_2 * period_2)
nrows_2 = round(maxfreq / hzperbin_2)
maxfreq_calc_2 = hzperbin_2 * nrows_2
spectra_2 = dblarr(ncols_2, nrows_2)
;====== Read in data (stored in database) from the stations selected
;to be plotted
option_0 = 0
option_1 = 0
option_2 = 1
;stop
if option_0 eq 1 then begin
num_panels = 5
station_names = 'VNA_HAL_SYO_0'
endif
if option_1 eq 1 then begin
num_panels = 3
station_names = 'VNA_HAL_SYO_1'
endif
if option_2 eq 1 then begin
num_panels = 4
station_names = 'VNA_HAL_SYO_MAW_2'
endif
plot_array_hal = dblarr(datacount_20)
plot_array_neu = dblarr(datacount_20)
plot_array_syo = dblarr(datacount_20)
plot_array_maw = dblarr(datacount_2)
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)
times = findgen(datacount) / 36000. + $
((long(start_hour)*3600) + (long(start_min)*60) + start_second) / 3600.
times_20 = findgen(datacount*2) / 72000. + $
((long(start_hour)*3600) + (long(start_min)*60) + start_second) / $
3600.
times_2 = findgen(datacount_2) / 7200. + $
((long(start_hour)*3600) + (long(start_min)*60) + start_second) / $
3600.
system_gain = 4.43 ;system gain = 4.43 V/(nT*Hz) --> applies only to UNH ULF system
station_index = 0
;read HAL riometer
file_path = '~/DAISY/riometer_file_Z1R' + date + 'abs.txt'
ulf_read_riometer, file_path, rio_data_count, rio_times, rio_data_arr
if HAL eq 1 then begin
filepath = '/mirl/ULF/cdf/HAL/' + $
strtrim(string(yyyy), 1) + '/' + $
strtrim(string(mm), 1) + '/' + $
'mica_ulf_hal_' + $
strtrim(string(date), 1) + $
'_v00.cdf'
file_found = file_test(filepath)
if file_found eq 0 then stop
temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate_20, yyyy, mm, dd)
for i_fudge = 0, datacount_20-1 do $
plot_array_hal[i_fudge] = temp_array[0, i_fudge]
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 0 then stop
temp_array = read_cdf(filepath, fromSecond_sp, toSecond_sp, sample_rate_20, yyyy, mm, dd)
for i_fudge = 0, datacount_20-1 do $
plot_array_neu[i_fudge] = temp_array[0, i_fudge]
endif
if SYO eq 1 then begin
filepath = '/mirl/ULF/incoming/database/SYO/' + $
'SYO_ULF_' + $
strtrim(string(date), 1) + $
'.csv'
file_found = file_test(filepath)
if file_found eq 0 then stop
temp_array = read_syowa(filepath, fromSecond_sp, toSecond_sp, 20)
; for i_fudge = 0, n_elements(temp_array[0, *])-1 do $
; plot_array_syo[i_fudge] = $
; temp_array[0, i_fudge]
;DELTA? bugger
for i_fudge = 0, n_elements(temp_array[0, *])-2 do $
plot_array_syo[i_fudge] = $
temp_array[0, i_fudge+1] - temp_array[0, i_fudge]
endif
;read SYOWA riometer
file_path = '~/DAISY/SYOWA_r_' + date
ulf_read_riometer_syowa, file_path, rio_data_count_s, rio_times_s, rio_data_arr_s
if MAW eq 1 then begin
filepath = '/home/mirl/mawson/mawpul_27feb2019.dat'
file_found = file_test(filepath)
if file_found eq 0 then stop
temp_array = read_maw(filepath, fromSecond_sp, toSecond_sp, 2)
for i_fudge = 0, datacount_2-1 do $
plot_array_maw[i_fudge] = temp_array[0, i_fudge]
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(num_panels)
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(num_panels, 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)
spec_array_20 = dblarr(2, datacount*2)
spec_array_2 = dblarr(2, datacount_2)
if option_0 eq 1 then begin
title = 'Global ULF ' + date
xyouts, 0.45, 0.95, title, /normal, alignment=0.5, charsize=charsize2
fromHour = fromSecond_sp / 3600.
toHour = toSecond_sp / 3600.
plot_index = 0
plot, [fromhour, tohour], [-.25, 1.5], /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= ' HAL Absorption '
oplot, rio_times[fromSecond_sp:toSecond_sp-1], $
rio_data_arr[0,fromSecond_sp:toSecond_sp-1]
;====== Plot FFT spectrogram HAL BX
plot_index = 1
plot, [0, duration], [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= 'HAL Bx Freq (Hz)'
spec_array_20(0, *) = times_20 - fromHour
spec_array_20(1, *) = plot_array_hal
calcspec, spec_array_20, spectra_20, duration, period_20, sph, $
nfft_20, pstart, missing
greyplot_multi, spectra_20(*,*), positions[plot_index,*], winsize, $
doy, times_20, pstart, duration, missing, $
upper, lower, /colorbar
;====== Plot FFT spectrogram VNA BX
plot_index = 2
plot, [0, duration], [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= 'VNA Bx Freq (Hz)'
spec_array_20(0, *) = times_20 - fromHour
spec_array_20(1, *) = plot_array_neu
calcspec, spec_array_20, spectra_20, duration, period_20, sph, $
nfft_20, pstart, missing
greyplot_multi, spectra_20(*,*), positions[plot_index,*], winsize, $
doy, times_20, pstart, duration, missing, $
upper, lower
plot_index = 3
plot, [fromhour, tohour], [-.3,-.15], /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= ' SYO V '
oplot, rio_times_s[fromSecond_sp:toSecond_sp-1], $
rio_data_arr_s[fromSecond_sp:toSecond_sp-1]
;====== Plot FFT spectrogram SYO BX
plot_index = 4
plot, [0, duration], [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= 'SYO Bx Freq (Hz)'
spec_array_20(0, *) = times_20 - fromHour
spec_array_20(1, *) = plot_array_syo
calcspec, spec_array_20, spectra_20, duration, period_20, sph, $
nfft_20, pstart, missing
syo_freq = [0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05]
syo_gain = [1.030, 3.050, 4.987, 6.890, 8.683, 10.290, $
11.890, 13.453, 14.776, 15.791, 17.084]
spectra_freq = findgen(nrows) / nrows
syo_cal = spline(syo_freq, syo_gain, spectra_freq)
;Apply calibration
;for i_time = 0, n_elements(spectra_20[*, 0])-1 do $
; spectra_20[i_time, *] = spectra_20[i_time, *] / syo_cal
;stop
greyplot_multi, spectra_20(*,*), positions[plot_index,*], winsize, $
doy, times_20, pstart, duration, missing, $
upper, lower
endif
if option_1 eq 1 then begin
title = 'Global ULF ' + date
xyouts, 0.45, 0.95, title, /normal, alignment=0.5, charsize=charsize2
plot_index = 0
;====== Plot time-series (Bx)
plot, times_20, plot_array_hal, /noerase, $
nsum = nsum_value, pos=positions[plot_index,*], $
xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
ystyle=1, yticklen=-0.01, ycharsize=charsize1, $
yrange=[-max_B_field, max_B_field], $
; yrange=[-max_B_field, max_B_field], $
ytitle = 'HAL Bx dB/dt (nT/sec)'
plot_index = 1
;====== Plot time-series (Bx)
plot, times_20, plot_array_neu, /noerase, $
nsum = nsum_value, pos=positions[plot_index,*], $
xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
ystyle=1, yticklen=-0.01, ycharsize=charsize1, $
yrange=[-max_B_field, max_B_field], $
ytitle = 'VNA Bx dB/dt (nT/sec)'
plot_index = 2
;====== Plot time-series (Bx)
plot, times_20, plot_array_syo, /noerase, $
nsum = nsum_value, pos=positions[plot_index,*], $
xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
ystyle=1, yticklen=-0.01, ycharsize=charsize1, $
yrange=[-max_B_field, max_B_field], $
ytitle = 'SYO Bx dB/dt (nT/sec)'
endif
if option_2 eq 1 then begin
title = 'Global ULF ' + date
xyouts, 0.45, 0.95, title, /normal, alignment=0.5, charsize=charsize2
fromHour = fromSecond_sp / 3600.
toHour = toSecond_sp / 3600.
;====== Plot FFT spectrogram HAL BX
plot_index = 0
plot, [0, duration], [0,1], /noerase, /nodata, $
pos=positions[plot_index,*],$
xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
xticks=xtick_count, xminor=xminor_count, xtickname=xtick_labels, $
; xtickv = [-1,25], xtickname=[' ',' '], $ ;xtick_labels_blank, $
yticks=freq_ytickcount, yminor= yminor_count, $
ytickname=freq_yticklabel, $
ystyle=1, yticklen=-0.01, ycharsize=charsize1,$
ytitle= 'HAL Bx Freq (Hz)'
spec_array_20(0, *) = times_20 - fromHour
spec_array_20(1, *) = plot_array_hal
calcspec, spec_array_20, spectra_20, duration, period_20, sph, $
nfft_20, pstart, missing
greyplot_multi, spectra_20(*,*), positions[plot_index,*], winsize, $
doy, times_20, pstart, duration, missing, $
upper, lower;, /colorbar
;====== Plot FFT spectrogram VNA BX
plot_index = 1
plot, [0, duration], [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= 'VNA Bx Freq (Hz)'
spec_array_20(0, *) = times_20 - fromHour
spec_array_20(1, *) = plot_array_neu
calcspec, spec_array_20, spectra_20, duration, period_20, sph, $
nfft_20, pstart, missing
greyplot_multi, spectra_20(*,*), positions[plot_index,*], winsize, $
doy, times_20, pstart, duration, missing, $
upper, lower
;====== Plot FFT spectrogram SYO BX
plot_index = 2
plot, [0, duration], [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= 'SYO Bx Freq (Hz)'
spec_array_20(0, *) = times_20 - fromHour
spec_array_20(1, *) = plot_array_syo
calcspec, spec_array_20, spectra_20, duration, period_20, sph, $
nfft_20, pstart, missing
syo_freq = [0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05]
syo_gain = [1.030, 3.050, 4.987, 6.890, 8.683, 10.290, $
11.890, 13.453, 14.776, 15.791, 17.084]
spectra_freq = findgen(nrows) / nrows
syo_cal = spline(syo_freq, syo_gain, spectra_freq)
;Apply calibration
;for i_time = 0, n_elements(spectra_20[*, 0])-1 do $
; spectra_20[i_time, *] = spectra_20[i_time, *] / syo_cal
;stop
greyplot_multi, spectra_20(*,*), positions[plot_index,*], winsize, $
doy, times_20, pstart, duration, missing, $
upper, lower
;====== Plot FFT spectrogram MAW BX
plot_index = 3
plot, [0, duration], [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= 'MAW Bx Freq (Hz)'
stop
spec_array_2(0, *) = times_2 - fromHour
spec_array_2(1, *) = plot_array_maw
calcspec, spec_array_2, spectra_2, duration, period_2, sph, $
nfft_2, pstart, missing
greyplot_multi, spectra_2(*,*), positions[plot_index,*], winsize, $
doy, times_2, pstart, duration, missing, $
upper, lower
endif
;==========================================================================
; Reset graphic device
print, file_path_base
device, /close
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