;+
;procedure: thm_cal_fgm
;
;Purpose: applies calibration to THEMIS fluxgate magnetometer data
;
; searches for the current calibration in ASCII file
; takes calibration that has time stamp before or at the start of the data
; converts data to nano Tessla
; applies the calibration
; corrects for phase shift from averaging before despinning
;
;keywords:
; probe = Probe name. The default is 'all', i.e., calibrate data for all
; available probes.
; This can be an array of strings, e.g., ['a', 'b'] or a
; single string delimited by spaces, e.g., 'a b'
; datatype = The type of data to be loaded, 'fge', 'fgh', or 'fgl'. default is
; 'all', to calibrate all variables.
; in_suffix = optional suffix to add to name of input data quantity, which
; is generated from probe and datatype keywords.
; out_suffix = optional suffix to add to name for output tplot quantity,
; which is generated from probe and datatype keywords.
; coord = apply thm_cotrans if output other
; than spinning spacecraft (SSL) is desired.
;
;optional parameters:
;
; name_thx_fgx_in --> input data (t-plot variable name)
; name_thx_fgx_hed --> header information for input data (t-plot variable name)
; name_thx_fgx_out --> name for output (t-plot variable name)
; pathfile --> path and filename of the calibration file
;
;keywords:
;Example:
; tha_cal_fgm, probe = 'a', datatype= 'fgl'
;
;Notes: under construction!!
;
;Written by Hannes Schwarzl.
; $LastChangedBy: jimm $
; $LastChangedDate: 2007-09-27 11:40:03 -0700 (Thu, 27 Sep 2007) $
; $LastChangedRevision: 1634 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/thmsoc/tags/tdas_2_02/idl/themis/spacecraft/fields/thm_cal_fgm.pro $
;Changes by Edita Georgescu
;eg 6/3/2007 - matrix multiplication
;eg 30/3/2007 - keyword deleted, fgmtype determined from: name_thx_fgx_in
; - calibrates fgl,fgh,fge
; - calibrate only mode=0 data
; 2007-03-30 Modified by Hannnes S to use multiple calibrations for a selected interval
; and use spinperiod from the statefile
; Search this string to see changes ;Hannes 30/3/2007
; 2007-03-30 Modified by Hannnes S correction for filter delay removed
; Search this string to see changes: Hannes 05/21/2007
;-
pro thm_cal_fgm,probe=probe, datatype=datatype, in_suffix=in_suf, out_suffix=out_suf, coord=coord, name_thx_fgx_in,name_thx_fgx_hed,name_thx_fgx_out,pathfile ;eg 30/3/2007;Hannes 30/3/2007 ;kb
;; 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']
vdatatypes = ['fgl', 'fgh', 'fge']
if keyword_set(valid_names) then begin
probe = vprobes
datatypes = vdatatypes
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 = vdatatypes $
else dts = thm_check_valid_name(strlowcase(datatype), vdatatypes, $
/include_all)
if not keyword_set(dts) then return
if not keyword_set(in_suf) then in_suf = ''
if not keyword_set(out_suf) then out_suf = ''
for i = 0, n_elements(probes)-1 do begin
thx = 'th' + probes[i]
cal_relpathname = thx+'/l1/fgm/0000/'+thx+'_fgmcal.txt'
cal_file = file_retrieve(cal_relpathname, _extra=!themis)
for j = 0, n_elements(dts)-1 do begin
raw_name = 'th'+probes[i]+'_'+dts[j]+in_suf
cal_name = 'th'+probes[i]+'_'+dts[j]+out_suf
hed_name = 'th'+probes[i]+'_'+dts[j]+'_hed'
thm_cal_fgm, raw_name, hed_name, cal_name, cal_file, coord=coord
endfor
endfor
return
endif else if n_params() ne 4 then begin
print, 'usage: thm_cal_fgm, probe=probe, datatype=datatype $'
print, ' [, in_suffix=in_suf, out_suffix=out_suf]'
print, 'or: thm_cal_fgm, raw_name, hed_name, cal_name, cal_file]'
return
endif
; what follows is exactly the original thm_cal_fgm, before 'standardization'
; (plus some metadata in dlimit.data_att) ;kb
; get the data using t-plot names
get_data,name_thx_fgx_in,data=thx_fgx_in,dl=thx_fgx_dl ;kb
get_data,name_thx_fgx_hed,data=thx_fgx_hed
; kb
if size(thx_fgx_in,/type) ne 8 then begin
print, name_thx_fgx_in, " does not exist."
return
endif
if size(thx_fgx_hed,/type) ne 8 then begin
print, name_thx_fgx_hed, " does not exist: "
message, " Support data not found."
endif
if thm_data_calibrated(thx_fgx_dl) then begin
message, name_thx_fgx_in + " has already been calibrated."
endif
; end kb
; make STRUCTRE thx_fgx; converting egineering units from LONG to DOUBLE
thx_fgx = CREATE_STRUCT('X',thx_fgx_in.X , 'Y', DOUBLE(thx_fgx_in.Y)) ;, 'V', thx_fgx_in.V) -mattd,5/2
; scale factor
kr=2.980232238769531E-3 ;kr=25000.0/2^23 --> see *CALPROC*.doc
; fillvalue for output vector
;fillvalue=1.e-31
fillvalue=!values.f_nan ;Hannes 30/3/2007
; check input ;eg 30/3/2007
fgt=['fgl','fgh','fge']
fgtype=strmid(name_thx_fgx_in,strpos(name_thx_fgx_in,'_fg')+1,3)
tm=where(fgtype eq fgt) & print,name_thx_fgx_in
;eg 30/3/2007
; get array sizes
count=SIZE(thx_fgx.X,/N_ELEMENTS)
countHed=SIZE(thx_fgx_HED.X,/N_ELEMENTS)
PRINT,'get header information ...'
; get the information from the header
fgmhed=thx_fgx_HED.Y[*,12:15] ; store relevant header information
PRINT,'get sampling frequencies ...'
; get fs [vector] ;eg 30/3/2007
case tm of
0 : fs=float(2^(2+(fgmhed[*,3] and '111'B)))
1 : fs=replicate(128.0,countHed)
2 : fs=replicate(8.0,countHed)
endcase
PRINT,'get filter-modes ...'
; get filtermode [vector]
if tm lt 2 then ff=ishft('F0'X and fgmhed[*,1],-6) else ff=replicate(0,countHed)
; get mode [vector] only mode=0 is science data and have to be calibrated ;eg 30/3/2007
if tm lt 2 then mode=ishft(fgmhed[*,1] and '1100'B,-2) else mode=replicate(0,countHed)
;eg 30/3/2007
PRINT,'done getting header information'
;read the calibration file
PRINT,'read calibration file:'
PRINT,pathfile
ncal=file_lines(pathfile)
;define variables
calstr=strarr(ncal)
spinperi=dblarr(ncal)
offi=dblarr(ncal,3)
cali=dblarr(ncal,9)
spinperii=dblarr(1)
offii=dblarr(3)
calii=dblarr(9)
utci='2006-01-01T00:00:00.000Z'
utc=dblarr(ncal)
utcStr=strarr(ncal)
openr,2,pathfile
readf,2,calstr
close,2
for i=0,ncal-1 DO BEGIN
calstri=calstr(i)
utci=strmid(calstr(i),0,25) ;eg 6/3/2007
reads,strmid(calstr(i),26),offii,calii,spinperii;,format="(a25,3f5,9f9,f9)" ;eg 6/3/2007
offi(i,*)=offii
cali(i,*)=calii
spinperi(i)=spinperii
utcStr(i)=utci
;translate time information
STRPUT, utci, '/', 10
utc(i)=time_double(utci)
ENDFOR
PRINT,'done reading calibration file'
;Hannes 30/3/2007
;get the spinperiod from the state file
preSpinper=strmid(name_thx_fgx_in,0,4)
name_thx_spinper=preSpinper+'state_spinper'
PRINT,'take spin period from state file ...'
get_data,name_thx_spinper,data=thx_spinper
if size(thx_spinper,/type) ne 8 then begin
message, 'aborted: must load spin data from state file'
endif
thx_spinper_interp=thm_interpolate_state(thx_xxx_in=thx_fgx_in,thx_spinper=thx_spinper) ;--> linear interpolation
;end Hannes 30/3/2007
PRINT,'search calibration for selected time interval ...'
calIndex=0
compTime=utc(0)
refTime=thx_fgx.X(0)
i=0
WHILE ((compTime lt reftime) && (i lt ncal-1)) DO BEGIN
i=i+1
compTime=utc(i)
IF (compTime gt reftime) THEN BEGIN
i=i-1
BREAK
ENDIF
ENDWHILE
;change 2007-03-29
istart=i
compTime=utc(i)
refTime=thx_fgx.X(count-1L)
;i=0
WHILE ((compTime lt reftime) && (i lt ncal-1)) DO BEGIN
i=i+1
compTime=utc(i)
IF (compTime gt reftime) THEN BEGIN
i=i-1
BREAK
ENDIF
ENDWHILE
istop=i
PRINT, 'Select calibrations from:'
FOR i=istart,istop DO BEGIN
PRINT, utcStr(i)
ENDFOR
;PRINT,'offsets: ',TRANSPOSE(offs)
;PRINT,'cal matrix: '
;PRINT,calm
;PRINT,'spin period: ',spinper
;PRINT, 'Select calibration from:'
;PRINT, utcStr(i)
;offs=offi(i,*)
;calim=TRANSPOSE(cali(i,*))
;calm=[[calim[0:2]],[calim[3:5]],[calim[6:8]]]
;spinper=spinperi(i) Hannes 30/3/2007
;if abs(spinper-mean(tspin.y))/spinper gt 1.e-4 then spinper=mean(tspin.y) ;eg 30/3/2007,Hannes 30/3/2007
;PRINT,'offsets: ',TRANSPOSE(offs)
;PRINT,'cal matrix: '
;PRINT,calm
;PRINT,'spin period: ',spinper Hannes 30/3/2007
;end Hannes 30/3/2007
PRINT,'apply scale factor ...'
;apply scale factor
thx_fgx.Y=thx_fgx.Y*kr
PRINT,'apply calibration ...'
;apply calibration matrix and offset
;lastfs=0
i=0L
ffi=ff[i]
fsi=fs[i]
modi=mode(i) ;eg 30/3/2007
timei=thx_fgx_HED.X[i]
;Hannes 30/3/2007
;starting calibration
offs=offi(istart,*)
calim=TRANSPOSE(cali(istart,*))
calm=[[calim[0:2]],[calim[3:5]],[calim[6:8]]]
;spinper=spinperi(istart)
iact=istart
IF (iact eq istop) THEN BEGIN
nextCalChangeTime=thx_fgx.X(count-1L)+10.d0 ;something greater than the last point in time
ENDIF ELSE BEGIN
nextCalChangeTime=utc(iact+1L)
PRINT, 'nextCalChangeTime:'
PRINT, utcStr(iact+1L)
ENDELSE
;end Hannes 30/3/2007
for j=0L,count-1L do begin
;Hannes 30/3/2007
;apply the right calibration at the right time
IF (thx_fgx.X[j] ge nextCalChangeTime) THEN BEGIN
iact=iact+1
IF (iact eq istop) THEN BEGIN
nextCalChangeTime=thx_fgx.X(count-1L)+10.d0 ;something greater than the last point in time
ENDIF ELSE BEGIN
nextCalChangeTime=utc(iact+1L)
PRINT, 'nextCalChangeTime:'
PRINT, utcStr(iact+1L)
ENDELSE
offs=offi(iact,*)
calim=TRANSPOSE(cali(iact,*))
calm=[[calim[0:2]],[calim[3:5]],[calim[6:8]]]
;spinper=spinperi(iact)
ENDIF
;end Hannes 30/3/2007
; update header information if necessary (WHILE loop should work for all possible cases)
while ((thx_fgx.X[j] gt timei) && (i lt countHed-1)) do begin
i=i+1
ffi=ff[i]
fsi=fs[i]
modi=mode(i)
timei=thx_fgx_HED.X[i]
endwhile
if modi eq 0 then begin ;eg 30/3/2007
; thx_fgx.Y[j,*]=MATRIX_MULTIPLY(calm, thx_fgx.Y[j,*], /BTRANSPOSE)-offs ;eg 6/3/2007
thx_fgx.Y[j,*]=calm ## thx_fgx.Y[j,*] - offs ;eg 6/3/2007
;correct for filter
if ffi eq 2 then begin
spinper=thx_spinper_interp.Y[j];Hannes 30/3/2007
arg=-!dpi/fsi/spinper
dfilt=128./fsi*sin(!dpi/128.d0/spinper)/sin(-arg)
dfm=identity(3,/double)
dfm[1,1]=dfilt
dfm[0,0]=dfilt
;sarg=sin(arg) /delay removed Hannes 05/21/2007
;carg=cos(arg) ;& print,arg,dfilt /delay removed Hannes 05/21/2007
;delaym=identity(3,/double) /delay removed Hannes 05/21/2007
;delaym[0,1]=-sarg ;eg 6/3/2007 /delay removed Hannes 05/21/2007
;delaym[1,0]=sarg ;eg 6/3/2007 /delay removed Hannes 05/21/2007
;delaym[0,0]=carg /delay removed Hannes 05/21/2007
;delaym[1,1]=carg /delay removed Hannes 05/21/2007
;mfilt=dfm#delaym /delay removed Hannes 05/21/2007
mfilt=dfm ;/apply only amplitude correction Hannes 05/21/2007
; thx_fgx.Y[j,*]=MATRIX_MULTIPLY(mfilt, thx_fgx.Y[j,*],/BTRANSPOSE) ;eg 6/3/2007
thx_fgx.Y[j,*]=mfilt ## thx_fgx.Y[j,*] ;eg 6/3/2007
endif ;correct for filter
endif else thx_fgx.Y[j,*]=fillvalue ;eg 30/3/2007 set vectors to fillvalue ???
endfor
PRINT,'calibration applied'
PRINT,'done'
;store data, where mode = 0
; kb
units = 'nT'
dl = thx_fgx_dl
str_element, dl, 'data_att', data_att, success=has_data_att
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', 'ssl', /add
str_element, data_att, 'units', units, /add
str_element, dl, 'data_att', data_att, /add
str_element, dl, 'ytitle', string( name_thx_fgx_out, units, format='(A,"!C!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
; end kb
store_data,name_thx_fgx_out,data=thx_fgx,dl=dl
; kb
if keyword_set(coord) && strlowcase(coord) ne 'ssl' then begin
thm_cotrans, name_thx_fgx_out, out_coord = coord
options, tplot_var, 'ytitle', /def, $
string( name_thx_fgx_out, units, format='(A,"!C!C[",A,"]")'), /add
endif
; end kb
;return,thx_fgx
end