E-mail: aschwanden@lmsal.com - Markus J.Aschwanden (Lockheed Martin Solar & Astrophysics Lab.) - Updated: 2008-Dec-18

Loop 3D-Geometry Modeling - Software Tutorial

This tutorial is intended for SolarSoftWare (SSW) users who want to perform 3D-geometry fitting to coronal loops. The IDL routines are available in the package $/SSW/package/mjastereo/idl/, which are automatically included in the IDL path with the command IDL>setenv SSW_INSTR "stereo".

(1) Tracing Loops in Images

The first step in loop modeling is to define the image coordinates of an observed loop, which can be done with the IDL procedure LOOP_TRACING.PRO. The input can be an arbitrary FITS image, that has an image array IMAGE(NX,NY) and a structure INDEX containing the FITS header, as it is obtained from standard image reading routines in SSW, such as for SOHO/EIT, TRACE, STEREO/EUVI, RHESSI, or AIA/SDO image data. The image may be pre-processed with highpass filtering or be a running-time difference image to bring out the finest loop structures of interest. The image is then displayed and a few loop coordinates can be manually selected. The image can be rebinned by a factor NZOOM and displayed enlarged on the screen. The image file can be read with the routine READFITS.PRO and the FITS header can be converted into an INDEX structure:


IDL> IMAGE =READFITS(FITSFILE,H)
IDL> INDEX =FITSHEAD2STRUCT(H)
IDL> HELP,IMAGE,INDEX
IDL> NSP =7 (number of loop points)
IDL> NZOOM =8 (zooming factor)
IDL> LOOP_TRACING,IMAGE,INDEX,NSP,NZOOM,X_SUN,Y_SUN,W_SUN,FLUX
IDL> SPLINE_P,X_SUN,Y_SUN,XSP,YSP

The output coordinates of this loop are stored in the arrays X_SUN(7),Y_SUN(7) for the clicked points, and in the arrays XSP(52),YSP(52) with about 8 times higher resolution (using a 2D spline interpolation). The output coordinates are in units of solar radii with respect to Sun center.

(2) Coplanar Circular Loop Fit with Two Constrained Footpoints

We define the 3D geometry of a coplanar circular loop with 6 parameters: The parameters h0 of the loop center offset and the curvature radius r_loop are visualized in the following figure, projected into the loop plane.

The parameters of loop baseline, azimuth angle az, and loop plane inclination angle th are visualized in the following figure, projected into an arbitrary line-of-sight direction of a solar image.

For a mathematical definitions of coordinate transformations between the (1) observers coordinate system (x,y,z), (2) the heliographic coordinate system (l,b) with longitude l and latitude b, and (3) the loop plane coordinate system (X,Y,Z), see Section 3.4.4. in the textbook:

Aschwanden,M.J. 2004, Physics of the Solar Corona. An Introduction, Praxis Publishing, Chichester, UK & Springer, Berlin, New York
URL = "../eprints/2004_book/chapter_03.pdf"

We fit a circular loop model to the traced loop positions, which is defined in terms of the 6 free parameters mentioned above (BASE,H0,AZ,TH,L1,B1), but with fixed footpoints at the first [XSP(0),YSP(0)] and last loop position [XSP(NSP-1),YSP(NSP)]. This fit can be performed with the IDL procedure LOOP_MODEL1_FIT.PRO. In the following figure we plot the input of the traced loops points (XSP,YSP) with diamond symbols, the spline fit of the traced loop (XSP,YSP) with a thin solid curve, and the fit of the coplananr circular loop geometry with a thick solid linestyle. In this example, the traced 6 loop positions are given by the values XOBS and YOBS, and we use an arbitrary date and time specified in DATEOBS, which is extracted from the index of a FITS-FILE. There are 7 data points traced in this loop, which we interpolate with the routine SPLINE_P to obtain a smoothly interpolated curve with a higher resolution (of 53 data points here). The radius of the Sun and the heliographic position of the center of the solar disk (L0,B0) are computed from the ephemerides based on the input date and time given in the parameter DATE_OBS. In this example we skip the loop tracing and define the traced loop coordinates in XOBS, YOBS:


IDL>
XOBS=[-929,-999,-1054,-1075,-1039,-989,-940]/944. ;in solar radii
YOBS=[-3.,-20,-60,-117,-123,-107,-84]/944. ;in solar radii
DATEOBS='2007-06-21T04:00:00'
SPLINE_P,XOBS,YOBS,XSP,YSP
LOOP_MODEL1_FIT,DATEOBS,XSP,YSP,BASE,H0,AZ,TH,L1,B1,RSUN,L0,B0,XFIT,YFIT,ZFIT
window,0,xsize=512,ysize=512
!x.range=minmax([xobs,xfit])+0.1*[-1,+1]
!y.range=minmax([yobs,yfit])+0.1*[-1,+1]
!x.style=1
!y.style=1
plot,xobs,yobs,psym=4 ;traced loop points
oplot,xsp,ysp ;interpolated tracing
oplot,xfit,yfit,thick=3 ;best-fit loop coordinates
coordsphere,1,0,0,0,0,2.
xyouts_norm,0.05,0.9,'base='+string(base,'(F7.4)')+' solar radii',1.5
xyouts_norm,0.05,0.85,'h0='+string(h0,'(F7.4)')+' solar radii',1.5
xyouts_norm,0.05,0.80,'az='+string(az,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.75,'th='+string(th,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.70,'l1='+string(l1,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.65,'b1='+string(b1,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.15,'l0='+string(l0,'(F5.1)')+' deg',1.5
xyouts_norm,0.05,0.10,'b0='+string(b0,'(F5.1)')+' deg',1.5
xyouts_norm,0.05,0.05,'rsun='+string(rsun,'(F6.1)')+' arcsec',1.5

(3) 3D Reconstruction of loop in RHESSI Image

We read a RHESSI data cube and select image #7, trace the hard X-ray loop with 5 points, forward-fit the coplanar 3D loop geometry, and measure at the same time the loop flux and width at the spline points:


IDL>
images =readfits(rhessi_datacube,h)
index =fitshead2struct(h)
image =images(*,*,7)
dateobs =index.date_obs
nsp =5
nzoom =8
LOOP_TRACING,IMAGE,INDEX,NSP,NZOOM,FOV_SUN,X_SUN,Y_SUN,W_SUN,FLUX
SPLINE_P,X_SUN,Y_SUN,XSP,YSP
LOOP_MODEL1_FIT,DATEOBS,XSP,YSP,BASE,H0,AZ,TH,L1,B1,RSUN,L0,B0,XFIT,YFIT,ZFIT
window,0,xsize=512,ysize=512
!x.range=[fov_sun(0),fov_sun(2)]
!y.range=[fov_sun(1),fov_sun(3)]
!x.title='x [solar radii]'
!y.title='y [solar radii]'
!x.style=1
!y.style=1
!p.position=[0.1,0.1,0.9,0.9]
plot,x_sun,y_sun,psym=4 ;traced loop points
oplot,xsp,ysp ;interpolated tracing
oplot,xfit,yfit,thick=3 ;best-fit loop coordinates
dlat=1.0 ;(deg) heliographic coordinate grid
coordsphere,rsun,0,0,0,0,dlat
xyouts_norm,0.05,0.9,'base='+string(base,'(F7.4)')+' solar radii',1.5
xyouts_norm,0.05,0.85,'h0='+string(h0,'(F7.4)')+' solar radii',1.5
xyouts_norm,0.05,0.80,'az='+string(az,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.75,'th='+string(th,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.70,'l1='+string(l1,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.65,'b1='+string(b1,'(F6.1)')+' deg',1.5
xyouts_norm,0.05,0.15,'l0='+string(l0,'(F5.1)')+' deg',1.5
xyouts_norm,0.05,0.10,'b0='+string(b0,'(F5.1)')+' deg',1.5
xyouts_norm,0.05,0.05,'rsun='+string(rsun,'(F6.1)')+' arcsec',1.5
for i=0,nsp-1 do xyouts,x_sun(i),y_sun(i),$
' f='+string(flux(i),'(F4.1)')+' DN, w='+string(w_sun(i),'(F7.4)'),size=1.5





E-mail: aschwanden@lmsal.com - Markus J.Aschwanden (Lockheed Martin Solar & Astrophysics Lab.)