;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: main_process, database_dir, '20161031', '000000', '240000', plot_dir, 'ps', NAL, LYR, HOR, SDY, POK, IQA, ISR, 1, 2048, 200, 0.2, 0, 1
pro main_process, $
database_dir, $
date, $
start_time, $
end_time, $
plot_dir, $
file_format, $
NAL, $
LYR, $
HOR, $
SDY, $
POK, $
IQA, $
ISR, $
maxfreq, $
nfft, $
step, $
max_B_field, $
f_low, $
f_high, $
upper, $
lower
;;====== The followings will be located in batchprocess.pro
database_dir = '/mirl/ULF/Database'
plot_dir = '/mirl/ULF/data_plot/test'
file_format = 'ps'
NAL = 1
LYR = 1
HOR = 0
SDY = 0
POK = 0
IQA = 0
ISR = 0
maxfreq = 1
nfft = 2^11
step = 300 ;sec
max_B_field = 0.2 ;nT/sec
f_low = 0 ;Hz
f_high = maxfreq ;Hz
;====== 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
;====== 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 POK 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
;====== 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.
;====== Set up X axis (time scale)
if duration_sec gt 3600. then begin
xtick_count = duration
if xtick_count ne fix(duration) then begin
print, 'Time interval error! Set up the time in hourly step.'
stop
endif
xminor_count = 12
xtick_labels = strarr(xtick_count + 1)
xtick_labels_temp = timegen(xtick_count + 1, units='hours', $
start = julday(mm, dd, yyyy, start_hour, start_min, start_second),$
final = julday(mm, dd, yyyy, end_hour, end_min, end_second))
caldat, xtick_labels_temp, month, day, year, hours, minutes, seconds
counts = size(xtick_labels_temp, /n_elements)
jump = fix(duration/4)
;The following lines are devided into (1) and (2) to debug the IDL's function, 'timegen'.
;This function generates only two vectors when the start time and end
;time are set at 13:00 and 15:00, respectively.
;On the other hand, if they are set at 14:00 and 16:00, the function generates three vectors.
;There are many other combinations that work or don't work. Hyomin 5/1/2008
if counts eq xtick_count+1 then begin ;----------------(1)
month = strtrim(string(month), 2)
day = strtrim(string(day),2)
year = strtrim(string(year),2)
hours = strtrim(string(hours),2)
minutes = strtrim(string(minutes),2)
seconds = strtrim(string(seconds),2)
counts = size(xtick_labels_temp, /n_elements)
if jump eq 0 then begin
jump = 1
endif
for x_index = 0, xtick_count do begin
if month(x_index) lt 10 then begin
month(x_index) = '0' + month(x_index) ;to make the label "09:00" rather than "9:00"
endif
if day(x_index) lt 10 then begin
day(x_index) = '0' + day(x_index)
endif
if hours(x_index) lt 10 then begin
hours(x_index) = '0' + hours(x_index)
endif
if minutes(x_index) lt 10 then begin
minutes(x_index) = '0' + minutes(x_index)
endif
if seconds(x_index) lt 10 then begin
seconds(x_index) = '0' + seconds(x_index)
endif
xtick_labels(x_index) = ' '
endfor
for x_index = 0, xtick_count, jump do begin
xtick_labels(x_index) = hours(x_index) + ':' + minutes(x_index)
end
endif else begin ;---------------- (2)
hours_fixed = intarr(counts+1)
minutes_fixed = intarr(counts+1)
for x_index = 0, counts-1 do begin
hours_fixed(x_index) = hours(x_index)
minutes_fixed(x_index) = minutes(x_index)
endfor
hours_fixed(counts) = hours(counts-1) + 1
minutes_fixed(counts) = minutes(counts-1)
hours_fixed = strtrim(string(hours_fixed),2)
minutes_fixed = strtrim(string(minutes_fixed),2)
for x_index = 0, counts do begin
if hours_fixed(x_index) lt 10 then begin
hours_fixed(x_index) = '0' + hours_fixed(x_index)
endif
if minutes_fixed(x_index) lt 10 then begin
minutes_fixed(x_index) = '0' + minutes_fixed(x_index)
endif
endfor
if jump eq 0 then begin
jump = 1
endif
for x_index = 0, xtick_count, jump do begin
xtick_labels(x_index) = hours_fixed(x_index) + ':' + minutes_fixed(x_index)
end
endelse
endif
if duration_sec eq 3600 then begin
xtick_labels_temp = timegen(units='minutes', $
start = julday(mm, dd, yyyy, start_hour, start_min, start_second),$
final = julday(mm, dd, yyyy, end_hour, end_min, end_second))
caldat, xtick_labels_temp, month, day, year, hours, minutes, seconds
month = strtrim(string(month), 2)
day = strtrim(string(day),2)
year = strtrim(string(year),2)
hours = strtrim(string(hours),2)
minutes = strtrim(string(minutes),2)
seconds = strtrim(string(seconds),2)
counts = size(xtick_labels_temp, /n_elements)
xtick_labels_all = strarr(counts)
for x_index = 0, counts-1 do begin
if month(x_index) lt 10 then begin
month(x_index) = '0' + month(x_index) ;to make the label "09:00" rather than "9:00"
endif
if day(x_index) lt 10 then begin
day(x_index) = '0' + day(x_index)
endif
if hours(x_index) lt 10 then begin
hours(x_index) = '0' + hours(x_index)
endif
if minutes(x_index) lt 10 then begin
minutes(x_index) = '0' + minutes(x_index)
endif
if seconds(x_index) lt 10 then begin
seconds(x_index) = '0' + seconds(x_index)
endif
xtick_labels_all(x_index) = hours(x_index) + ':' + minutes(x_index)
endfor
xtick_count = 6
xminor_count = 10
jump = 10
xtick_labels = strarr(xtick_count+1)
for x_index = 0, xtick_count do begin
xtick_labels(x_index) = xtick_labels_all(x_index*jump)
end
endif
if duration_sec eq 1800 then begin
xtick_labels_temp = timegen(units='minutes', $
start = julday(mm, dd, yyyy, start_hour, start_min, start_second),$
final = julday(mm, dd, yyyy, end_hour, end_min, end_second))
caldat, xtick_labels_temp, month, day, year, hours, minutes, seconds
month = strtrim(string(month), 2)
day = strtrim(string(day),2)
year = strtrim(string(year),2)
hours = strtrim(string(hours),2)
minutes = strtrim(string(minutes),2)
seconds = strtrim(string(seconds),2)
counts = size(xtick_labels_temp, /n_elements)
xtick_labels_all = strarr(counts+1)
for x_index = 0, counts-1 do begin
if month(x_index) lt 10 then begin
month(x_index) = '0' + month(x_index) ;to make the label "09:00" rather than "9:00"
endif
if day(x_index) lt 10 then begin
day(x_index) = '0' + day(x_index)
endif
if hours(x_index) lt 10 then begin
hours(x_index) = '0' + hours(x_index)
endif
if minutes(x_index) lt 10 then begin
minutes(x_index) = '0' + minutes(x_index)
endif
if seconds(x_index) lt 10 then begin
seconds(x_index) = '0' + seconds(x_index)
endif
xtick_labels_all(x_index) = hours(x_index) + ':' + minutes(x_index)
endfor
if counts eq 30 then begin
hours_added = strtrim(string(hours(counts-1) + 1),2)
if hours_added lt 10 then begin
hours_added = '0' + hours_added
endif
minutes_added = '00'
xtick_labels_all(counts) = hours_added +':' + minutes_added
endif
xtick_count = 3
xminor_count = 10
jump = 10
xtick_labels = strarr(xtick_count+1)
for x_index = 0, xtick_count do begin
xtick_labels(x_index) = xtick_labels_all(x_index*jump)
end
endif
;====== 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_POK = database_dir + '/POK_' + 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_POK = dblarr(3, datacount)
data_array_IQA = dblarr(3, datacount)
data_array_ISR = dblarr(3, datacount)
;====== 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
nsum_value = duration * 5 ;for reduction of data size of time series plot
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!!!!!!!'
ncols = ncols/skip
hzperbin = 1/ (nfft * period)
nrows = round(maxfreq / hzperbin)
maxfreq_calc = hzperbin * nrows
spectra = dblarr(ncols,nrows)
;====== Read in data (stored in database) from the stations selected to be plotted
no_stations = NAL + LYR + HOR + SDY + POK + IQA + ISR
plot_array_Bx = dblarr(no_stations*2, datacount)
plot_array_By = dblarr(no_stations*2, datacount)
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)
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
file_array_NAL = readmud(database_path_NAL, fromSecond, toSecond, sample_rate)
temp_array_NAL = reform(file_array_NAL(0,*))
file_array_NAL(0,*) = convol(temp_array_NAL, coeff) ;digital filtering
temp_array_NAL = reform(file_array_NAL(1,*))
file_array_NAL(1,*) = convol(temp_array_NAL, coeff) ;digital filtering
plot_array_Bx(station_index,*) = times
plot_array_Bx(station_index+1, *) = file_array_NAL(0, *)/system_gain ;x signal in terms of nT
plot_array_By(station_index,*) = times
plot_array_By(station_index+1, *) = file_array_NAL(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.
station_code(station_index/2) = 'NAL'
station_index = station_index + 2
endif
if LYR eq 1 then begin
file_array_LYR = readmud(database_path_LYR, fromSecond, toSecond, sample_rate)
temp_array_LYR = reform(file_array_LYR(0,*))
file_array_LYR(0,*) = convol(temp_array_LYR, coeff) ;digital filtering
temp_array_LYR = reform(file_array_LYR(1,*))
file_array_LYR(1,*) = convol(temp_array_LYR, coeff) ;digital filtering
plot_array_Bx(station_index,*) = times
plot_array_Bx(station_index+1, *) = file_array_LYR(0, *)/system_gain ;x signal in terms of nT
plot_array_By(station_index,*) = times
plot_array_By(station_index+1, *) = file_array_LYR(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.
station_code(station_index/2) = 'LYR'
station_index = station_index + 2
endif
if HOR eq 1 then begin
file_array_HOR = readmud(database_path_HOR, fromSecond, toSecond, sample_rate)
temp_array_HOR = reform(file_array_HOR(0,*))
file_array_HOR(0,*) = convol(temp_array_HOR, coeff) ;digital filtering
temp_array_HOR = reform(file_array_HOR(1,*))
file_array_HOR(1,*) = convol(temp_array_HOR, coeff) ;digital filtering
plot_array_Bx(station_index,*) = times
plot_array_Bx(station_index+1, *) = file_array_HOR(0, *)/system_gain ;x signal in terms of nT
plot_array_By(station_index,*) = times
plot_array_By(station_index+1, *) = file_array_HOR(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.
station_code(station_index/2) = 'HOR'
station_index = station_index + 2
endif
if SDY eq 1 then begin
file_array_SDY = readmud(database_path_SDY, fromSecond, toSecond, sample_rate)
temp_array_SDY = reform(file_array_SDY(0,*))
file_array_SDY(0,*) = convol(temp_array_SDY, coeff) ;digital filtering
temp_array_SDY = reform(file_array_SDY(1,*))
file_array_SDY(1,*) = convol(temp_array_SDY, coeff) ;digital filtering
plot_array_Bx(station_index,*) = times
plot_array_Bx(station_index+1, *) = file_array_SDY(0, *)/system_gain ;x signal in terms of nT
plot_array_By(station_index,*) = times
plot_array_By(station_index+1, *) = file_array_SDY(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.
station_code(station_index/2) = 'SDY'
station_index = station_index + 2
endif
if POK eq 1 then begin
file_array_POK = readmud(database_path_POK, fromSecond, toSecond, sample_rate)
temp_array_POK = reform(file_array_POK(0,*))
file_array_POK(0,*) = convol(temp_array_POK, coeff) ;digital filtering
temp_array_POK = reform(file_array_POK(1,*))
file_array_POK(1,*) = convol(temp_array_POK, coeff) ;digital filtering
plot_array_Bx(station_index,*) = times
plot_array_Bx(station_index+1, *) = file_array_POK(0, *)/system_gain ;x signal in terms of nT
plot_array_By(station_index,*) = times
plot_array_By(station_index+1, *) = file_array_POK(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.
station_code(station_index/2) = 'POK'
station_index = station_index + 2
endif
if IQA eq 1 then begin
file_array_IQA = readmud(database_path_IQA, fromSecond, toSecond, sample_rate)
temp_array_IQA = reform(file_array_IQA(0,*))
file_array_IQA(0,*) = convol(temp_array_IQA, coeff) ;digital filtering
temp_array_IQA = reform(file_array_IQA(1,*))
file_array_IQA(1,*) = convol(temp_array_IQA, coeff) ;digital filtering
plot_array_Bx(station_index,*) = times
plot_array_Bx(station_index+1, *) = file_array_IQA(0, *)/system_gain ;x signal in terms of nT
plot_array_By(station_index,*) = times
plot_array_By(station_index+1, *) = file_array_IQA(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.
station_code(station_index/2) = 'IQA'
station_index = station_index + 2
endif
if ISR eq 1 then begin
file_array_ISR = readmud(database_path_ISR, fromSecond, toSecond, sample_rate)
temp_array_ISR = reform(file_array_ISR(0,*))
file_array_ISR(0,*) = convol(temp_array_ISR, coeff) ;digital filtering
temp_array_ISR = reform(file_array_ISR(1,*))
file_array_ISR(1,*) = convol(temp_array_ISR, coeff) ;digital filtering
plot_array_Bx(station_index,*) = times
plot_array_Bx(station_index+1, *) = file_array_ISR(0, *)/system_gain ;x signal in terms of nT
plot_array_By(station_index,*) = times
plot_array_By(station_index+1, *) = file_array_ISR(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.
station_code(station_index/2) = 'ISR'
station_index = station_index + 2
endif
;====== Specify graphic format
file_format = file_format ;'WIN': Microsoft Windows
;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
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
;====== Set up plot positions
position1 = [x0,y0,x1,y1]
position2 = [x2,y2,x3,y3]
position3 = [x4,y4,x5,y5]
position4 = [x6,y6,x7,y7]
doy = 100 ; fixed (temp)
winsize = convert_coord(xwinsize, dataysize, /normal, /to_device) ;Convert coordinates for the window size from normal to device coords
;print,'Step size = ',(duration/ncols) * spm * mph,' sec
;print,'nrows = ', nrows, ' , ncols = ', ncols
;print,'f_max = ', maxfreq_calc,' Hz'
;====== 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
;====== Plot graphs (time-series and spectrogram; one stations (4 panels) per page)
page_index = station_index/2 ;determine number of page to be plotted in PS format. Prints data from up to 3 stations in one page.
spec_array = dblarr(2, datacount)
for plot_index = 0, station_index-1, 2 do begin
;====== Plot time-series (Bx)
plot, plot_array_Bx(plot_index, *), plot_array_Bx(plot_index+1,*), nsum = nsum_value, pos=position1, $
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 = strtrim(station_code(plot_index/2),2) + ' Bx dB/dt (nT/sec)'
;====== Plot FFT spectrogram (Bx)
plot, [0,1], [0,1], /noerase, /nodata, pos=position2,$
xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
xticks=xtick_count, xminor=xminor_count, 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),2) + ' Bx Freq (Hz)'
spec_array(0, *) = plot_array_Bx(plot_index, *)
spec_array(1, *) = plot_array_Bx(plot_index+1, *)
calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
greyplot, spectra(*,*), position2, winsize, doy, times, pstart, duration, missing, upper, lower
;====== Plot time-series (By)
plot, plot_array_By(plot_index, *), plot_array_By(plot_index+1,*), nsum = nsum_value, pos=position3, /noerase, $
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 = strtrim(station_code(plot_index/2),2) + ' By dB/dt (nT/sec)'
;====== Plot FFT spectrogram (By)
plot, [0,1], [0,1], /noerase, /nodata, pos=position4, $
xstyle=1, xticklen=-0.05, xcharsize=charsize1, $
xticks=xtick_count, xminor=xminor_count, 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),2) + ' By Freq (Hz)'
spec_array(0, *) = plot_array_By(plot_index, *)
spec_array(1, *) = plot_array_By(plot_index+1, *)
calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing
greyplot, spectra(*,*), position4, winsize, doy, times, pstart, duration, missing, upper, lower
;======Print main titles
main_title1 = yyyy +'/'+mm+'/'+dd+' '+start_hour+':'+start_min+'-'+end_hour+':'+end_min+' UT'
xyouts, 0.45, .92, main_title1, /normal, alignment=0.5, charsize=charsize2
main_title2 = 'Page ' + strtrim(string(plot_index/2+1),2) + '/' + strtrim(string(page_index),2)
xyouts, 0.82, .92, main_title2, /normal, alignment=0.5, charsize=charsize2
x_title = 'Universal Time (hh:mm)'
xyouts, 0.45, y6-0.05, x_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