;+
;Procedure: ttrace_crib
;
;Purpose: crib demonstrating the usage of tracing routines
;
;Notes: tracing may take a while to complete, don't worry it is
;not hung.
;Haje Korth's IDL/Geopack DLM must be installed for this
;       to work
;
; $LastChangedBy: jimm $
; $LastChangedDate: 2008-05-13 18:17:46 -0700 (Tue, 13 May 2008) $
; $LastChangedRevision: 3082 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/idl_socware/tags/tdas_5_11/external/IDL_GEOPACK/trace/ttrace_crib.pro $
;-

;this makes all int declarations long by default and doesn't

;allow ambiguous array access (ie a(i) is illegal only a[i] is

;acceptable)

compile_opt idl2

;make sure a color table is loaded that will display all the lines

thm_init
;BASIC USAGE

;footprint trace in gsm

timespan,'2007-03-23'

thm_load_state,probe='c',coord='gsm',suffix='_gsm'

;input to the t89 model(par value) is the kp index(can be an array of

;kp indices

ttrace2iono,'thc_state_pos_gsm',newname='thc_ifoot_gsm',external_model=$

't89',par=2.0D,/km

ttrace2equator,'thc_state_pos_gsm',newname='thc_efoot_gsm',$

external_model='t89',par=2.0D,/km

;increase the x margin so that you can see the labels on the left

tplot_options, 'xmargin', [15, 5]

tplot,['thc_state_pos_gsm','thc_ifoot_gsm','thc_efoot_gsm'], title = $

'Themis C position & footprints'

stop

;MAP PROJECTIONS OF TRACES

;now load in geo

;for projection of ionospheric footprint on map

del_data,'*'

thm_load_state,probe='c',coord='geo'

ttrace2iono,'thc_state_pos',newname='thc_ifoot',external_model='t89',$

par=2.0D,/km,in_coord='geo',out_coord='geo'

xyz_to_polar,'thc_ifoot'

get_data,'thc_ifoot_phi',data=d

i_lon=d.y

get_data,'thc_ifoot_th',data=d

i_lat=d.y

;generate a map,(centered on some random point in the projection)

;if you want to see the projection with ground station data, see 

;themis/examples/thm_crib_asi.pro

;look at the plain map_set routine to see other types of maps that can

;be generated for plotting

thm_map_set,central_lon=i_lon[400],central_lat=i_lat[400],xsize=700,ysize=400

;ionospheric footpoint trace plotted on the map

plots,i_lon,i_lat

stop

del_data,'*'

;reset the colors because the map routine switches the color table

device,decomposed=0

loadct2,43

;2D TRACE PLOTS

;load data

thm_load_state,probe='c',coord='gsm'

tKm2Re,'thc_state_pos',/replace

;generate the trace

ttrace2equator,'thc_state_pos',newname='thc_efoot',external_model='t89'$

,par=2.0D,in_coord='gsm',out_coord='gsm'

;plot y vs x with y having its maximum on the left

;NOTE: if the plot goes off the edge of your screen you should
;decrease the xsize and ysize keywords...the default in this crib
;is set for high resolution
tplotxy,'thc_state_pos',multi='2 2',versus='xryr', xtitle = 'Re', ytitle = 'Re',title='x vs y, pos and e-foot',window=0,xsize=1000,ysize=1000,wtitle='position & footprint plots'

tplotxy,'thc_efoot',versus='xryr',/overplot,linestyle=2,/isotropic

;plot y vs z with y having its maximum on the left

tplotxy,'thc_state_pos',/add,versus='yrz', xtitle = 'Re', ytitle = 'Re',/$

isotropic,title='y vs z, pos and e-foot'

tplotxy,'thc_efoot',versus='yrz',/overplot,linestyle=2,/isotropic

;plot x vs z with x having its maximum on the left

tplotxy,'thc_state_pos',/add,versus='xrz', xtitle = 'Re', ytitle = 'Re',/$

isotropic,title='x vs z, pos and e-foot'

tplotxy,'thc_efoot',versus='xrz',/overplot,linestyle=2,/isotropic

stop

;2D PLOT FIELD LINES WITH TRACES

;generate some dummy data to generate field line traces

t = replicate(time_double('2007-03-23/11:00:00'), 18)

x = -1 * (dindgen(18)+2)^2

y = replicate(0D, 18)

z = replicate(0D, 18)

v = [[x], [y], [z]]

str = {x:t, y:v}

;this stores the dummy data

store_data, 'trace_data', data = str

;get the traces

ttrace2iono,'trace_data',trace_var_name = 'trace_n', newname=$

'thc_ifoot',external_model='t89',par=2.0D,in_coord='gsm',out_coord=$

'gsm'

ttrace2iono,'trace_data',trace_var_name = 'trace_s', newname=$

'thc_ifoot',external_model='t89',par=2.0D,in_coord='gsm',out_coord=$

'gsm', /south

;generate the 2-d plot

timespan,'2007-03-23/10:00:00',2,/hour

tplotxy,'thc_state_pos',/add, versus = 'xrz', xtitle = 'Re', ytitle = 'Re',$

title = 'fields x vs z',xrange=[-10,-4],yrange=[-6,4]

tplotxy,'thc_efoot', versus = 'xrz', /overplot,linestyle=5

tplotxy, 'trace_n', versus = 'xrz', /overplot, linestyle = 2

tplotxy, 'trace_s', versus = 'xrz', /overplot, linestyle = 2

stop

;HOW TO USE THE OTHER MODELS

;these models are quite a bit slower than t89 so don't worry they will

;finish (-:

timespan,'2008-02-02',1,/day

thm_load_state,probe='c',coord='gsm'


;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/YOURUSERNMAE/My Documents' ;example  if you don't want to use the default windows location (C:/data/ or E:/data/)

;load kyoto dst
kyoto_dst_load

;load wind data
wi_mfi_load
wi_3dp_load

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','T96'

;arguments to t96 are:

;solar wind pressure nP

;dst index nT

;y component of interplanetary magnetic field

;z component of interplanetary magnetic field

;par must be a 10 element array(or Nx10) so all other values should be

;0d

;longer timespans take a long time to run

timespan,'2008-02-02/12:00:00',30,/minute

time_clip,'thc_state_pos','2008-02-02/12:00:00','2008-02-02/12:30:00',$

newname='pos_clipped'

ttrace2iono,'pos_clipped',newname='thc_ifoot_gsm',external_model=$

't96',par='t96_par',/km

tplot,['pos_clipped','thc_ifoot_gsm'], title = $

'Themis C position & footprints'

stop

tplotxy,'pos_clipped',multi='2,1',title='thc_state_pos',versus='xrz',xtitle='X(km)',ytitle='Z(km)',xsize=1000,ysize=500,window=2

tplotxy,'thc_ifoot_gsm',/add,title='thc_ifoot_gsm',versus='xrz',xtitle='X(nT)',ytitle='Z(nT)'

stop

;arguments to t01 are:

;solar wind pressure nPa

;dst index nT

;y component of interplanetary magnetic field

;z component of interplanetary magnetic field

;g1,g2 indices describing solar wind conditions

;generate g1,g2 with geopack_getg

;par must be a 10 element array(or Nx10) so all other values should be

;0d

;par = [2.0D,-30.0D,0.0D,-5D,g[0,0],g[0,1],0D,0D,0D,0D]

;parameters will be generated by get_tsy_params

get_tsy_params,'kyoto_dst','wi_b3gsm',$

'wi_3dp_k0_ion_density','wi_3dp_k0_ion_vel','T01'

ttrace2iono,'pos_clipped',newname='thc_ifoot_gsm',external_model='t01',$

par='t01_par',/km

tplot,['thc_ifoot_gsm','pos_clipped']

stop

tplotxy,'pos_clipped',multi='2,1',title='thc_state_pos',versus='xrz',xtitle='X(km)',ytitle='Z(km)',xsize=1000,ysize=500,window=2

tplotxy,'thc_ifoot_gsm',/add,title='thc_ifoot_gsm',versus='xrz',xtitle='X(nT)',ytitle='Z(nT)'

stop

;geopack_getw generates arguments for the ts04 model

;arguments are:

;solar wind density cm^-3

;solar wind speed km/s

;z component of interplanetary magnetic field nT

;if you pass in N length arrays for each of the arguments it will

;produce an Nx2 length output

;All arrays of data must be sampled at 5 minute intervals

;these parameters represent integrals over the timeseries

;so history is important(thus we use a whole series of ace params)

get_tsy_params,'kyoto_dst','wi_b3gsm',$

'wi_3dp_k0_ion_density','wi_3dp_k0_ion_vel','T04s'

;arguments to ts04 are:

;solar wind pressure nPa

;dst index nT

;y component of interplanetary magnetic field

;z component of interplanetary magnetic field

;w1,w2,w3,w4,w5,w6 storm time integral indices

;par can be a 10 element or an Nx10 element array

ttrace2iono,'pos_clipped',newname='thc_ifoot_gsm',external_model='t04s'$

,par='t04s_par',/km

;same as above w/r/t time range

tplot,['thc_ifoot_gsm','pos_clipped']

stop

tplotxy,'pos_clipped',multi='2,1',title='thc_state_pos',versus='xrz',xtitle='X(km)',ytitle='Z(km)',xsize=1000,ysize=500,window=2

tplotxy,'thc_ifoot_gsm',/add,title='thc_ifoot_gsm',versus='xrz',xtitle='X(nT)',ytitle='Z(nT)'

stop

;here is an example using ace data
;it may take a little while to complete

ace_mfi_load
ace_swe_load

;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','T96',/speed

ttrace2iono,'thc_state_pos',newname='thc_ifoot_gsm',external_model=$

't96',par='t96_par',/km

tplot,['thc_state_pos','thc_ifoot_gsm'], title = $

'Themis C position & footprints'

stop

tplotxy,'pos_clipped',multi='2,1',title='thc_state_pos',versus='xrz',xtitle='X(km)',ytitle='Z(km)',xsize=1000,ysize=500,window=2

tplotxy,'thc_ifoot_gsm',/add,title='thc_ifoot_gsm',versus='xrz',xtitle='X(nT)',ytitle='Z(nT)'

stop

;here is an example using omni_hro_data
;it may take a little while to complete
;also omni data has nan flags in any data gaps
;so you may want to use tdeflag to interpolate
;across data gaps

omni_hro_load

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','T96',/speed,/imf_yz

ttrace2iono,'thc_state_pos',newname='thc_ifoot_gsm',external_model=$

't96',par='t96_par',/km

tplot,['thc_state_pos','thc_ifoot_gsm'], title = $

'Themis C position & footprints'

end