;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: gak_cdf, '20170321', '000000', '240000', plot_dir, 'NAL', 2048, 200, 0.2 pro gak_cdf, $ date, $ start_time, $ end_time, $ plot_dir, $ station_in, $ nfft, $ step, $ max_B_field ;;====== The followings will be located in batchprocess.pro database_dir = '/mirl/ULF/Database' ;plot_dir = '/mirl/ULF/data_plot/' + station_in plot_dir = '/home/mirl/idl/' maxfreq = 2 nfft = 2^11 step = 60 ;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 'NAL' then begin station_code_one_letter = 'N' endif if station_in eq 'LYR' then begin station_code_one_letter = 'L' endif if station_in eq 'HOR' then begin station_code_one_letter = 'H' endif if station_in eq 'SDY' then begin station_code_one_letter = 'S' endif if station_in eq 'POK' then begin station_code_one_letter = 'P' endif if station_in eq 'IQA' then begin station_code_one_letter = 'I' endif if station_in eq 'ISR' then begin station_code_one_letter = 'R' endif if station_in eq 'GAK' then begin station_code_one_letter = 'G' 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. ;====== Specify MUD sources database_path = database_dir + '/' + station_in + '_' + $ strtrim(string(yyyy), 1) + '_' + $ strtrim(string(mm), 1) + '.MUD' ;====== Initialize data variables sample_rate = 20. ;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 = 4.43 ;system gain = 4.43 V/(nT*Hz) --> applies only to UNH ULF system file_array = readmud_gak(database_path, fromSecond, toSecond, sample_rate) 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 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 file_path_plot = plot_dir + '/' + file_name_plot file_path_base = plot_dir + '/' + $ date + '_' + start_time + '_' + end_time +'_' + $ station_code_one_letter spec_array = dblarr(2, datacount) ;set freq in skt file ;for i=0, 204 do print, ' [',strcompress(string(i+1), /remove_all),'] = ', 2*i / 205. + 1./205./2. year_str = strmid(date, 0, 4) month_array = ["Jan","Feb","Mar","Apr","May","Jun",$ "Jul","Aug","Sep","Oct","Nov","Dec"] month_str = month_array[fix(strmid(date, 4, 2))-1] day_str = strmid(date, 6, 2) start_epoch = day_str + '-' + month_str + '-' + year_str + ' 00:00:00.000' start_time_long64 = cdf_parse_tt2000(start_epoch) stop_time_long64 = start_time_long64 + 86400000000000LL skeleton_base = ' -delete ./' cdf_file = 'ulf' param = 'ulf' cmd_base = '/usr/local/cdf/bin/skeletoncdf -cdf ' + cdf_file cmd = cmd_base + skeleton_base + param spawn, cmd spawn, 'chmod 644 ' + cdf_file + '.cdf' cdf_id = cdf_open(cdf_file) spec_times = findgen(ncols) * (duration * spm * mph) / ncols cdf_varput, cdf_id, 'Epoch', $ start_time_long64 + long64(spec_times) * 1000000000LL spec_array[0, *] = times spec_array[1, *] = plot_array_Bx calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing cdf_varput, cdf_id, 'ulf_x', transpose(spectra), /zvariable spec_array[1, *] = plot_array_By calcspec, spec_array, spectra, duration, period, sph, nfft, pstart, missing cdf_varput, cdf_id, 'ulf_y', transpose(spectra), /zvariable cdf_close, cdf_id end