Table of Contents
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.
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)
Total brightness structures containing the image, the simulation parameters and a FITS header.
Polarized brightness structures containing the image, the simulation parameters and a FITS header.
Integrated Ne structures containing the image, the simulation parameters and a FITS header.
Size of the image in pixels.
LOS integration range in R_Sun, with the origin at the impact distance. Use the
frontinteg keyword to change the origin to the observer.
Number of integration points along the LOS, within the LOS range.
ID of the model.
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.
Set to the predefined LASCO C3 FOV.
Rotation order convention for the variables
neang is: 1st: Oz, 2nd: Oy, 3rd: Ox.
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
Array of parameter corresponding to the model 54. See
models51to60.cpp here for more details.
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):
x size (sx)
y size (sy)
z size (sz)
xc Sun center in pix
yc Sun center in pix
zc Sun center in pix
(0,0,0) is at the first vertice of the density cude.
voxel size in rsun, same for the 3 directions of space
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 )
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)
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
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=0.4 mp=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)
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
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.
The following projection types are implemented: ARC, SIN, TAN, AZP. See [b2006SPD370307T] and [b2002AA3951077].
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)
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.
The example described hereafter can be reproduced simply by using the demo mode of
In this example, we will use the CME event of May 15 2007. The fist step is to prepare the images with
; ---- 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
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.
Structure containing the rendered images.
Structure containing all the parameters of the program, especially the GUI parameters.
Structure containing the rendered wireframes, as well as all the parameters used to compute the rendered views.
Contains the 3D point positions of the wireframe.
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.
The tab Sensit. allows you to perform a sensitivity analysis on the different model parameters.
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%.
SCRaytrace measures the distances in Solar radii, and the angles in radians. One solar radius is taken to be 6.955x108 meters.
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.
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.
Points to Solar North.
Points to Carrington longitude and latitude (270°,0).
Points to Carrington longitude and latitude (0,0).
Defines the position and orientation of a model, like an electron density cube for example, with respect to the absolute coordinate system ABS.
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:
Vertical axis of the camera detector.
Horizontal axis of the camera detector.
Optical axis of the telescope, going out of the camera.
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.
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.
The software can use two types of model representations: geometrical functions and density cubes.
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.
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.
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
Name the name of the physics.
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.