RHESSI Command Line Snippets

 

Listed below are snippets of IDL code to do various things with the RHESSI objects.  These snippets show examples of the kinds of things that can be done.  This is by no means a complete list of options  - just samples showing some of the syntax and methods that might help you use the RHESSI data from the command line. 

In most cases, the routines/methods shown below have additional calling arguments that you can find by looking at other more complete documentation, or in the header of the routine.  The full table of RHESSI object parameters is here.

Contents:
Flare List
Observing Summary Object
Image Object
Spectrum Object
Eventlist
Aspect Solution
Binned Eventlist
Calibrated Eventlist
Point Spread Function
Monitor Rates
Simulations

GOES Object
 


FLARE LIST

Read flare list into structure and find flares with energies > 7000:

a = hsi_read_flarelist()
index = where (a.energy_hi[1] gt 7000.)
print,a[index].id_number

or

print, hsi_select_flare(energy_hi_range=[7001,1.e7])

flare_struct = hsi_select_flare(energy_hi_range=[7001,1.e7], /struct)    ; returns structure

Find flares that occurred between select times that have a peak countrate > 500:

flare_ids = hsi_whichflare(['2002/9/2 00:00', '2002/9/2 12:00'])
flare_str = hsi_getflare(flare_ids)
index = where (flare_str.peak_countrate gt 500)
print,flare_str[index].id_number

or

print, hsi_select_flare(peak_time_range = ['2002/9/2 00:00', '2002/9/2 12:00'], $
    peak_countrate_range=[501,1.e9])

Get array of structures for flares that were in attenuator state 3 at the peak of the flare:

flare_struct = hsi_select_flare(flag_incl='A3', /struct)

Print list of flare flags (list contains name of flag in FITS file, 2 letter code, and description):

print,hsi_flare_flag_code(/expand,/sort)
 


OBSERVING SUMMARY OBJECT

Plot observing summary rates:

obs_object = hsi_obs_summary()
obs_object -> set, obs_time_interval= $
    [' 20-feb-2002 10:56:00.000', ' 20-feb-2002 11:21:00.000']

obs_object -> plot    ; plots rates
obs_object -> plotman, /corrected    ; plots corrected rates in plotman

Plot modulation variance with night, flare flags:

obs_object -> plotman, class='mod_var', /eclipse, /flare

Notes: Set class_name keyword to 'mod_var', 'roll_period', 'roll_angle', or 'ephem', 'pointing' to plot modulation variance, roll_period, roll_angle, ephemeris, or pointing data.  You can use most of the standard plot keywords (yrange, ylog, ystyle, etc), and either plot or plotman for all of these plots, e.g.,

obs_object -> plotman, yrange=[1,200], /ylog, ystyle=1
obs_object -> plot, class='roll_period'
obs_object -> plotman, class='ephem'

Plot corrected observing summary rates and roll angle manually:

time = obs_object -> getaxis(/ut)
data = obs_object -> getdata(/corrected)
utplot, anytim(time, /ext), data.countrate[0]    ; plots lowest energy band vs time

time = obs_object -> getaxis(/ut, class='roll_angle')
data = obs_object -> getdata(class='roll_angle')
utplot, anytim(time, /ext), data.roll_angle

Print night time intervals on 20-feb-2002:

obs_object -> set, obs_time_interval=[' 20-feb-2002 00:00', ' 21-feb-2002 00:00']
changes = obs_object -> changes()
dn = changes.eclipse_flag
q = where (dn.state eq 1, nq)
for i=0,nq-1 do ptim, [dn.start_times[q[i]], dn.end_times[q[i]]]

Note: type help, dn, /struct to see names of all flags in list

Get attenuator states for a time interval:

t='21-apr-2002 ' + ['00:00','02:00']
obs=hsi_obs_summary(obs_time=t)
changes = obs->changes()
atten = changes.attenuator_state
help,atten,/st


Plot flare and attenuator state change bars (see hsi_show_flags__define for other flags to display) near top of your existing utplot plot (must be last plot drawn) :

t='21-apr-2002 ' + ['00:00','02:00']
show_flags_obj = hsi_show_flags()
show_flags_obj  -> set,/flare,/atten
show_flags_obj  -> display,t

 


IMAGE OBJECT

Single image:

image_object = hsi_image()
image_object -> set, im_time_interval = ['20-feb-2002 11:06:18', '20-feb-2002 11:06:22']
image_object -> set, im_energy_binning = [3,12]
data = image_object -> getdata()

image_object -> plotman, colortable=5
loadct,3 & image_object -> plot, xrange=[850,950], yrange=[200,300], /limb,/cbar, legend_loc=0, title='', xtitle='', ytitle='', /no_timestamp

image_object -> fitswrite, this_out_filename='hsi_image_test.fits'

Image cube with 7 time bins, 6 energy bins:

image_object -> set, im_time_interval = ['20-feb-2002 11:06:00', '20-feb-2002 11:07:00']
image_object -> set, im_time_bin = 4.
image_object -> set, im_energy_binning=[3,6,12,25,50,100,300]
image_object -> fitswrite
image_object -> plot, t_idx=3, e_idx=2        ; plot 3rd time, 2nd energy
image_object -> plotman, /choose

Read image FITS file into object:

image_object -> set, im_input_fits = 'hsi_image_test.fits'
data = image_object -> getdata()

Read image FITS file without using object:

images = hsi_image_fitsread(fitsfile='hsi_image_test.fits')

Get image x, y, time, energy axes midpoints (use /edges_2, /edges_1, /width for different formats):

xaxis = image_object -> getaxis(xaxis)
yaxis = image_object -> getaxis(/yaxis)
taxis = image_object -> getaxis(/ut)
eaxis = image_object -> getaxis(/energy)


SPECTRUM OBJECT

Set time by specifying width, and energy using code: 

spectrum_object = hsi_spectrum()
spectrum_object-> set, obs_time_interval= [' 20-feb-2002 10:56:00.000', ' 20-feb-2002 11:21:00.000']
spectrum_object-> set, sp_time_interval= .2    ; time bins are .2 seconds wide
spectrum_object-> set, seg_index_mask= [1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]    ; use front segments 1,3,4,6,8,9,
spectrum_object -> set, sp_energy_binning=1    ; energy code 1
data = spectrum_object -> getdata()

List what the RHESSI binning codes mean:

ebfile = hsi_loc_file('energy_binning.txt', path='$HSI_SPEC', count=count)
more, rd_ascii (ebfile[0])

Set time and energy binning explicitly:

spectrum_object -> set, sp_time_interval = anytim(' 20-feb-2002 11:05:00.000') + findgen(100)*4.
spectrum_object-> set, sp_energy_binning = findgen( 398 ) + 3    ; 1-keV bins from 3 to 400 keV

Get midpoints of time and energy bins (use /edges_2, /edges_1, /width for different formats):

time_bins = spectrum_object -> getaxis (/ut)
energy_bins = spectrum_object -> getaxis (/energy)

Plot using spectrum object plot method:

spectrum_object -> plotman    ; plot spectrum
spectrum_object -> plot, /pl_time, /dim1_sum    ; plot time profile, summing over energy bins
spectrum_object -> plotman, /pl_spec    ;plot spectrogram

Plot flux spectrum and livetime manually:

data = spectrum_object -> getdata(sp_data_unit='flux', /sp_data_structure)     ; other options are 'counts', 'rate'
emid = spectrum_object -> getaxis(/energy)
plot, emid, data[0].flux, /xlog, /ylog, psym=10, yrange=[.001,1.]    ; first time interval only

utplot, anytim ( average(data.ut,1), /ext), data.ltime

Set to native energies:

a2d = 4 ; index for detector 5
e1 = 1. & e2=100.   ; min, max native energies to use
spectrum_object->native_energy, a2d, this_energy_range=[e1,e2]

Compute native energies at a given time for a detector:

en = hsi_get_e_edges(gain_time_wanted='2-apr-2011 23:00', a2d=1)  ; for detector 2

Write spectrum and SRM FITS file:

spectrum_object-> filewrite, /buildsrm, all_simplify=0  ;full SRM, default names

Note: use /namedialog to enable widget for picking output file names and locations.  Use srmfile='xxx.fits', specfile='yyy.fits' to name files explicitly.

Write separate-detector spectrum and SRM FITS files (NON-NATIVE energy bins):

hsi_spectrum_sep_det_files, flare = [3111313, 5082502], /not_native, spec_dir = 'c:\rhessi\', det=[4, 8]  ; uses default energy binning

hsi_spectrum_sep_det_files, flare=[11040103, 11040212], det=[1, 3, 4],  /not_native, sp_energy_binning=29  ; uses energy binning code 29

hsi_spectrum_sep_det_files, obj=obj, det=2, /not_native ; uses what's set in object, writes files for Det 2

hsi_spectrum_sep_det_files, times='7-apr-2012 ' + ['09:24', '09:47'], /not_native, sp_energy_binning=10  ; write files for 9 dets, specified times, energy binning code 10

hsi_spectrum_sep_det_files, times='7-apr-2012 ' + ['09:24', '09:47'], /not_native, sp_energy_binning=findgen(500)+5.)  ; specify energy bins explicitly

Write separate-detector spectrum and SRM FITS files with NATIVE energy bins:

hsi_spectrum_sep_det_files, flare = [3111313, 5082502], max_energy=100, spec_dir = 'c:\rhessi\', /plots, det=[4,8]

generates spectrum and srm files for RHESSI detectors 4 and 8 for events with flare number 3111313 and 5082502 in the RHESSI flare catalog, sends plots to screen of count flux vs. time, uses max_energy 0f 100 kev, saves files in directory c:\rhessi (will create directory if it does not already exist)

hsi_spectrum_sep_det_files, times = ['13-Nov-2003 04:00:00','13-Nov-2003 06:00:00'], max_energy = 400, /extend_times

generates spectrum and srm files for each RHESSI detector for selected time interval extended to include surrounding night, uses max_energy of 400 kev, saves files in current directory

hsi_spectrum_sep_det_files, obj=obj, max_energy=300

generates spectrum and srm file for time interval and detectors set in spectrum object called obj, with max_energy of 300


EVENTLIST

Read and plot eventlist data:

ev_object = hsi_eventlist(time_range=[0.,0.], obs_time_interval=
dall = ev_object ->getdata()
q = where(dall.a2d_index eq 5) ; get Detector 6 data
d5 = dall[q]
ut_ref = ev_object->get(/ut_ref)
tarr = d5.time / 2d^20+ ut_ref
utplot, tarr-tarr[0], d5.channel, tarr[0], psym=3

or plot in plotman:

utobj = obj_new('utplot', tarr-tarr[0], d5.channel, utbase=tarr[0])
utobj->plotman, psym=3

Create an eventlist file:

ev_object = hsi_eventlist(time_range=[0.,0.], obs_time_interval=
ev_object -> set, energy_band=[1800, 2200]
ev_object -> write

Use an eventlist file as input:

o = hsi_image(ev_filename='hsi_eventlist...
o->set, im_energy_binning=[1800,2200]
image = o->getdata()
 


ASPECT SOLUTION

Get aspect solution used in image object:

o = image_object -> getdata(class='hsi_aspect_solution')

Create standalone aspect solution for specified time interval:

o = hsi_aspect_solution(obs_time_interval=['20-feb-2002 11:06:00', '20-feb-2002 11:07:00'])

Plot pointing axis in xy plot, or get aspect data and plot pointing axis vs time yourself:

o -> plot

data = o -> getdata()
x = data.pointing[*,0]
y = data.pointing[*,1]
t = hsi_sctime2any(data.t0) + data.time*2.^(-7)
utplot, anytim(t,/ext), sqrt(x^2 + y^2)

Get spin axis position for each rotation:

spin_axis = o -> get_spin_axis()
spin_axis_ave = o -> get_spin_axis(spin_period=spin_period, /average) ; also retrieve period


BINNED EVENTLIST

Binned eventlist returns [9,3] array of pointers. For a valid pointer (in this case [0,0]), get structure and show gaps and livetime:

be = image_object -> getdata(class='hsi_binned_eventlist)
print,  ptr_valid(be)    ;  show which pointers are valid
be00 = *be[0,0]
plot,be00.time, b.gap, yrange=[0,2]
plot,be00.time, b.livetime


CALIBRATED EVENTLIST

Calibrated eventlist returns [9,3] array of pointers.  For valid pointers,

cbe = image_object -> getdata(class='hsi_calib_eventlist')
print, ptr_valid(cbe)


POINT SPREAD FUNCTION

Get point spread function  used in CLEAN image algorithm and plot it:

image_object -> set, image_alg='clean'
psf = image_object -> getdata(class='hsi_psf')
loadct,5 & tvscl, congrid(psf,512,512)
 


MONITOR RATES

monitor_object = hsi_monitor_rate()
monitor_object ->set, obs_time_interval=$
    ['20-feb-2002 11:06:00', '20-feb-2002 11:07:00']
d=monitor_object -> getdata()


SIMULATIONS

hsi_switch, /sim, /aspect_flight

s=hsi_eventlist(sim_ut_ref='2-sep-2002 00:07:36', $
    sim_time_range=[0,44.], sim_energy_band=[23,52], sim_photons_per_coll=1e5)
s->set, sim_model=ima, sim_xyoffset=[939,0.]
s->write,'simtest.fits'
 


GOES Object

Start GOES GUI with command line access to the GOES object:

goes,g    ;start GOES GUI, and return object in variable g

or

g = ogoes()
g -> gui

Create text file with information other than what's in the text file created by the GOES object.  in this case the output file will have columns for time at center of bin, energy loss rate, the integrated energy loss, and the derivative of the flux in the two GOES channels:

dintegr = g->getdata(/struct,/integrate)
lrad_integr = dintegr.lrad

d = g->getdata(/struct)

times = anytim(anytim(d.utbase)+d.tarray, /vms)
deriv0 = deriv(d.tarray, d.yclean[*,0])
deriv1 = deriv(d.tarray, d.yclean[*,1])
str = string (transpose([[d.lrad], [lrad_integr], [deriv0], [deriv1]]), format='(4g15.5)')
out =times + ' ' + str
prstr,file='lrad.txt',out
 


Plotman

After plotting an image in plotman, to change the ticks to point outward:

If pobj is your plotman object (or if in hessi GUI, get the reference via hessi_data, plotman=pobj),
panel = pobj->get(/current_panel_struct)
data_obj = *panel.saved_data.data
data_obj->plotman, plotman_obj=pobj, ticklen=-.02, nodup=0

or if you want to give the panel a new name:
data_obj->plotman, plotman_obj=pobj, ticklen=-.02, panel_desc='ticks'

The outward ticks will persist with that panel, with the overlay of another image on this image, and the writing of a plot file. Note: this assumes the panel you're modifying to have outward ticks was made from an object that has a plotman method (e.g. map, xyplot, utplot, etc.)


Kim Tolbert
Last Modified: 12 December 2018