;This routine plots data from SPA and MCM ;Because of timing gaps, daily ASCII files are read ;Modified: ;May 30, 2018 MWC now use ulf_tick_labels to make labels for X axis ;example: ulf_make_plot_sp, database_dir, '20170501', '000000', '240000', plot_dir, 'ps', 'SPA', 1, 2048, 200, 0.2, 0, 1 pro ulf_make_plot_sp, $ database_dir, $ date, $ start_time, $ end_time, $ plot_dir, $ file_format, $ station_in, $ maxfreq, $ nfft, $ step, $ max_B_field, $ f_low, $ f_high, $ upper, $ lower ;;====== The followings will be located in batchprocess.pro database_dir = '/mirl/ULF/incoming/database/' + station_in + '_converted/' plot_dir = '/mirl/ULF/data_plot/' + station_in file_format = 'ps' 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 if station_in eq 'SPA' then begin station_code_one_letter = 'A' endif if station_in eq 'MCM' then begin station_code_one_letter = 'M' endif if station_in eq 'JBS' then begin station_code_one_letter = 'J' 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) ulf_tick_labels, double(start_time_sec), double(end_time_sec), xtick_count, tick_values, xtick_labels ;====== Specify MUD sources database_path = database_dir + '/' + $ strtrim(string(yyyy), 1) + '_' + $ strtrim(string(mm), 1) + '_' + $ strtrim(string(dd), 1) + '_' + $ station_in + '_dbdt_v2.txt' ;====== Initialize data variables sample_rate = 10. ;Hz period = 1./sample_rate ;sec datacount = long(duration * sph * sample_rate);double(float(duration * sph * sample_rate)) times = findgen(datacount) * period/sph ;====== 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 plot_array_Bx = dblarr(datacount) plot_array_By = dblarr(datacount) fromSecond = (long(dd) - 1) * 86400 + (long(start_hour)*3600) + (long(start_min)*60) + start_second toSecond = long(fromSecond) + long(duration_sec) system_gain = 1;4.43 ;system gain = 4.43 V/(nT*Hz) --> applies only to UNH ULF system if station_in eq 'SPA' then $ ulf_read_sp, database_path, data_count, times, file_array, 3 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 ;temp_array = reform(file_array[0,*]) ;file_array[0,*] = convol(temp_array, coeff) ;digital filtering ;temp_array = reform(file_array[1,*]) ;file_array[1,*] = convol(temp_array, coeff) ;digital filtering plot_array_Bx = file_array[0, *] / system_gain ;x signal in terms of nT plot_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. ;====== 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_code_one_letter + graphic_type file_path_plot = plot_dir + '/' + file_name_plot file_path_base = plot_dir + '/' + $ date + '_' + start_time + '_' + end_time +'_' + $ station_code_one_letter ;====== 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., yoffset=0., xsize=8.5, ysize=11;was .75, .5, 7 and 9.5 MWC 20170615 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 ;====== 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 = 1 ;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) ;====== Plot time-series (Bx) plot, times, plot_array_Bx, 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 = station_in + ' 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= station_in + ' Bx Freq (Hz)' spec_array[0, *] = times spec_array[1, *] = plot_array_Bx 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, times, plot_array_By, 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 = station_in + ' 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= station_in + ' By Freq (Hz)' spec_array[0, *] = times spec_array[1, *] = plot_array_By 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(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 ;========================================================================== ; 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