Simulating Images

From RHESSI Wiki

(Difference between revisions)
Jump to: navigation, search
m (Flare loops)
Line 137: Line 137:
  cgg2.curvature = 0.0
  cgg2.curvature = 0.0
  model = [cgg1,cgg2]
  model = [cgg1,cgg2]
-
  im = obj_im -> gatdata(/plot,sim_model=model)
+
  im = img_obj -> getdata(sim_model = model)
-
 
+
  img_obj -> plot
-
Note that the amplitude for cgg2 is much less than for cgg1. This is necessary so that the point source (cgg2) does not overwhelm the loop (cgg1) in the image.  This possibility results from the fact that the amplitudes specify the total flux from the given component, so for large sources which are spread out over a substantial area, a larger relative amplitude is required to give them an appreciable maximum in the image.  To get a good image of the above source, you will likely have to improve the spatial sampling of the simulation and reconstructed image if you have not already done so.
+
 +
Note that the amplitude for cgg2 is much less than for cgg1.  This is necessary so that the point source (cgg2) does not overwhelm the loop (cgg1) in the image.  This possibility results from the fact that the amplitudes specify the total flux from the given component, so for large sources which are spread out over a substantial area, a larger relative amplitude is required to give them an appreciable maximum in the image.  To get a good image of the above source, you will likely have to improve the spatial sampling of the simulation and reconstructed image if you have not already done so.
== new hsi_switch ==
== new hsi_switch ==

Revision as of 15:20, 12 July 2012

This document gives some examples of how to simulate RHESSI images and eventlist files. The control parameters used are documented in the simulation control document . Read about eventlist files in Using Eventlist Files Note that it is not ncessary to create an eventlist file in order to do a simulation; you are still able to simulate images directly using the image object, and all of the simulation control parameters documented in the simulation control document can still be used directly from the hsi_image object. You do need to use the eventlist object if you want to save the data. The old method of using simulated telemetry packets, via the routines HSI_SIM_FLARE and HSI_SIM_FILE, is obsolete. The pre-launch simulated data files are obsolete; the packet reader does not read these files correctly any more.

Contents

Simple Simulation

Here is a simple example. First, call hsi_switch and set /sim,/aspect_sim

hsi_switch, /sim, /aspect_sim

Then, create the object, if you want to save the data in a file, create an hsi_eventlist object. If you don't want to save the data, then create an image or spectrum, or lightcurve object, depending on what you want. (Note that if you want to use a true aspect solution, and not a simulated one, you must save the eventlist in a file). You can set control parameters here or by using separate calls to that objects set method. Here we will simulate 4 seconds of data, the default, for the default source, which is a double gaussian source located at position [600.0, 200.0] but change the reference time.

o = hsi_eventlist(sim_ut_ref='2-jul-2003 02:20:00')

The getdata method returns the simulated eventlist.

ev = o -> getdata(sim_xyoffset = [600.0, 200.0])

To save the data, use the write method, with time_range=[0,0], energy_band=[0,0] to insure that all of the data is written to the file:

o -> Write, 'hsi_test.fits', time_range=[0, 0], energy_band=[0, 0]

To recover the image, or spectrum, create the object, set the filename with the ev_filename control parameter, and use getdata():

img_obj = hsi_image(ev_filename = 'hsi_test.fits')
image = img_obj -> getdata(xyoffset = [600.0, 200.0])

Note that we set the xyoffset in the image object. This is always necessary when using simulated data, otherwise the software will set the xyoffset to the nearest flare value.

Note also that you never want to set any sim_ control parameter for this image object. Setting a sim_ control parameter will cause the object to start a new simulation, even though you passed in a filename. This may change, but for now you need to be careful.

To see your image, try the plot method.

img_obj -> plot

All of the hsi_image methods and control parameters can be used as usual. For example, if you want to use the CLEAN algorithm, just set the image_alg control parameter.

img_obj -> getdata(image_alg = 'CLEAN')

Using the true aspect solution

In the above example, we used a simulated aspect solution for the simulation and for the image. For many applications, you might want to use the actual aspect solution obtained directly from the data. Here's how to do that: First set hsi_switch, /sim, or hsi_switch, /sim, /aspect_flight, then create the object. This ensures that the hsi_eventlist object knows that it wants to simulate the data, and not look for it in the level 0 data, but the aspect solution will know that it needs to obtain the solution from the level 0 data.

hsi_switch, /sim, /aspect_flight
o = hsi_eventlist(sim_ut_ref = '2-jul-2003 02:20:00')

Next call getdata, and write the eventlist file

ev = o -> getdata()
o -> Write, 'hsi_test_1.fits', time_range = [0, 0], energy_band = [0, 0]

Getting the image back out is the same as before, but as of 15-oct-2004, you need a call to hsi_switch in order for the aspect solution to be found correctly. You also will get a bunch of error messages of the form:

 HSI_EVENTLIST::SET: You are setting simulation parameters even though you
                have not selected Simulations (via $HSI_USE_SIM).  Not
                setting parameters.

Ignore those messages. This should be fixed in the future.

Now simply make the image.

img_obj = hsi_image(ev_filename = 'hsi_test_1.fits')
image = img_obj -> getdata()

Models

A model is passed through the sim_model parameter. For example, if you are simulating eventlists then

obj = hsi_eventlist(sim_ut_ref='2-jul-2003 02:20:00', sim_model = model)

where the variable model is the model represented as an array or a structure.

Point Source

If you want a 1 pixel size source to be located at sim_xyoffset, and pass in a (64,64) image, then use the following:

model = fltarr(64,64) 
model[31,31] = 1.0

If you pass in a (65,65) image then use:

model = fltarr(65,65)
model[32,32] = 1.0

It is also possible to define a point source using the gaussian_source_str (see next section) definition which defines a double gaussian for input into a sim_model, by setting the widths of the two gaussians to close to zero and placing them on top of each other.

model = {gaussian_source_str}              ;Input to hsi_modulate_point_source
model.amplitude = 1.0
model.xysigma = [0.0001, 0.0001]
model.xypos = [0.0, 0.0]

Gaussian Source

Most models of RHESSI data use Gaussian sources of various sizes and aspect ratios. Curved Gaussian sources can also be used to simulate flare loops. The suite of RHESSI software has a convenient means built in to simulate Gaussian sources. One must specify the source amplitude (normalized to unit maximum - the overall flux normalazation is handled in the number of photons specified), location with respect to the center of the simulation box (note: this is a relative position with respect to the center of the simulation box specified by the sim_xyoffset keyword), width in the X and Y directions, and the tilt angle. In the case of curved elliptical Gaussians, a curvature parameter must also be specified. This information is entered into a source structure for use by the software. For example, suppose you wish to specify a single Gaussian source offset +10 arcseconds in the X-direction and +2 arcseconds in the Y-direction from the center of the image, with a X width (sigma) of 5.0 arcseconds and a Y width of 2.0 arcseconds, and a tilt angle of 45.0 degrees, the following commands should be entered:

gg = {gaussian_source_str}
gg.amplitude = 1.0
gg.xysigma = [5.,2.]
gg.xypos = [10.,2.]
gg.tilt_angle_deg = 45.0

You can specify multiple Gaussian sources as well. Suppose you wish to specify 2 sources, both with widths of [5.0,2.0], with one located at an offset [20.0,-3.0] and one half as strong as this at [-5.0,15.0], both with 0.0 tilt angles, you would type:

gg1 = {gaussian_source_str}
gg1.amplitude = 1.0
gg1.xysigma = [5.,2.]
gg1.xypos = [20.,-3.]
gg1.tilt_angle_deg = 0.
gg2 = {gaussian_source_str}
gg2.amplitude = 0.5
gg2.xysigma = [5.,2.]
gg2.xypos = [20.,-3.]
gg2.tilt_angle_deg = 0.
model = [gg1,gg2]

Flare loops

The RHESSI simulation software can also be used to specify a curved, elliptical Gaussian which is a good proxy for a solar flare loop. Instead of declaring your source as a {gaussian_source_str}, it is declared as a {curved_gaussian_source_str}. In this case there is one additional structure tag (.curvature) which must be specified. For example, suppose you want to declare a curved Gaussian with the same basic parameters as given in the first Gaussian source example above. You would simply type:

cgg = {curved_gaussian_source_str}
cgg.amplitude = 1.0
cgg.xysigma = [10.,2.]
cgg.xypos = [10.,2.]
cgg.tilt_angle_deg = 45.
cgg.curvature = 0.2

This last parameter specifies the curvature by giving the inverse of the radius of curvature of the source, normalized by the larger of the two Gaussian widths. Basically, the radius of curvature of the source is larger of the two Gaussian widths divided by cgg.curvature. In the example above then, the radius of curvature is 25 pixels, which should be quite noticeable in RHESSI data.

f you wish to specify a model with a curved, elliptical Gaussian plus some number of point like sources, all components should be specified using the {curved_gaussian_source_str} system, with the point sources given a curvature of 0. For example, if you wish to make a model with the curved source specified above plus a weaker, symmetric point source (not curved) above the loop, you can type:

cgg1 = {curved_gaussian_source_str}
cgg1.amplitude = 1.0
cgg1.xysigma = [10.,2.]
cgg1.xypos = [10.,2.]
cgg1.tilt_angle_deg = 45.
cgg1.curvature = 0.2
cgg2 = {curved_gaussian_source_str}
cgg2.amplitude = 0.2
cgg2.xysigma = [2.,2.]
cgg2.xypos = [30.,25.]
cgg2.tilt_angle_deg = 0.
cgg2.curvature = 0.0
model = [cgg1,cgg2]
im = img_obj -> getdata(sim_model = model)
img_obj -> plot

Note that the amplitude for cgg2 is much less than for cgg1. This is necessary so that the point source (cgg2) does not overwhelm the loop (cgg1) in the image. This possibility results from the fact that the amplitudes specify the total flux from the given component, so for large sources which are spread out over a substantial area, a larger relative amplitude is required to give them an appreciable maximum in the image. To get a good image of the above source, you will likely have to improve the spatial sampling of the simulation and reconstructed image if you have not already done so.

new hsi_switch

The major new development in the sim software is the ability to use the true aspect solution, instead of the original simulated aspect solution. This is controlled using the routine hsi_switch; if the call is:

hsi_switch, /sim, /aspect_sim

Then the eventlist and the aspect solution are simulated.

To use a simulated eventlist and a true aspect solution, call

hsi_switch, /sim, /aspect_flight

Or

hsi_switch, /sim

To reset to the default behavior, which is eventlist from packets, and aspect solution from packets, call:

hsi_switch, /flight, /aspect_flight

Or

hsi_switch

Note that a call to hsi_switch always resets the environment, and the default for both eventlist and aspect data is to use the flight data, i.e., from the packets. This means that if you make two calls:

hsi_switch, /sim
hsi_switch, /aspect_sim

You'll end up with a simulated aspect solution, but a non-simulated eventlist, and the programs will crash. So please be sure to pass keywords for both the aspect solution and eventlist in the hsi_switch calls.


Originally written 16-Oct-2004, jmm Imported --Schriste 21:52, 5 November 2008 (UTC)

References

Personal tools
Namespaces
Variants
Actions
Navigation
Toolbox