;+ ;NAME: ; thm_cal_scm ;PURPOSE: ; calibrate THEMIS SCM data ; ; The processing is done in 8 successive steps ; ----------- by default, processing stops after step 5. ; ; # 0: counts, NaN inserted into each gap for proper tplotting ; # 1: Volts, spinning sensor system, with DC field ; # 2: Volts, spinning sensor system, without DC field[, xy DC field in nT] ; # 3: nTesla, spinning sensor system, without DC field ; # 4: nTesla, spinning SSL system, without DC field ; # 5: nTesla, fixed DSL system, without DC field, filtered input data (t-plot variable name) ; thx_scx_hed --> header information for input data (t-plot variable name) ; ;Example: ; tha_cal_fgm, probe = 'a', datatype= 'scp' ; ;HISTORY: ; 13-mar-2007, jmm, jimm@ssl.berkeley.edu ; June-2007, krb, Based on Patrick Robert's cowave_THEscmwf.f ;$LastChangedBy: jimm $ ;$LastChangedDate: 2007-07-20 13:04:26 -0700 (Fri, 20 Jul 2007) $ ;$LastChangedRevision: 1193 $ ;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/thmsoc/tags/tdas_2_02/idl/themis/spacecraft/fields/thm_cal_scm.pro $ ;- pro thm_scm_sensor_to_SSL, xfo, yfo, zfo, nbp r2 = xfo*xfo + yfo*yfo + zfo*zfo x= -0.537691*xfo + 0.843131*yfo +0.004316*zfo y= -0.843039*xfo - 0.537696*yfo +0.012899*zfo z= 0.013196*xfo + 0.003295*yfo +0.999913*zfo p2 = x*x + y*y + z*z ;; test if no modulo difference diff = sqrt(p2) - sqrt(r2) epsi = (sqrt(p2)+sqrt(r2))/2.e4 bad = where(abs(diff) gt epsi, nbad) if nbad gt 0 then $ print, '*** sensor to SSL : modulo has changed, diff= ', max(abs(diff)) xfo = x yfo = y zfo = z end Pro thm_cal_scm, probe=probe, datatype=datatype, $ in_suffix=in_suffix, out_suffix=out_suffix, $ trange=trange, $ ;time_range nk = k_nk, $ ;N points of the cal-convolution kernel mk = k_mk, $ ; set nk to multiple of sample frequency fdet = k_fdet, $ ;detrend frequency despin = k_despin, $ ; classic despin 0:off 1:on n_spinfit = k_n_spinfit, $ ;; n spins to fit for misalignment, ;; dc field calculation and despin. fcut = k_fcut, $ ;Frequency cut-off for calibration fmin = k_fmin, $ ;Min frequency for filtering fmax = k_fmax, $ ;Max frequency for filtering step = k_step, $ ;Highest Processing step to complete dircal = k_dircal, $ ;directory of calibration files calfile = k_calfile, $ ;override automatic cal file name fsamp = fsamp, $ ; sample rate. output. edge_truncate=edget, edge_wrap=edgew, edge_zero=edgez, $ no_download=no_download, $ ; don't look on web for cal file coord = k_coord, verbose=verbose, valid_names=valid_names, $ progobj = progobj, $ ;20-jul-2007, jmm thx_scx, thx_scx_hed thm_init ; If verbose keyword is defined, override !themis.verbose vb = size(verbose, /type) ne 0 ? verbose : !themis.verbose ;; implement the 'standard' interface as a wrapper around the original ;; interface if no positional parameters are present. if n_params() eq 0 then begin vprobes = ['a','b','c','d','e'] ; Valid datatypes include the raw data in Volts, dc field, axis misalignment, ; and iano (data quality flag). Created with '_dc', ; '_misalign', and '_iano' suffixes, respectively. vdatatypes = ['scf', 'scp', 'scw', $ 'scf_dc', 'scf_misalign', 'scf_iano', $ 'scp_dc', 'scp_misalign', 'scp_iano', $ 'scw_dc', 'scw_misalign', 'scw_iano'] defdatatypes = ['scf', 'scp', 'scw'] if keyword_set(valid_names) then begin probe = vprobes datatypes = vdatatypes thm_cotrans, out_coord=k_coord, /valid_names, verbose=0 if keyword_set(vb) then begin message, /info, string(strjoin(vprobes, ','), $ format = '( "Valid probes:",X,A,".")') message, /info, string(strjoin(vdatatypes, ','), $ format = '( "Valid datatypes:",X,A,".")') message, /info, string(strjoin(k_coord, ','), $ format = '( "Valid coords:",X,A,".")') endif return endif if not keyword_set(probe) then probes = vprobes $ else probes = thm_check_valid_name(strlowcase(probe), vprobes, /include_all) if not keyword_set(probes) then return if not keyword_set(datatype) then dts = defdatatypes $ else dts = thm_check_valid_name(strlowcase(datatype), vdatatypes, $ /include_all) if not keyword_set(dts) then return for i = 0, n_elements(probes)-1 do begin p_i = probes[i] If(obj_valid(progobj)) Then progobj -> update, 0.0, text = 'Processing Probe: '+ p_i print, 'Processing Probe: ', p_i for j = 0, n_elements(defdatatypes)-1 do begin dt_j = defdatatypes[j] dt_j_types = strfilter(dts,dt_j+'*',count=n_dt_j_types) if n_dt_j_types gt 0 then begin in_name = 'th'+p_i+'_'+dt_j hed_name = 'th'+p_i+'_'+dt_j+'_hed' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = 'Calibrating: '+ dt_j_types print, 'Calibrating: ', dt_j_types thm_cal_scm, in_name, hed_name, in_suffix = in_suffix, $ out_suffix = out_suffix, datatype = dt_j_types, $ trange=trange, $ nk = k_nk, $ mk = k_mk, $ fdet = k_fdet, $ despin = k_despin, $ n_spinfit = k_n_spinfit, $ fcut = k_fcut, $ fmin = k_fmin, $ fmax = k_fmax, $ step = k_step, $ dircal = k_dircal, $ calfile = k_calfile, $ fsamp = fsamp, $ edge_truncate=edget, edge_wrap=edgew, $ edge_zero=edgez, $ no_download=no_download, coord=k_coord, $ verbose=verbose, valid_names=valid_names, $ progobj = progobj endif endfor endfor return endif else if n_params() lt 2 then begin print, 'for usage, type:' print, "doc_library, 'thm_cal_scm'" return endif ; --------------------------------------------------------------------------- ; Here begins the code that gets called if positional parameters are provided ; --------------------------------------------------------------------------- if not keyword_set(in_suffix) then in_suff = '' ELSE in_suff = in_suffix if not keyword_set(out_suffix) then out_suff = '' ELSE out_suff = out_suffix ;instead of reading the rff_header, create the variables using the IDL ;create rff routine, I doubt if any is needed. ;Figure out what probe and mode, from the filename get_data, thx_scx+in_suff, time_scx, data_scx, dlimit = dlim_scx ;unpack the header data get_data, thx_scx_hed, time_scx_hed, val_scx_hed ; test if data is present and in raw form. if size(dlim_scx,/type) ne 8 then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** '+thx_scx+in_suff+' does not exist, skipping.' print, '*** ', thx_scx+in_suff, " does not exist, skipping." return endif if thm_data_calibrated(dlim_scx) then begin in_step = dlim_scx.data_att.step print, '*** Input data already calibrated to step: ', in_step message, '*** Aborted: Currently can only use raw data as input' return endif ; get info needed to unpack the state data temp = file_basename(dlim_scx.cdf.filename) temp = strsplit(temp, '_', /extract) probe = strlowcase(strmid(temp[0], 2, 1)) thx = 'th'+probe probe_n = byte(probe)-96 probe_n = probe_n[0] str_probe_n = string(probe_n, format='(I1)') mode = temp[2] ;unpack the state variables get_data, 'th'+probe+'_state_spinphase', time_thx_state, val_thx_spinpha get_data, 'th'+probe+'_state_spinper', time_thx_state, val_thx_spinper ; get_data, 'th'+probe+'_state_spinras', time_thx_state, val_thx_spinras ; get_data, 'th'+probe+'_state_spindec', time_thx_state, val_thx_spindec ; test if auxiliary data is present if size(val_scx_hed, /n_dim) ne 2 then begin print, '*** ',thx_scx_hed, " does not exist: " message, " Aborted: Support data not found." endif if size(val_thx_spinpha, /n_dim) ne 1 || $ size(val_thx_spinper, /n_dim) ne 1 $ ; l-vector orientation not currently needed - 2007-6-21 ; size(val_thx_spinras, /n_dim) ne 1 || $ ; size(val_thx_spindec, /n_dim) ne 1 $ then begin message, '*** Aborted, must load support data from state file' endif ; ---------------------------------------------------- ; Process calibration parameters, set default vaules. ; ---------------------------------------------------- ;These keywords replace the subroutine r_inpara_rts If size(k_mk, /type) ne 0 Then mk = k_mk Else mk = 4 If size(k_nk, /type) ne 0 Then nk = k_nk Else nk = 0 If size(k_fdet, /type) ne 0 Then fdet = k_fdet Else fdet = 2.0 if size(k_despin, /type) ne 0 Then despin = k_despin Else despin = 1 if size(k_n_spinfit, /type) ne 0 Then $ n_spinfit = k_n_spinfit else n_spinfit = 2 If size(k_fcut, /type) ne 0 Then fcut = k_fcut Else fcut = 0.1 If size(k_fmin, /type) ne 0 Then fmin = k_fmin Else fmin = 0.0 If size(k_fmax, /type) ne 0 Then fmax = k_fmax Else fmax = 0.0 if size(k_step, /type) ne 0 then step = k_step Else step = 5 ;; check value of fcut for lower frequency limit of calibration if fcut lt 0.001 then begin fcut = 0.001 print, 'fcut < 0.001, set to ', fcut endif if keyword_set(k_coord) then begin thm_cotrans, out_coord=vcoord, /valid_names, verbose=0 coord = thm_check_valid_name(strlowcase(k_coord), vcoord) if not keyword_set(coord) then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** invalid coord specification' print, '*** invalid coord specification' return endif if n_elements(coord) ne 1 then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** coord must specify a single coordinate system' print, '*** coord must specify a single coordinate system' return endif if step eq 5 then begin if strlowcase(coord) eq 'ssl' then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** Warning: step 5 requested, with coord=SSL, setting step = 4' print, '*** Warning: step 5 requested, with coord=SSL, setting step = 4' step = 4 endif else if strlowcase(coord) ne 'dsl' && fmin lt fcut then begin fminnew=fcut print If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** Warning: for step 5 output in coord sys other than DSL, ' print, '*** Warning: for step 5 output in coord sys other than DSL, ' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'fmin must be >fcut, inital fmin = '+strcompress(/remove_all, string(fmin))+$ ' set to fcut = '+strcompress(/remove_all, string(fminnew)) print, ' fmin must be >fcut' print, ' inital fmin=',fmin, ' set to fcut =',fminnew fmin=fminnew endif endif endif ; * warning for step > 6 if step gt 6 then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** Warning : step requested = '+strcompress(/remove_all, string(step))+$ ' Step 6 is maximum allowed, set to 6' print, '*** Warning : step requested =', step print, ' Step 6 is maximum allowed, set to 6' print, ' use thm_cotrans to get output in GSE or GSM' step = 6 endif print, 'thm_cal_scm parameters:' print, ' nk: ', nk print, ' mk: ', mk print, ' fdet: ', fdet print, ' despin: ', despin print, 'n_spinfit: ', n_spinfit print, ' fcut: ', fcut print, ' fmin: ', fmin print, ' fmax: ', fmax print, ' step: ', step if keyword_set(coord) then $ print, ' coord: ', coord print, '----------------------------------------------------------' ; generate name of calibration file ; the cal files are per spacecraft IF size(/type, k_dircal) EQ 0 Then begin ;; read calibration files from data repository, rather than source cal_relpathname = thx+'/l1/scmcal/THEMIS_SCM'+str_probe_n+'.cal' calfile = file_retrieve(cal_relpathname, _extra=!themis, $ no_download=no_download) Endif else begin if size(/type, k_dircal) EQ 7 and keyword_set(k_dircal) then begin dircal = k_dircal endif else begin dirscm = routine_info('thm_cal_scm',/source) dirscm = file_dirname(dirscm.path) PRINT, ' dirscm: ',dirscm dircal= filepath(root_dir=dirscm,'scm/cal') ; cal_relpathname = thx'+thx+'_fgmcal.txt' endelse PRINT, ' dircal: ',dircal calfile = 'THEMIS_SCM'+str_probe_n+'.cal' calfile = filepath(root_dir=dircal, calfile) endelse PRINT, ' calfile: ',calfile ; initialize gainant ( ensures that cal file is read once and only once ) void = thm_scm_gainant_vec(/init) ; this code replaces scm_reduce_period_V3 ; no need to reduce state data, so scm_reduce_period{S1,S4} are omitted. if keyword_set(trange) and n_elements(trange) eq 2 then begin t1 = time_double(trange[0]) t2 = time_double(trange[1]) wave_in_range = where(time_scx ge t1 and time_scx le t2, n) if n gt 0 then begin data_scx = data_scx(wave_in_range, *) time_scx = time_scx(wave_in_range) endif else begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'no data in trange'+time_string(trange) print, 'no data in trange', time_string(trange) return endelse endif ;Get SCM sample frequency from the header. This is scm_get_sample_rate. val_scx_apid = reform(uint((32b*val_scx_hed[*,0]/32b))*uint(256) $ + uint(val_scx_hed[*,1])) TMrate = 2.^(reform(val_scx_hed[ *, 14]/16b)+1) ;;dec2hex apid_str = strtrim(string( val_scx_apid(0,0), format='(Z)'), 2) PRINT, 'apid=',apid_str ;; Conversion of TM data in Volts ; do the stuff scm_create_rff does, without writing to a new file: ; interpolation of sample rate aligned on scm time If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Sample rate interpolation...' print, 'Sample rate interpolation...' thm_scm_rate_interpol, time_scx_hed, TMrate, time_scx, SA_rate ; interpolation of spinras, spindec.... not needed ; interpolation of spin phase ; print, 'Phase interpolation...' ; print, ' starting at ', systime(0) ; thm_spin_interpol, time_thx_state, val_thx_spinpha, val_thx_spinper, $ ; time_scx, SP_phase ; print, ' ending at ', systime(0) If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Phase interpolation using sunpulse times...' print, 'Phase interpolation using sunpulse times...' thm_spin_phase, time_scx, SP_phase, probe=probe iano = intarr(n_elements(time_scx)) ;; process continuous stretches of data in separate batches: ;; they require differing kernel sizes, based on sample rate, spin rate ;; --------------------- ;; Step 0 ;; --------------------- ;; find sample rate changes, data gaps and time reversals ;; find data rates in this file. rates = SA_rate[UNIQ(SA_rate, SORT(SA_rate))] if nk EQ 0 then nks = rates*mk ;; nks is the set of various kernel sizes fsamp = rates n_rates = n_elements(rates) If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ ' Sample frequencies found: '+strcompress(/remove_all, string(rates)) print, ' Sample frequencies found: ', rates ;; find data gaps, working with one sample rate at a time. for r = 0, n_rates-1 do begin fe = rates[r] If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ ' Searching for discontinuities in data at samp. freq.: '+$ strcompress(/remove_all, string(fe))+' dt: '+strcompress(/remove_all, string(1.0/fe)) print, ' Searching for discontinuities in data at samp. freq.: ', fe print, ' dt: ', 1.0/fe ind_r = where(SA_rate eq fe) dt = time_scx[ind_r[1:*]]-time_scx[ind_r] reverse = where(dt lt 0, n_reverse) if n_reverse gt 0 then $ iano[ind_r[reverse]]=16 ; time reverse discontinuity = where(abs(dt-1.0/fe) gt 2.0e-5, n_discont) if n_discont gt 0 then $ iano[ind_r[discontinuity]] = 17 ; discontinuity in time ;; discontinuity in sample rate (Minor bug: some changes in sample rate ;; can be mislabled as discontinuities in time) iano[ind_r[n_elements(ind_r)-1]] = 22 endfor ;; Prepare to process data in batches, applying steps 1-6 on each batch coord_sys = 'scs' units = 'counts' errs = where(iano ge 15, n_batches) If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Number of continuous batches: '+strcompress(/remove_all, string(n_batches)) print, 'Number of continuous batches: ', n_batches ;; size of output arrays includes a space for one NaN after each batch nout = n_elements(time_scx)+n_batches out_scx = fltarr(nout,3) out_scx_dc = fltarr(nout,2) out_scx_misalign = fltarr(nout) iano_out=fltarr(nout) data_scx_out = fltarr(nout,3) time_scx_out = dblarr(nout) ind0 = 0 for batch = 0, n_batches-1 do begin ind1 = errs[batch] nbp = ind1-ind0+1 print, '----------------------------------------------------------' print, '----------------------------------------------------------' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Batch #'+strcompress(/remove_all, string(batch))+$ ' Nbp of the batch = '+strcompress(/remove_all, string(nbp))+', duration is '+$ strcompress(/remove_all, string(float(nbp)/fe))+ ' sec.' print, 'Batch #', batch print, 'Nbp of the batch=', nbp, ', so duration is ', float(nbp)/fe, ' sec.' ; ; The processing is done in 8 successive steps ; ; # 1: Volts, spinning sensor system, with DC field ; # 2: Volts, spinning sensor system, without DC field, xy DC field in nT ; # 3: nTesla, spinning sensor system, without DC field ; # 4: nTesla, spinning SSL system, without DC field ; # 5: nTesla, fixed DSL system, without DC field ; # 6: nTesla, fixed DSL system, with xy DC field ; # 7: nTesla, fixed GSE system, without DC field, filtered update, 0.0, text = $ 'step 1: TM data in spinning system [Volt] with DC' print,'step 1: TM data in spinning system [Volt] with DC' print calfactor = (10.04/2.^16) xfo *= calfactor yfo *= calfactor zfo *= calfactor units = 'V' ;;thm_scm_testsatura endif if step ge 2 then begin print,'-----------------------------------------------------' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'step 2: TM data in spinning system [Volt] without DC' print ;; get antenna response at spin frequency thm_scm_modpha, thm_scm_gainant_vec(fs, 1, calfile, fe), rfsx, phafsx thm_scm_modpha, thm_scm_gainant_vec(fs, 2, calfile, fe), rfsy, phafsy thm_scm_modpha, thm_scm_gainant_vec(fs, 3, calfile, fe), rfsz, phafsz print print, 'Spin Frequency: ', fs print, 'Transfer function at spin frequency:' print, 'rfsx (V/nT),phafsx (d.) =', rfsx, phafsx print, 'rfsy (V/nT),phafsy (d.) =', rfsy, phafsy print, 'rfsz (V/nT),phafsz (d.) =', rfsz, phafsz ;; spin phase sph = SP_phase[ind0:ind1] ;; check for abnormally incrementing spin phase print If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'number of spins in this batch= '+strcompress(/remove_all, string(nbp*fs/fe)) print, 'number of spins in this batch= '+strcompress(/remove_all, string(nbp*fs/fe)) dph = SP_phase[ind0+1:ind1]-SP_phase[ind0:ind1-1] dph_neg = where(dph lt 0, n_dph_neg) if n_dph_neg gt 0 then dph[dph_neg]+=360.0 bad_sph = where(abs(dph - 360.0*fs/fe) gt 360.0*fs/fe * 0.01, n_bad_sph) if n_bad_sph gt 0 then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** n bad spin phase values: '+strcompress(/remove_all, string(n_bad_sph))+$ ' out of '+ strcompress(/remove_all, string(nbp)) print, '*** n bad spin phase values: ', n_bad_sph, ' out of ', nbp endif If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Initial centering of the waveforms: ' print, 'Initial centering of the waveforms: ' xavg = total(xfo)/nbp yavg = total(yfo)/nbp zavg = total(zfo)/nbp xfo -= xavg yfo -= yavg zfo -= zavg ;; calculate sine fit for calculation of misalignment angle and ;; classic_despin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Spin fit: batch #'+strcompress(/remove_all, string(batch))+$ ' number of spins = '+strcompress(/remove_all, string(nbp*fs/fe)) print, 'Spin fit: ' thm_scm_casinus_vec, xfo, fe, fs, n_spinfit, $ axvo, phaxvo, n, nbi, xsub thm_scm_casinus_vec, yfo, fe, fs, n_spinfit, $ ayvo, phayvo, n, nbi, ysub thm_scm_casinus_vec, zfo, fe, fs, n_spinfit, $ azvo, phazvo, n, nbi, zsub if despin gt 0 then begin ;; * Remove the sine signal If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Removing sine signal: batch #'+strcompress(/remove_all, string(batch)) print, 'Removing sine signal: ' xfo = xsub yfo = ysub zfo = zsub endif if nbp lt nbi then begin iano[ind0:ind1] = iano[ind1] ; flag all points endif ;; Determine antenna misalignment. ;; this is cal_depoint (vectorized) If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Antenna Misalignment Calculation: batch #'+strcompress(/remove_all, string(batch)) print, 'Antenna Misalignment Calculation: ' ;; compute spin sine in nT bpx=axvo/rfsx bpy=ayvo/rfsy bpz=azvo/rfsz ;; calculate the misalignment of Z axis from spin axis deno = sqrt(bpx*bpx + bpy*bpy + bpz*bpz) null_deno = where(deno le 1.e-10, n_null_deno) if n_null_deno gt 0 then deno[null_deno] = !values.f_nan sd = bpz*sqrt(2.)/deno misalign = sd + !values.f_nan sd_good = where(abs(sd) LE 1., n_sd_good) if n_sd_good gt 0 then $ misalign[sd_good] = asin(sd[sd_good])*180/!dpi ;; if fdet is not zero, then detrend to remove low frequency if fdet gt 0 then begin nsmooth=long(fe/fdet+0.5) If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ ' Sample rate, fdetrend='+strcompress(/remove_all, string(fe))+', '+$ strcompress(/remove_all, string(fdet))+$ ' Number of points for smoothing : ' +strcompress(/remove_all, string(nsmooth)) print, ' Sample rate, fdetrend=',fe,fdet print, ' Number of points for smoothing : ',nsmooth if nsmooth gt nbp then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** batch '+strcompress(/remove_all, string(batch))+': Detrend frequency too small' print, '*** batch ',batch, 'Detrend frequency too small' iano[ind0:ind1] = iano[ind1] ; flag all points endif if nsmooth lt 2 then begin If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** batch '+strcompress(/remove_all, string(batch))+': Detrend frequency too high' print, '*** batch ',batch, 'Detrend frequency too high' iano[ind0:ind1] = iano[ind1] ; flag all points endif xfo -= smooth(xfo, nsmooth, /edge_truncate) yfo -= smooth(yfo, nsmooth, /edge_truncate) zfo -= smooth(zfo, nsmooth, /edge_truncate) endif ;; final centering of waveform xavg = total(xfo)/nbp yavg = total(yfo)/nbp zavg = total(zfo)/nbp xfo -= xavg yfo -= yavg zfo -= zavg ; *** Computation of the X-Y values of the DC field in fixed DSL system ; * X,Y amplitude and phase are computed in sensor system ; * Transformation to SSL is a rotation of about 45 + 12.4 degres ; * and must be consistent with the Sensor_to_SSL matrix ; * Tranformation from SSL to DSL is a rotation of sun pulse phase. print If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'Computation of the DC field in the spin plane:' print, 'Computation of the DC field in the spin plane:' rotdeg= 45. + 12.45 + sph phaxc= -(phaxvo-phafsx-rotdeg) phayc= -(phayvo-phafsy-rotdeg) depha=phayc-phaxc depha = (phayvo-phaxvo) mod 360 gt180 = where(depha gt 180, n_gt180) if n_gt180 gt 0 then depha[gt180] -= 360 lt_180 = where(depha lt -180, n_lt_180) if n_lt_180 gt 0 then depha[lt_180] += 360 ;;print, 'In sensor system, bpx,phaxvo,phafsx =',bpx,phaxvo,phafsx ;;print, 'In sensor system, bpy,phayvo,phafsy =',bpy,phayvo,phafsy ;;print, 'In DSL, after cal, phaxc,phayc,diff=',phaxc,phayc,depha bdcx=bpx*sin((phaxc)*!dpi/180.) bdcy=bpy*sin((phayc)*!dpi/180.) ;;print, 'X,Y of DC field in DSL=',bdcx,bdcy ;; save values of dc field and misalignment for storage in tplot ;; insert a NaN in the gap out_scx_misalign[ind0+batch:ind1+batch] = misalign out_scx_misalign[ind1+batch+1] = !values.f_nan out_scx_dc[ind0+batch:ind1+batch, 0] = bdcx out_scx_dc[ind0+batch:ind1+batch, 1] = bdcy out_scx_dc[ind1+batch+1, *] = !values.f_nan endif if step ge 3 then begin print,'-----------------------------------------------------' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'step 3: Calibrated data in sensor spinning system [nT] without DC' print, 'step 3: Calibrated data in sensor spinning system [nT]'+$ ' without DC' print print, 'Deconvolution-calibration ...' print, 'Sample rate =', fe print, 'Size of FIR kernel =', nk thm_scm_deconvo_vec, xfo, nbp, nk, fcut, fnyq, fe,'thm_scm_gainant_vec',$ 1, calfile, 0., $ edge_t=edget, edge_w=edgew, edge_z=edgez thm_scm_deconvo_vec, yfo, nbp, nk, fcut, fnyq, fe,'thm_scm_gainant_vec',$ 2, calfile, 0., $ edge_t=edget, edge_w=edgew, edge_z=edgez thm_scm_deconvo_vec, zfo, nbp, nk, fcut, fnyq, fe,'thm_scm_gainant_vec',$ 3, calfile, 0., $ edge_t=edget, edge_w=edgew, edge_z=edgez units = 'nT' endif if step ge 4 then begin print,'-----------------------------------------------------' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'step 4: Calibrated data in SSL spinning system [nT] without DC' print,'step 4: Calibrated data in SSL spinning system [nT]'+$ ' without DC' print print, 'Performing rotation: ' thm_scm_sensor_to_SSL, xfo, yfo, zfo, nbp coord_sys = 'ssl' endif if step ge 5 then begin print,'-----------------------------------------------------' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'step 5: Calibrated data in DSL system [nT] without DC '+$ strcompress(/remove_all, string(f1))+' '+$ strcompress(/remove_all, string(f2)) print, 'step 5: Calibrated data in DSL system [nT]'+$ ' without DC', f1, f2 print print, 'Performing rotation in spin plane: ' sinphi = sin(sph/180.0*!dpi) cosphi = cos(sph/180.0*!dpi) xo = cosphi*xfo - sinphi*yfo yo = sinphi*xfo + cosphi*yfo xfo = xo yfo = yo ;; filter in fixed system between f1 and f2 (optional) if (abs(f1) le 1.e-6 && abs(f2-fnyq) le 1.e-6) then begin print, 'no need to apply filtering' endif else begin ;; deconvo w/o gain just does filtering: print, 'Applying filter in DSL system' thm_scm_deconvo_vec, xfo, nbp, nk, f1, f2, fe, '', $ 1, void, 0., $ edge_t = edget, edge_w = edgew, edge_z = edgez thm_scm_deconvo_vec, yfo, nbp, nk, f1, f2, fe, '', $ 2, void, 0., $ edge_t = edget, edge_w = edgew, edge_z = edgez endelse coord_sys = 'dsl' endif if step ge 6 then begin print,'-----------------------------------------------------' If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ 'step 6: Calibrated data in DSL system [nT] with DC '+$ strcompress(/remove_all, string(f1))+' '+$ strcompress(/remove_all, string(f2)) print,'step 6: Calibrated data in DSL system [nT]'+$ ' with DC', f1, f2 if start_step gt 2 then begin print, "*** Must start with step 2 in order to do step 6" endif else begin print, " Adding DC components of B in spin plane." xfo += bdcx yfo += bdcy endelse endif ;; place data into output arrays, placing NaN in gap time_scx_out[ind0+batch:ind1+batch] = time_scx[ind0:ind1] time_scx_out[ind1+batch+1] = time_scx[ind1] + 1.0/fe data_scx_out[ind0+batch:ind1+batch,*] = data_scx[ind0:ind1, *] data_scx_out[ind1+batch+1,*] = !values.f_nan out_scx[ind0+batch:ind1+batch, 0] = xfo out_scx[ind0+batch:ind1+batch, 1] = yfo out_scx[ind0+batch:ind1+batch, 2] = zfo out_scx[ind1+batch+1, *] = !values.f_nan iano_out[ind0+batch:ind1+batch] = iano[ind0:ind1] iano_out[ind1+batch+1] = !values.f_nan ind0 = ind1+1 endfor ;------------------ ; outputs to tplot ;------------------ thx_scx_out = thx_scx+out_suff ; diagnostic outputs: if strfilter(datatype, '*_misalign') ne '' then $ store_data, thx_scx+'_misalign', data={x:time_scx_out, y:out_scx_misalign} if strfilter(datatype, '*_dc') ne '' then $ store_data, thx_scx+'_dc', data={x:time_scx_out, y:out_scx_dc} if strfilter(datatype, '*_iano') ne '' then $ store_data, thx_scx+'_iano', data={x:time_scx_out, y:iano_out} ; store the calibrated output, along with metadata dl = dlim_scx str_element, dl, 'data_att', data_att, success=has_data_att str_Nk = string(Nks ,format='(i6)') str_Despin= string(Despin,format='(i1)') str_N_spinfit= string(N_spinfit,format='(i1)') str_Fdet = string(Fdet ,format='(f7.2)') str_Fcut = string(Fcut ,format='(f6.2)') str_Fmin = string(Fmin ,format='(f7.2)') str_Fmax = string(Fmax ,format='(f7.2)') str_step = string(step ,format='(i1)') str_Fsamp = string(Rates ,format='(f5.0)') str_param = 'Nk='+str_Nk[0]+ $ ', Step='+str_step+' ,Despin='+str_Despin+$ ', N_spinfit='+str_n_spinfit+$ ', Fdet='+str_Fdet+'Hz, Fcut='+str_Fcut+'Hz, Fmin=' $ +str_Fmin+ 'Hz, Fmax='+str_Fmax+'Hz' if has_data_att then begin str_element, data_att, 'data_type', 'calibrated', /add endif else data_att = { data_type: 'calibrated' } str_element, data_att, 'coord_sys', coord_sys, /add str_element, data_att, 'units', units, /add str_element, data_att, 'fsamp', Rates, /add str_element, data_att, 'nk', Nks, /add str_element, data_att, 'despin', Despin, /add str_element, data_att, 'n_spinfit', N_spinfit, /add str_element, data_att, 'fdet', Fdet, /add str_element, data_att, 'fcut', Fcut, /add str_element, data_att, 'fmin', Fmin, /add str_element, data_att, 'fmax', Fmax, /add str_element, data_att, 'step', step, /add str_element, data_att, 'str_cal_param', str_param, /add str_element, dl, 'data_att', data_att, /add str_element, dl, 'ytitle', string( 'th'+probe+' '+mode, str_Fsamp[0], units, $ format='(A,"!C",A,"!C[",A,"]")'), /add str_element, dl, 'labels', [ 'x', 'y', 'z'], /add str_element, dl, 'labflag', 1, /add str_element, dl, 'colors', [ 2, 4, 6], /add str_element, dl, 'color_table', 39, /add store_data, thx_scx_out, data={x:time_scx_out, y:out_scx}, dl=dl if keyword_set(coord) then begin if step eq 6 && strlowcase(coord) ne 'dsl' then begin print If(obj_valid(progobj)) Then progobj -> update, 0.0, text = $ '*** Warning: for step 6 data can only be provided in DSL, coord keyword ignored' print, '*** Warning: for step 6 data can only be provided in DSL, coord' print, ' keyword ignored' endif else begin thm_cotrans,thx_scx_out, out_coord = coord options, thx_scx_out, 'ytitle', /def, $ string( 'th'+probe+' '+mode, str_Fsamp[0], units, $ format='(A,"!C",A,"!C[",A,"]")'), /add endelse endif end