;+ ;***************************************************************************************** ; ; FUNCTION : rbsp_remove_dcfield (OBSOLETE: SEE rbsp_efw_DCfield_removal_crib.pro) ; PURPOSE : Subtract off DC magnetic field from Bw values. Saves as new tplot files ; ; CALLED BY: see rbsp_efw_DCfield_removal_crib ; ; ; REQUIRES: Need: SSL THEMIS software and IDL Geopack DLM both found at ; http://themis.ssl.berkeley.edu/software.shtml ; ; ; INPUT: ; magname -> tplot name of the [n,3] magnetic field data in GSM ; posname -> tplot name of the [n,3] GSM position data ; probe -> 'a' or 'b' ; model -> Any of the Tsyganenko models or IGRF. Defaults to 't96' ; ; ; ; EXAMPLES: see remove_dcfield_crib.pro ; ; Produces the following variables (ex for magname = 'rbspa_mag_gsm', model='t96') ; rbspa_mag_gsm_t96 ; rbspa_mag_gsm_t96_dif ; ; If the model is not t89 then produces: ; rbspa_mag_gsm_t96_wind ; rbspa_mag_gsm_t96_wind_dif ; rbspa_mag_gsm_t96_ace ; rbspa_mag_gsm_t96_ace_dif ; rbspa_mag_gsm_t96_omni ; rbspa_mag_gsm_t96_omni_dif ; ; ; ; NOTES: ruthlessly pilfered from THEMIS crib sheets ; Subtracts off Tsyganenko values using input from ACE, Wind and OMNI, as ; well as user defined input. ; ; ; ; CREATED: 2012/04/13 ; CREATED BY: Aaron W. Breneman ; LAST MODIFIED: 12/10/2012 v1.0.0 - major update. Much simplified ; MODIFIED BY: AWB ; ;***************************************************************************************** ;- pro rbsp_subtract_dcfield,$ magname,$ posname,$ probe,$ model=model if ~keyword_set(model) then model = 't96' model = strlowcase(model) rbspx = 'rbsp'+probe ;In order for the TDAS model subtraction routines to work we have to have the ;mag data named in the following convention: rbsp[a,b]_fgl_gsm get_data,magname,data=dat,dlimits=dlim,limits=lim store_data,rbspx + '_fgl_gsm',data=dat,dlimits=dlim,limits=lim get_data,rbspx+'_fgl_gsm',data=mag_gsm ;GET ACTUAL AND MODEL DATA ;actual kp values can be found at: http://www.ngdc.noaa.gov/stp/GEOMAG/kp_ap.html if model eq 't89' then par = 2.0D ;call the appropriate magnetic field model if model ne 't89' and model ne 'igrf' then begin call_procedure,'t'+model,posname,pdyn=2.0D,dsti=-30.0D,$ yimf=0.0D,zimf=-5.0D endif else begin if model eq 't89' then call_procedure,'t'+model,posname,kp=2.0 if model eq 'igrf' then call_procedure,'tt89',posname,kp=2.0,igrf_only=1 endelse ;Interpolate the model to the number of data points of actual data tinterpol_mxn,posname + '_b'+model,$ mag_gsm.x,$ newname=magname + '_' + model store_data,[posname + '_b'+model],/delete ;now substract model from input data dif_data,magname,$ magname + '_' + model,$ newname=magname + '_' + model + '_dif' ; tplot,[magname,$ ; magname + '_' + model,$ ; magname + '_' + model + '_dif'] ;AUTO PARAMETER DETERMINATION FROM ACTUAL DATA if model ne 't89' then begin ;-------------- ;WIND ;-------------- ;you may have to set the default download directory manually ;here are some examples: ;setenv,'ROOT_DATA_DIR=~/data' ;good for single user unix/linux system ;setenv,'ROOT_DATA_DIR=C:/Documents and Settings/YOURUSERNAME/My Documents' ;example if you don't want to use the default windows location (C:/data/ or E:/data/) kyoto_load_dst ;load wind data wi_mfi_load,tplotnames=tn wi_3dp_load,tplotnames=tn2 if (tn[0] ne '') and (tn2[0] ne '') then begin cotrans,'wi_h0_mfi_B3GSE','wi_b3gsm',/GSE2GSM get_tsy_params,'kyoto_dst','wi_b3gsm','wi_3dp_k0_ion_density','wi_3dp_k0_ion_vel',strupcase(model) ;Call the model with the Wind parameters if model eq 'igrf' then call_procedure,'igrf',posname,parmod=model+'_par'$ else call_procedure,'t'+model,posname,parmod=model+'_par' get_data,posname + '_b'+model,data=dat store_data,posname + '_b'+model+'_wind',data=dat ;Interpolate the model to the number of data points of actual data tinterpol_mxn,posname + '_b'+model+'_wind',$ mag_gsm.x,$ newname=magname + '_'+ model + '_wind' dif_data,magname,$ magname + '_' + model + '_wind',$ newname=magname + '_' + model + '_wind_dif' ; tplot,[magname,$ ; magname + '_' + model + '_wind',$ ; magname + '_' + model + '_wind_dif'] endif else print,'==> NO WIND DATA AVAILABLE' ;----------------------------------- ;ACE (only available from 2011 on) ;----------------------------------- ace_mfi_load,tplotnames=tn ace_swe_load,tplotnames=tn2 if (tn[0] ne '') and (tn2[0] ne '') then begin ;load_ace_mag loads data in gse coords cotrans,'ace_k0_mfi_BGSEc','ace_mag_Bgsm',/GSE2GSM get_tsy_params,'kyoto_dst','ace_mag_Bgsm','ace_k0_swe_Np','ace_k0_swe_Vp',strupcase(model),/speed if model eq 'igrf' then call_procedure,'igrf',posname,parmod=model+'_par' $ else call_procedure,'t'+model,posname,parmod=model+'_par' get_data,posname + '_b'+model,data=dat store_data,posname + '_b'+model+'_ace',data=dat ;Interpolate the model to the number of data points of actual data tinterpol_mxn,posname + '_b'+model+'_ace',$ mag_gsm.x,$ newname=magname + '_' + model + '_ace' dif_data,magname,$ magname + '_' + model + '_ace',$ newname=magname + '_' + model + '_ace_dif' ; tplot,[magname,$ ; magname + '_' + model + '_ace',$ ; magname + '_' + model + '_ace_dif'] endif else print,'==> NO ACE DATA AVAILABLE' ;--------- ;OMNI ;--------- ;omni data example ;NOTE: you may want to degap and deflag the data(using tdegap and tdeflag) ;to remove gaps and flags in the tsyganemo parameter data, especially ;if you find that there are large gaps in the result omni_hro_load,tplotnames=tn if tn[0] ne '' then begin store_data,'omni_imf',data=['OMNI_HRO_1min_BY_GSM','OMNI_HRO_1min_BZ_GSM'] get_tsy_params,'kyoto_dst','omni_imf','OMNI_HRO_1min_proton_density','OMNI_HRO_1min_flow_speed',strupcase(model),/speed,/imf_yz if model eq 'igrf' then call_procedure,'igrf',posname,parmod=model+'_par' $ else call_procedure,'t'+model,posname,parmod=model+'_par' get_data,posname + '_b'+model,data=dat store_data,posname + '_b'+model+'_omni',data=dat ;Interpolate the model to the number of data points of actual data tinterpol_mxn,posname + '_b'+model+'_omni',$ mag_gsm.x,$ newname=magname + '_' + model + '_omni' dif_data,magname,$ magname + '_' + model + '_omni',$ newname=magname + '_' + model + '_omni_dif' ; tplot,[magname,$ ; magname + '_' + model + '_omni',$ ; magname + '_' + model + '_omni_dif'] endif else print,'==> NO OMNI DATA' endif ;dipole tilt example ;add one degree to dipole tilt ;Can also add time varying tilts, or replace the default dipole tilt with a user defined value ; tt96, 'th'+probe+'_state_pos',pdyn=2.0D,dsti=-30.0D,yimf=0.0D,zimf=-5.0D,get_tilt='tilt_vals',add_tilt=1 ; tplot, ['th'+probe+'_state_pos_bt96', 'th'+probe+'_fgs_gsm','tilt_vals'] ;Remove, rename stuff... store_data,['*OMNI_HRO*'],/delete store_data,['*omni_imf*'],/delete store_data,['*ace_k0*','ace_mag_Bgsm'],/delete store_data,['*wi_3dp*','*wi_h0*','wi_b3gsm'],/delete store_data,['*kyoto_dst*'],/delete store_data,['t96_par','par_out'],/delete store_data,[rbspx+'_fgl_gsm'],/delete options,magname + '_' + model,'ytitle',rbspx + '!C'+strupcase(model)+' mag!CGSM!C[nT]' options,magname + '_' + model+'_wind','ytitle',rbspx + '!C'+strupcase(model)+' mag!CGSM!Cwith Wind input!C[nT]' options,magname + '_' + model+'_ace','ytitle',rbspx + '!C'+strupcase(model)+' mag!CGSM!Cwith ACE input!C[nT]' options,magname + '_' + model+'_omni','ytitle',rbspx + '!C'+strupcase(model)+' mag!CGSM!Cwith OMNI input!C[nT]' options,magname + '_' + model + '_dif','ytitle',rbspx + '!C'+strupcase(model)+' mag-model!CGSM!C[nT]' options,magname + '_' + model+'_wind_dif','ytitle',rbspx + '!C'+strupcase(model)+' mag-model!CGSM!Cwith Wind input!C[nT]' options,magname + '_' + model+'_ace_dif','ytitle',rbspx + '!C'+strupcase(model)+' mag-model!CGSM!Cwith ACE input!C[nT]' options,magname + '_' + model+'_omni_dif','ytitle',rbspx + '!C'+strupcase(model)+' mag-model!CGSM!Cwith OMNI input!C[nT]' options,magname + '_' + model,'ysubtitle','' options,magname + '_' + model+'_wind','ysubtitle','' options,magname + '_' + model+'_ace','ysubtitle','' options,magname + '_' + model+'_omni','ysubtitle','' options,magname + '_' + model + '_dif','ysubtitle','' options,magname + '_' + model+'_wind_dif','ysubtitle','' options,magname + '_' + model+'_ace_dif','ysubtitle','' options,magname + '_' + model+'_omni_dif','ysubtitle','' ;change name from T89 to IGRF if appropriate. ;Since the IGRF model is called with the T89 routine and a keyword it was easier to change ;IGRF -> T89 to get the code to work. Here I change it back. if model eq 'igrf' then begin get_data,rbspx+'_mag_bt89_original',data=dd store_data,rbspx+'_mag_bigrf_original',data={x:dd.x,y:dd.y} store_data,[rbspx+'_mag_bt89_original'],/delete endif end ;------------------------------------------------- ;TEST THE DIFFERENCE B/T THE IGRF AND T89 MODELS ;------------------------------------------------- ; call_procedure,'t'+model,rbspx+'_state_pos',kp=2.0 ; get_data,rbspx+'_state_pos_bt89',data=ddd ; store_data,'t89_tmp',data={x:ddd.x,y:ddd.y} ; call_procedure,'tt89',rbspx+'_state_pos',kp=2.0,igrf_only=1 ; get_data,rbspx+'_state_pos_bt89',data=ddd ; store_data,'igrf_tmp',data={x:ddd.x,y:ddd.y} ; get_data,'igrf_tmp',data=igrf ; get_data,'t89_tmp',data=t89 ; split_vec,'igrf_tmp' ; split_vec,'t89_tmp' ; dif_data,'igrf_tmp_x','t89_tmp_x',newname='diff_x' ; dif_data,'igrf_tmp_y','t89_tmp_y',newname='diff_y' ; dif_data,'igrf_tmp_z','t89_tmp_z',newname='diff_z' ; ; tplot,['igrf_tmp_x','t89_tmp_x','diff_x'] ; tplot,['igrf_tmp_y','t89_tmp_y','diff_y'] ; tplot,['igrf_tmp_z','t89_tmp_z','diff_z'] ;if model eq 'igrf' then call_procedure,'igrf',rbspx+'_state_pos',$ ; pdyn=2.0D,dsti=-30.0D,yimf=0.0D,zimf=-5.0D $ ;--------------------------- ;Rename the first model by adding "original" to name b/c it will otherwise be ;overwritten later when the Wind, ACE, OMNI models are called. ; get_data,rbspx+'_state_pos',dlimits=dlm ; get_data,rbspx+'_state_pos_b'+model,data=dat ; store_data,rbspx+'_state_pos_b'+model+'_original',data=dat,dlimits=dlm