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.