Return to instructions_visual.txt CVS log

Up to [Development] / JSOC / proj / globalhs / sosh

bug fixes

Visualizing the Sun:
An Introduction to the SoshPy Visualization Package

a more extensive discussion of helioseismology with added graphics can 
be found at http://solar-center.stanford.edu/SoSH/ .  short instructions 
can be found in quickstart_visual.txt .

the study of oscillations inside the sun is called helioseismology.  in 
particular, we are here concerned with acoustic waves for which the sun 
provides a resonant cavity.  the Sonification of Solar Harmonics (SoSH) 
Project is a collaborative research initiative that seeks to transform 
helioseismology into a listening experience for the scientist and 
nonscientist alike.  it has developed the SoSH Tool, which is able to 
sonify solar acoustic data by filtering and transposing it up to the 
range of human hearing.  the SoshPy Visualization Package was created to 
provide complementary graphical representations of the solar harmonics 
being sonified.

in what follows we provide a brief discussion of helioseismology, which 
augments the discussion in the documentation provided with the SoSH Tool 
(instructions_audio.txt).  we then move on to a detailed description of 
the included python module and the various scripts that use it.


Section 1: Helioseismology

consider that we are here concerned with sound, which is simply the 
coherent vibration of matter.  now imagine a plasma element (a chunk of 
matter) at the surface of the sun.  in response to a given mode, which 
inhabits the entire sun, this single plasma element will oscillate at 
the mode frequency.  in other words, it moves back and forth; some part 
of this motion will be in the radial direction, some in the latitudinal 
direction, and some in the longitudinal direction.  a model can tell us 
how much of the motion lies in each direction.  we call the radial 
motion the vertical displacement, and we put the latitudinal and 
longitudinal components of the motion together to make the horizontal 
displacement.

now, how are we to measure this displacement?  actually, it turns out to 
be much easier to measure the speed of our plasma element, which is 
simply the rate at which the displacment is changing.  in mathematical 
terms, we say the velocity of the plasma element is the derivative of 
its position.  since the plasma element is moving in response to one of 
the sun's harmonics, it can be modeled as a simple harmonic oscillator.  
in other words, its motion is mathematically equivalent to the motion of 
a weight attached to a spring.  in particular, the model tells us that 
the magnitude of the the plasma element's velocity is equal to the 
magnitude of its displacement multiplied by its frequency, which is 
simply the mode frequency.

it is the velocity of the plasma element that we are able to measure.  
the details of this measurement are beyond the scope of this discussion, 
but in short we are actually measuring the frequency of spectral lines 
of atoms at the solar surface.  for plasma at rest, these spectral lines 
have precise frequencies that can be measured in a laboratory.  the 
motion of the plasma causes these frequencies to be shifted slightly.  
this is called the doppler effect, and it is exactly the same phenomenon 
that you can observe with sound waves when a siren drives past: when the 
siren is approaching it sounds higher, and when it is moving away it 
sounds lower.  we use the doppler effect to create velocity images of 
the sun's surface, and for this reason we call them dopplergrams.  each 
pixel of the image corresponds to one of our plasma elements, and the 
pixel value tells us how fast that element is moving toward or away from 
the observer.

here we run into a fundamental limitation of our method of observation. 
we would like to know both the vertical and horizontal components of the 
velocity.  unfortunately, a dopplergram only gives us the velocity along 
the observer's line of sight.  at the center of the solar image, this is 
exactly the vertical component, and at the very edge, it is the 
horizontal component.  everywhere else we measure a mix of the two, and 
nowhere can we measure both of them.  this is one of the reasons we have 
to rely on a model.

so what quantity should we plot to represent a mode?  it turns out that 
for a large majority of modes, the horizontal component is very small 
compared to the vertical component at the surface, and for many purposes 
can be neglected entirely.  hence, by default we will show the vertical 
component of velocity only.  the scripts, however, are also able to plot 
the latitudinal and longitudinal components, as well as these two 
combined into the total horizontal component.  they can also plot the 
magnitude of the total velocity, as well this magnitude squared.

we are also able to plot either a surface view or an interior view.  in 
the latter case, we will show a plane containing the solar rotation 
axis.  also in this case we plot the velocity scaled by the square root 
of the background density.  this is for two reasons.  firstly, since the 
velocities drop off rapidly with depth while the density increases 
rapidly, the quantity plotted shows visible variations throughout the 
interior.  secondly, the square of this scaled velocity is actually 
equal to the energy density of the mode.

note that for every quantity plotted, the colormap will be scaled to use 
the full range of the data.  therefore, in a surface view, the velocity 
squared is indistinguishable from the energy density, because the 
background mass density is constant over the solar surface.  similarly, 
for plots of single modes, we do not scale the displacement amplitudes 
by frequency to get a velocity, because this would also be invisible.  
when we begin adding modes together, however, we must properly account 
for their relative amplitudes, so in this case the frequency scaling is 
necessary because the frequency will be different for different modes.


Section 2: Surface Views - drawharmonics.py and addharmonics.py

for instructions on running these scripts, see the quickstart guide.  in 
what follows i give a more detailed account of how they work.  we start 
with drawharmonics.py, which plots surface views of single modes.

but first we need a small amount of mathematical background.  an 
oscillation mode on the sun, also known as a harmonic, is given by the 
product of a function of latitude and longitude and another function of 
radius only.  the first function is called a spherical harmonic, and for 
these we have explicit analytical expressions.  the function of radius 
is called a radial eigenfunction, and these must be computed numerically 
by solving a system of different equations.

for a surface view, we only need the spherical harmonic, which is given 
by Ylm = Plm(cos(theta))*exp(im*phi), where the Plm's are associated 
legendre polynomials.  in the coordinate system that we use, theta is 
the angle from the solar rotation axis, so it different from the 
latitude by 90 degrees.  the angle phi is exactly the longitude.  each 
spherical harmonic is described by two integers, l and m.  the spherical 
harmonic degree l is >= 0, and the azimuthal order m ranges from -l to 
l.  the associated legendre polynomial, however, depends on l and |m|, 
the absolute value of m. more details can be found in 
instructions_audio.txt .

at the beginning of the script, the arrays containing theta and phi are 
generated by the function image2sphere().  this function is exactly the 
same as that used to generate artificial spherical harmonics for 
scientific analysis.  it takes as arguments the number of pixels in the 
x and y directions (default 1000) as well as three parameters describing 
the orientation of the solar image.  the first of these is bangle, which 
is defined as the latitude of the subobserver point.  in other words, it 
is the tilt of the solar rotation axis towards the observer.  similarly, 
the parameter pangle is the tilt in the plane of the image.  these 
angels default to 30 degrees and 0 respectively.  a third parameter, 
distobs, gives the distance of the observer in solar radii.  it defaults 
to 220 (the actual distance between the sun and the earth is about 215 
solar radii), but this can be set lower or higher to zoom in or out of 
the solar image.

the script then enters an input loop wherein the user is prompted to 
enter the l, m, and n of the desired mode.  this final integer, the 
radial order n, describes the behavior of the mode along the solar 
radius.  for a surface view, it serves only to determine the ratio of 
the vertical and horizontal components.

next we must calculate the associated legendre polynomial.  this is 
ultimately done by the function setplm(), which again is exactly the 
same routine used for scientific analysis.  to find the Plm for a given 
l and m, we start by using an explicit expression for the Plm with l=m, 
and then use a recursion in l to find the Plm with the requested l.  
hence, a total of l-m Plm's are actually calculated, along with their 
corresponding derivatives.  to save on computing, these are all saved 
for potential use in subsequent iterations of the loop.

next we construction the spherical harmonic, Ylm, as well as its 
derivative in both the theta and phi directions.  the derivative in the 
theta direction gives the theta component, and dividing the derivative 
in the phi direction by sin(theta) gives the phi component.  the radial 
component is given directly by the Ylm.  in all cases, we must take the 
real part to find the physically meaningful quantity.

next, we read the vertical to horizontal ratio from a table.  this table 
was extracted from a model using the function writesurfacemodel().  the 
table contains both the magnitudes of the vertical and horizontal 
components of the displacement at the solar surface as well as a 
theoretical estimate of their ratio.  more details of the model are 
provided in the next section.  the script drawharmonics.py arbitrarily 
gives the vertical component (rc) a magnitude of 100 and sets the 
magnitude of the horizontal component (hc) according to the theoretical 
ratio.

finally, the python dictionary varlist is populated with the six 
plotting options: 
Vr - radial or vertical component 
Vt - theta or latitudinal component 
Vp - phi or longitudinal component 
Vh - magnitude of horizontal component, =sqrt(Vt^2+Vp^2) 
Vmag - magnitude of total velocity, =sqrt(Vsq) 
Vsq - magnitude squared, =(rc*Vr)^2+(hc*Vh)^2 
note that Vr, Vt and Vp are signed quantities, while Vh, Vmag and Vsq 
are strictly positive.  once the script terminates, all six quantities 
will be available in varlist.

unless the script was started with animate=0, the plot will be animated.  
to do so, it uses the two functions defined at the top of the script.  
the first, calcylmt(i), serves to evolve the mode in time.  to do so, it 
require the arrays ylm, dylmt and dylmp to already exist.  the time 
dependence is given by exp(-iwt), but rather than using the mode's 
actual frequency, we simply set w=2pi/nframes so that the animation will 
loop seamlessly. nframes is the total number of frames in the animation 
and can be set when the script is called.  note that calcylmt(i) only 
calculates the requested quantity for each frame, rather than the full 
six which were calculated for the first frame.

the second function needed for the animation, ylmanimate(i), simply 
takes the array returned by calcylmt(i) and uses it to set the data for 
a previously existing image instance.  this is the function that must be 
passed to the matplotlib animation utility, and it is also called to 
save frames to disk.

we now move on to the script addharmonics.py, which as the name suggests 
plots sums of modes.  it is a relatively straightforward extention of 
the previous script, but the first notable difference is that the input 
loop does not accept the list of modes to plot.  rather, this list is 
set once at the beginning of the script.  that is, the user first enters 
the number of modes, followed by the l, m and n of the desired modes.  
the script then enters its input loop, wherein plotting parameters may 
be changed.  in other words, to plot a new combination of modes, the 
script must be rerun.

for each mode, the Ylm and its derivatives are calculated as before and 
the three resulting components are stored in lists.  the same table is 
read to get the model values, but now the vertical (rc) and horizontal 
(hc) magnitudes are used directy.  we will also scale the displacement 
amplitudes by frequency (freq) to get the proper relative velocity 
amplitudes.  these three parameters are likewise stored in lists.

next the sums are performed and the six plotting variables are stored in 
varlist as before.  now, however, when we animate the plot, we give each 
mode in the sum its proper relative frequency.  in particular, we use 
w=2pi*freq*freqscale/nframes, where freqscale is a parameter that can be 
set when the script is called.  because the frequencies have already 
been scaled to be in units of millihertz, the default value of freqscale 
is 1/3.  in general, the animation will no longer loop seamlessly.  to 
see the sum evolve further in time, one may increase freqscale.  one may 
also wish to simultaneously increase nframes.  note that the animation 
functions in this script are named calcsumt(i) and sumanimate(i).

here we must insert a short discussion of colormaps.  for the signed 
quantities Vr, Vt and Vp, we would typically like for the value of zero 
to map to the center of the colormaps.  for a single mode, this is 
almost always accomplished automatically, because the data range will 
usually be symmetric around zero.  depending upon the orientation of the 
solar image, however, certain modes (often those with m=0) will not 
satisfy this requirement, perhaps leading to unexpected results.  to 
correct this situation, by default we force the color mapping to be 
symmetric around zero for the signed quantities.  this means that the 
full range of the colormap will not be used.  to revert to automatic 
scaling, one may specify colorshift=0 when starting the script.  the 
default value is 1.  in either case, the color mapping is determined by 
the data in the first frame.  this is adequate for single modes.  for 
sums, however, it may happen that a subsequent frame will exceed the 
data range of the first frame, in which case the colormap will saturate.  
this is usually not a problem, but if one wishes to avoid this behavior, 
one may specify colorshift=2.  in this case, the script will step 
through calcsumt(i) for i in range(nframes) and find the full range of 
possible data values in advance.  we may also add the capability to 
change the range as the animation proceeds, but this is not yet fully 
implemented.


Section 3: Interior Views - drawradial.py and addradial.py

these two scripts operate in close analogy to their surface 
counterparts.  the most important difference is that these require a 
much more extensive model for input.  in particular, the full radial 
eigenfunctions must be read from a file, along the mesh on which they 
were calculated.  the full eigenfunction actually consists of two 
functions, one giving the amplitude of the displacement in the vertical 
direction, and one giving the amplitude of the displacement in the 
horizontal direction.  we call these two functions xir and xih 
respectively, and they each depend on both spherical harmonic degree l 
and radial order n.  finally, we will also need the background mass 
density used to construct the model.

all of this data has been packaged into an hdf file available from the 
SoSH website: http://solar-center.stanford.edu/SoSH/#mods .  at the 
beginning of each of these scripts, we call the function loadmodel(), 
which reads the mesh and mass density used for the model, as well as the 
l, n and frequency of all the modes in the model, and puts them in 
global variables.  a second function, getradial(l, n), will return the 
actual displacement eigenfunctions for the modes we wish to plot.

after loadmodel(), the scripts call image2rtheta(), which returns the 
arrays containing the r and theta coordinates for each point in the 
image.  like its surface counterpart, this function can take as 
arguments the number of pixels in the x and y directions (default 1000).  
the only other paramter it can take, however, is distobs, the observer 
distance in solar radii.  in other words, the orientation of the r-theta 
plane is fixed.  distobs should in general be set to the same value used 
for the surface views; in both cases it defaults to 220.  note, however, 
that these scripts can take another parameter, rsurf, which is the 
fractional image radius at which to truncate the model (default 1.0).  
in some cases you may wish to remove the surface values if you find they 
dominate the plots.  this will reduce the size of the area occupied by 
the plot.

as mentioned above, the radial eigenfunctions will be scaled by the mass 
density.  for the value of phi in drawradial.py, we have arbitrarily 
chosen pi/4m for the right half of the image, so that the real and 
imaginary parts of the spherical harmonic (recall the factor of 
exp(imphi)) will be the same.  this means that all three components will 
be comparable.  of course, the value of phi on the left half of the 
image will be pi/4m + pi.  for m=0, we set phi=0.  in addradial.py, we 
set phi according to the first mode entered.