Chapter 4. Mini Tutorial

Table of Contents

1. Using the Command Line
1.1. raytracewl
1.2. Using WCS information from a SECCHI data header: rtraytracewcs
1.3. Electron density cube generation: buildcloud
2. Fit of a flux rope on SECCHI data using the IDL-GUI
2.1. Demo mode
2.2. Preparing the data
2.3. Manual fitting of the flux rope
2.4. Generating Thomson scattering views
2.5. Saving the parameters and images of the fit
2.6. Automatic fit of the model position.
2.7. Sensitivity analysis.
3. Software concept overview
3.1. Units
3.2. Coordinate Systems
3.3. The Sun
3.4. Observer
3.5. Model Representation
3.6. Physical process

1. Using the Command Line

The raytracing can be performed via the command line of IDL. It is of course less friendly than a graphic interface front end but it provides a way more flexibility, and simulations can be included in a batch script for example. The two main routines are raytracewl and rtraytracewcs. The following few examples show how to use these routines. Please refer to the program headers, they should be up to date with all the parameters and keyword description.

1.1. raytracewl

1.1.1. No parameter passing

Example 4.1. Slab model simulation.

; -- simulation of a slab, edge-on
raytracewl,sbt1,sbp1,sne1,imsize=[256,256],losrange=[-30,30],$
losnbp=128L,modelid=14,neang=[0.,-39.5,-90-12]*!dtor,/c3
wnd,0,alog10(sbt1.im > 1e-14)

; -- simulation of a slab, face-on
raytracewl,sbt2,sbp2,sne2,imsize=[256,256],losrange=[-30,30],$
losnbp=128L,modelid=14,neang=[90.,-39.5,-90-12]*!dtor,/c3
wnd,1,alog10(sbt2.im > 1e-14)


Inputs

sbt1

Total brightness structures containing the image, the simulation parameters and a FITS header.

sbp1

Polarized brightness structures containing the image, the simulation parameters and a FITS header.

sne1

Integrated Ne structures containing the image, the simulation parameters and a FITS header.

Outputs

imsize

Size of the image in pixels.

losrange

LOS integration range in R_Sun, with the origin at the impact distance. Use the frontinteg keyword to change the origin to the observer.

losnbp

Number of integration points along the LOS, within the LOS range.

modelid

ID of the model.

neang

Rotation angles applied to the density model, in radian. Rotation order is Z, Y then X. X points up in the image, Y points on the right, and Z is the optical axis. (X,Y,Z) is direct.

/c3

Set to the predefined LASCO C3 FOV.

1.1.2. Rotation order convention

Rotation order convention for the variables obsang and neang is: 1st: Oz, 2nd: Oy, 3rd: Ox.

1.1.3. With parameter passing

Example 4.2. CME model simulation

raytracewl,sbt1,sbp1,sne1,imsize=[256,256],losrange=[-30,30],$
losnbp=128L,modelid=54,modparam=[2.,30.*!dtor,12.,0.2,1e6,0,0,0,0.8,0.5],neang=[50.,-39.5,-90-12]*!dtor,/c3
wnd,0,bytscl(alog10(sbt1.im > 1e-20),-13,-9),/tv


Inputs
modparam

Array of parameter corresponding to the model 54. See models51to60.cpp here for more details.

1.1.4. With a Density Cube

The model 26 is useful to do a ray-tracing through an electron density cube. The format of the modparam parameter array is the following (see also the source code):

  1. x size (sx)

  2. y size (sy)

  3. z size (sz)

  4. xc Sun center in pix

  5. yc Sun center in pix

  6. zc Sun center in pix

    [Note]Note

    (0,0,0) is at the first vertice of the density cude.

  7. voxel size in rsun, same for the 3 directions of space

  8. data cube in lexicographical order (x,y,z)

We can build a fake electron density cube and make the raytracing for the example. Here we use a parallelopiped slab. Note that model 26 uses trilinear interpolation between neighbor voxels. For no interpolation, use model 25.

; ---- build the fake density cube
cube=fltarr(64,64,64)
cube[32:*,32:*,32:38]=1e4
; ---- build the parameter array
modparam=[64,64,64,32,32,32,0.8,reform(cube,64L*64L*64L)]
; ---- generate the view
raytracewl,sbt,imsize=[256,256],losnbp=64L,losrange=[-20,20],$$
modelid=26,modparam=modparam,neang=[90,-80,30]*!dtor,/cor2
wnd,0,alog10(sbt.im >1e-12 )

1.1.5. With Source Surface Field Map

For model #11, a Carrington map with the position of the neutral line have to be used. Electron density is placed depending on the distance to the neutral line and radial distance to the Sun. This can be used to reproduce the shape of the streamer belt.

; ---- download the CR2012 SSFM from WSO
rdtxtmagmap,nsheetmap,crot=2012
; ---- format parameter array in single row vector
modparam=reform(nsheetmap,n_elements(nsheetmap))
; ---- run raytracing
raytracewl,sbt,imsize=[256,256],losnbp=64L,losrange=[-20,20],modelid=11,modparam=modparam,neang=[0,90,0.]*!dtor,/cor2
; ---- display total brightness
wnd,0,alog10(sbt.im > 1e-14)

1.1.6. Overploting Radial Distance

In this example, we show how to plot the radial distance on simulated images. The raytracing program returns an extra parameter rho which is the impact distance for each LOS/pixel of the image. oplotimpactgrid is then used to overplot that distance on the simulated image.

raytracewl,sbt1,sbp1,sne1,imsize=[256,256],losrange=[-30,30],$$
losnbp=128L,modelid=33,modparam=[0.7,2.55,30.,10,.2],neang=[0.,-39.5,-90-12]*!dtor,/c3,/fake,rho=rho
wnd,0,alog10(sbt1.im > 1e-14)
oplotimpactgrid,rho,levels=[1,(findgen(7)+2)*4],$$
c_labels=replicate(1,10),ystyle=5,xstyle=5,c_linestyle=1

1.2. Using WCS information from a SECCHI data header: rtraytracewcs

1.2.1. Simulation of a CME seen in SECCHI COR2-A

For this example, we take the SECCHI COR2-A event of 2007/05/15 23:52:30. A pre-event image taken on 2007/05/15 18:52:30 is subtracted and secchi_prep is used to calibrate both the pre-event and the event images from a triplet of polarized images. The IDL routine rtraytracewcs takes a SECCHI FITS header in input to set the simulation parameters: detector FOV, position and attitude. In the example code shown bellow, we use the GCS CME model, modelid #54. The 3 parameters of the variable neang allow to set respectively the Carington longitude, latitude and the tilt angle of the CME. The data image and the corresponding simulated image is given in Figure 4.1, “Simulation of a CME as seen from SECCHI COR2-A”

eventtriplet = path + ['secchi/lz/L0/a/seq/cor2/20070515/20070515_235230_s4c2A.fts',$
'secchi/lz/L0/a/seq/cor2/20070515/20070515_235243_s4c2A.fts',$
'secchi/lz/L0/a/seq/cor2/20070515/20070515_235256_s4c2A.fts']

preevtriplet = path + ['secchi/lz/L0/a/seq/cor2/20070515/20070515_215230_s4c2A.fts',$
'secchi/lz/L0/a/seq/cor2/20070515/20070515_215243_s4c2A.fts',$
'secchi/lz/L0/a/seq/cor2/20070515/20070515_215256_s4c2A.fts']

secchi_prep,eventtriplet,hdrevent,imevent,/polariz_on,/rotate_on,/PRECOMMCORRECT_ON,/rotinterp_on,/silent
secchi_prep,preevtriplet,hdrpreev,impreev,/polariz_on,/rotate_on,/PRECOMMCORRECT_ON,/rotinterp_on,/silent
; ---- display the data image
m=get_smask(hdrevent)
wnd,0,alog10(m*(imevent-impreev) > 1e-12 < 1e-10),.25

; -- get the default parameters for model 54
rtraytracewcs,sbt,sbp,sne,imsize=[512,512],modelid=54,scchead=hdrevent,losnbp=2,losrange=[-20,20],$
modparam=mp,/usedefault

; -- replace parameters with those from fit
mp[1]=0.4
mp[2]=4.

; ---- generate image of model corresponding to the event
rtraytracewcs,sbt,sbp,sne,imsize=[512,512],modelid=54,scchead=hdrevent,losnbp=128,losrange=[-15,15],$
modparam=mp,neang=[80.,13.,60]*!dtor

; ---- display the result
wnd,1,alog10(sbt.im > 1e-11)

Figure 4.1. Simulation of a CME as seen from SECCHI COR2-A

Simulation of a CME as seen from SECCHI COR2-A
Simulation of a CME as seen from SECCHI COR2-A


1.2.2. Simulation of a comet nucleus using WCS and SPICE

In this example, we show how to position a celestial body in space and simulate how it is viewed from a given instrument. In this case, we use comet Encke, which was visible in SECCHI HI1-A images, like in 20070425_121000_s4h1A.fts for example.

To access the ephemeris of comet Encke, we use a SPICE kernel that can be downloaded at the NAIF web site. The code sequence needed to initialized the SPICE software, load the Encke ephemeris kernel and calculate the position of its nucleus in Carrington coordinates is shown bellow:

load_stereo_spice, errmsg=message, _extra=_extra
; ---- load the solar system small body ephemeris kernel
cspice_furnsh,'/user/albert/comets/SPK/wld4826.15'
; ---- convert image data image date into SPICE compatible format
utc = anytim2utc(date, /ccsds)
; ---- define the frames
origin_to='Sun' 
frame_from='IAU_SUN' 
correction='None' 
origin_from='Encke'
; ---- compute position
cspice_spkezr, origin_to, et, frame_from,correction, origin_from, origin, ltime
; ---- compute position in spherical coordinates
cspice_reclat,-origin[0:2],rad,lon,lat

; ---- format output
lonlatrad=[lon,lat,rad/oneau('km')*oneau('rsun')]
	

The model number 58 used here is a simple sphere of density that can positioned whereever needed, using the Carrington position of its center (here passed in the variable lonlatrad). So we simulate the position of the comet with a sphere of electron density and using Thomson scattering: this does not make sense from a physics point of view but this is just to illustrate how to position and render small features with the software.

; ---- prep the data image and get its header
secchi_prep,'20070425_121000_s4h1A.fts',hdrhi1a,imhi1a,/PRECOMMCORRECT_ON
; ---- set the parameters of the model
mp=[1.,0.5,lonlatrad]
; ---- simulate comet nucleus image as seen from SECCHI HI1A
rtraytracewcs,sbt,sbp,sne,imsize=[512,512],modelid=58,scchead=hdrhi1a,losnbp=500,losrange=[-130,130],modparam=mp
; ---- display
wnd,0,sne.im

Figure 4.2. Simulation of the comet Encke nucleus image as seen from HI1-A

Simulation of the comet Encke nucleus image as seen from HI1-A


1.2.3. FITS header of the simulated data

A FITS header corresponding to the simulated image can be returned by the software. This header is (normally) compliant with the FITS-WCS standard [b2006SPD370307T]. Note that an observation date must be provided in order for the program to be able to calculate the position in different coordinate systems.

1.2.4. Projection Type

The following projection types are implemented: ARC, SIN, TAN, AZP. See [b2006SPD370307T] and [b2002AA3951077].

1.2.5. Setting positions in Carrington coordinates

It is possible to set the position of the electron density and the observer using the Carrington coordinate system. hlonlat is used to set the position of the electron density model and obslonlat is used to set the position of the observer.

1.3. Electron density cube generation: buildcloud

buildcloud generates a density cube for a given model. The cube can be saved in a text file or a binnary file. The example below shows how to build a density cube for model 14. It creates a 64 x 64 x 64 cube of 60 x 60 x 60 R_Sun, the Sun center being at the center of the cube.

; ---- generate the density cube
buildcloud,14,cubesidenbpix=64L,cubesidersun=60.,outputtype=2
; ---- read the cube just created
rtreadbincube,'cube14.dat',c,szs,orig
; ---- raytrace it and display
raytracewl,sbt,imsize=[256,256],losnbp=64L,losrange=[-20,20],modelid=25,modparam=c,neang=[90,-80,30]*!dtor,/c3
wnd,0,alog10(sbt.im > 1e-14)

2. Fit of a flux rope on SECCHI data using the IDL-GUI

We present here rtsccguicloud.pro, an IDL routine that allows to fit manually the GCS CME model to an event observed by STEREO SECCHI instruments.

2.1. Demo mode

The example described hereafter can be reproduced simply by using the demo mode of rtsccguicloud:

rtsccguicloud,/demo

2.2. Preparing the data

In this example, we will use the CME event of May 15 2007. The fist step is to prepare the images with secchi_prep.

; ---- init the event image filenames
; -- A
eventtripa=['20070515_235230_s4c2A.fts','20070515_235243_s4c2A.fts','20070515_235256_s4c2A.fts']
preevtripa=['20070515_182230_s4c2A.fts','20070515_182243_s4c2A.fts','20070515_182256_s4c2A.fts']

; -- B
eventtripb=['20070515_235230_s4c2B.fts','20070515_235243_s4c2B.fts','20070515_235256_s4c2B.fts']
preevtripb=['20070515_185230_s4c2B.fts','20070515_185243_s4c2B.fts','20070515_185256_s4c2B.fts']

; -- EUVI images
euvia=['20070515_235215_n4euA.fts']
euvib=['20070515_235215_n4euB.fts']

; ---- prep the images
secchi_prep,eventtripa,hdreventa,imeventa,/polariz_on,/rotate_on,/PRECOMMCORRECT_ON,/rotinterp_on,/silent
secchi_prep,preevtripa,hdrpreeva,impreeva,/polariz_on,/rotate_on,/PRECOMMCORRECT_ON,/rotinterp_on,/silent

secchi_prep,eventtripb,hdreventb,imeventb,/polariz_on,/rotate_on,/PRECOMMCORRECT_ON,/rotinterp_on,/silent
secchi_prep,preevtripb,hdrpreevb,impreevb,/polariz_on,/rotate_on,/PRECOMMCORRECT_ON,/rotinterp_on,/silent

; -- get nice mask
ma=get_smask(hdreventa)
mb=get_smask(hdreventb)

secchi_prep,eveuvia,heuvia,imeuvia,/PRECOMMCORRECT_ON
secchi_prep,eveuvib,heuvib,imeuvib,/PRECOMMCORRECT_ON

; ---- format image for display
ima=bytscl(rebin(alog10(ma*(imeventa-impreeva) > 1e-12 < 1e-10),512,512))
imb=bytscl(rebin(alog10(mb*(imeventb-impreevb) > 1e-12 < 1e-10),512,512))

imea=alog10(rebin(imeuvia,512,512) > 1)
imeb=alog10(rebin(imeuvib,512,512) > 1)

; ---- run the GUI
rtsccguicloud,ima,imb,hdreventa,hdreventb,imeuvia=imea,hdreuvia=heuvia,imeuvib=imeb,hdreuvib=heuvib

2.3. Manual fitting of the flux rope

Figure 4.3. CME event of 2007/05/15 viewed in COR2-A and COR2-B.

Images of the CME viewed in COR2-A and COR2-B with the CME cloud overplotted.
The green cloud represent the GCS CME model.


Figure 4.4. Slider GUI tab.

Sliders GUI.

Various sliders allow to change the model parameters in order to fit the observations.


Figure 4.5. Source region position on EUVI-A and B.

Source region position on EUVI-A and B.

The source region position corresponding to the modeled CME is indicated by the star symbol. The plus signs show the orientation of the CME, and its angular extent. The points are plotted in green when they are on the visible side of the Sun, and in gray when they are located in the back side.


Figure 4.6. Cloud parameters GUI tab.

Cloud parameters GUI.

This GUI allows to change the number of points displayed in the cloud.


2.4. Generating Thomson scattering views

Figure 4.7. Model parameters GUI tab.

Model parameters GUI.

This GUI allows to set the parameters of the Thomson scattering renderer. The integration range is in units of Rsun, and its origin is at the impact distance. The LOS number of points are the number of integration points used for each LOS in the integration range: the larger that number, the slower the rendering, but the better the image quality. Finally, the electron density, and the model shell inner and outer thickness allow to set the constant electron density in the CME skin, and the thickness of the skin.


Figure 4.8. Simulated total brightness Thomson scattering images for COR2-A and COR2-B.

Simulated total brightness Thomson scattering images for COR2-A and COR2-B..


2.5. Saving the parameters and images of the fit

Different keywords of the rtsccguicloud routine allow to access the parameters of the fit, as well as the wireframe and simulated images. Here is an example of call to the program followed by a description of the different keywords.

rtsccguicloud,ima,imb,hdra,hdrb,ssim=ssim,sgui=sgui,swire=swire,ocout=oc

ssim

Structure containing the rendered images.

sgui

Structure containing all the parameters of the program, especially the GUI parameters.

swire

Structure containing the rendered wireframes, as well as all the parameters used to compute the rendered views.

ocout

Contains the 3D point positions of the wireframe.

2.6. Automatic fit of the model position.

It is posible to determine semi-automaticaly the position and direction of the CME, assuming the GCS morphology of the CME. The first step consists in defining the contour of the observed CME leading edge. This is the purpose of the Contour tab of the GUI. Press the Draw Contour button in this tab and then start by drawing the contour of the CME seen in A, using the mouse. Once the contour is satisfiying, press the right button of the mouse to exit. You can then draw the contour of the CME viewed from B now. Press the right button to exit once done.

The Eval. Fit button computes the merit function that compares the model contour to the user drawn contour. If there is a perfect match, the merit function gives 100%, and if there is no match, it gives 0%. The result is displayed below the button and two windows are open showing the relative positions of the different contours.

The Auto Fit tab allows you to run the optimizer. Note that the optimizer plays only on the longitude, latitude and height of the model, the other parameters remain fixed. Indeed, the tilt angle, aspect ratio and half angle parameters are not sensitive enough (large deviation -> small change of the merit function) so the optimizer will not converge quickly enough and/or the solution will not be unique.

2.7. Sensitivity analysis.

The tab Sensit. allows you to perform a sensitivity analysis on the different model parameters.

Figure 4.9. GUI tab used to perform sensitivity analysis.

GUI tab used to perform sensitivity analysis.


The edit box allows you to enter the range of variation for the analysis. In the example of Figure 4.9, “GUI tab used to perform sensitivity analysis.”, it is set to plus or minus 10%. Each following buttons allows you to run the analysis for each of the model parameters. In the example, it has been ran for the longitude and the latitude. The outcome of the analysis are given below. The first and second displayed numbers are the deviation, in units of the corresponding parameter. In the case of the longitude and the latitude it is then in degrees. The two following numbers are the number of steps. The step is hard wired in the software, it cannot be changed, for now. For the longitude and latitude, the step is set to 0.1 degrees. The two following numbers are the value of the merit function, in each direction. It should be equal aproxymately to the max of the merit function minus the percentage set in the sensitivity range. Finally, the last number shows the sensitivity range, which in this case is 10%.

3. Software concept overview

3.1. Units

SCRaytrace measures the distances in Solar radii, and the angles in radians. One solar radius is taken to be 6.955x108 meters.

3.2. Coordinate Systems

All the coordinate systems are direct and orthonormal. Figure 4.10, “Overview of the three coordinate systems used in SCRaytrace. ” shows a summary of the 3 coordinate systems used in the software. ABS is the absolute coordinate system. It is placed at the center of the Sun, represented by the yellow circle here. NPS is the coordinate system of the model and OBS is the one for the observer.

Figure 4.10. Overview of the three coordinate systems used in SCRaytrace.

Coordinate systems.

3.2.1. Abolute coordinate system: ABS

This is the coordinate system that is the reference for all the other systems described bellow: translation and rotation of the other coordinate systems will be calculated with respect to this one. The origin is at the center of the Sun.

Ox

Points to Solar North.

Oy

Points to Carrington longitude and latitude (270°,0).

Oz

Points to Carrington longitude and latitude (0,0).

3.2.2. Model coordinate system: NPS

Defines the position and orientation of a model, like an electron density cube for example, with respect to the absolute coordinate system ABS.

3.2.3. Observer coordinate system: OBS

Defines the position and orientation of the observer with respect to the absolute coordinate system. We use the convention of optics for the axis orientation:

Ox

Vertical axis of the camera detector.

Oy

Horizontal axis of the camera detector.

Oz

Optical axis of the telescope, going out of the camera.

3.3. The Sun

No spectrum is used and its radiance is supposed to be unity. For the Thomson scattering physics the default limb darkening coefficient is set to 0.58, but it can be overwritten by the user.

3.4. Observer

The main intrinsic parameters of the observing instrument are the size in pixel of the virtual detector and the angular resolution of one pixel or plate-scale. A projection type can also be chosen. The position and the orientation in space of the instrument also have to be defined.

3.5. Model Representation

The software can use two types of model representations: geometrical functions and density cubes.

3.5.1. Geometric Model

The model is defined by a geometrical function that gives the electron density (in the case of Thomson Scattering physics for example) for a given coordinate position x,y,z of space (Equation 4.1, “Geometric Model”). The geometrical function generally will have some extra parameters, noted p, in order to be able modify the shape and density of the model for example.

Equation 4.1. Geometric Model

Ne = f(x,y,z,p)

3.5.2. Density Cube

The electron density is given by a three-dimensional density cube. Each voxel of the cube contains the value of the electron density for that region of space. The advantage of the density cube is that it can contain the representation of a complex structure which could not be simply defined by a geometrical function. On the other hand, if a higher resolution is needed, to represent the fine details of a coronal structure for example, the amount of memory necessary to store the density cube could be a problem.

3.6. Physical process

Thomson scattering physics needs the electron density in electrons.cm-3. For the other types of physics, please refer to the source code documentation available at http://secchi.nrl.navy.mil/synomaps/scraytrace/doxy/annotated.html. The name of the classes for the various physical processes implemented are PhysicsName, with Name the name of the physics.

3.6.1. Thomson Scattering

The K corona component we observe in white light is the result of Thomson scattering of the photospheric light by the electrons of the corona [Minaert], [Van De Hulst], [Billings]. That light is polarized. The equations Equation 4.2, “Total brightness” and Equation 4.3, “Polarized brightness” give the total and polarized brightness scattered from the solar photosphere by a localized electron density and integrated along the lines of sight (LOS) of an observer. SCRaytrace provides an implementation of these equations.

Equation 4.2. Total brightness


Equation 4.3. Polarized brightness