1 tplarson 1.1 Visualizing the Sun:
2 An Introduction to the SoshPy Visualization Package
3
4 a more extensive discussion of helioseismology with added graphics can
5 be found at http://solar-center.stanford.edu/SoSH/ . short instructions
6 can be found in quickstart_visual.txt .
7
8 the study of oscillations inside the sun is called helioseismology. in
9 particular, we are here concerned with acoustic waves for which the sun
10 provides a resonant cavity. the Sonification of Solar Harmonics (SoSH)
11 Project is a collaborative research initiative that seeks to transform
12 helioseismology into a listening experience for the scientist and
13 nonscientist alike. it has developed the SoSH Tool, which is able to
14 sonify solar acoustic data by filtering and transposing it up to the
15 range of human hearing. the SoshPy Visualization Package was created to
16 provide complementary graphical representations of the solar harmonics
17 being sonified.
18
19 in what follows we provide a brief discussion of helioseismology, which
20 augments the discussion in the documentation provided with the SoSH Tool
21 (instructions_audio.txt). we then move on to a detailed description of
22 tplarson 1.1 the included python module and the various scripts that use it.
23
24
25 Section 1: Helioseismology
26
27 consider that we are here concerned with sound, which is simply the
28 coherent vibration of matter. now imagine a plasma element (a chunk of
29 matter) at the surface of the sun. in response to a given mode, which
30 inhabits the entire sun, this single plasma element will oscillate at
31 the mode frequency. in other words, it moves back and forth; some part
32 of this motion will be in the radial direction, some in the latitudinal
33 direction, and some in the longitudinal direction. a model can tell us
34 how much of the motion lies in each direction. we call the radial
35 motion the vertical displacement, and we put the latitudinal and
36 longitudinal components of the motion together to make the horizontal
37 displacement.
38
39 now, how are we to measure this displacement? actually, it turns out to
40 be much easier to measure the speed of our plasma element, which is
41 simply the rate at which the displacment is changing. in mathematical
42 terms, we say the velocity of the plasma element is the derivative of
43 tplarson 1.1 its position. since the plasma element is moving in response to one of
44 the sun's harmonics, it can be modeled as a simple harmonic oscillator.
45 in other words, its motion is mathematically equivalent to the motion of
46 a weight attached to a spring. in particular, the model tells us that
47 the magnitude of the the plasma element's velocity is equal to the
48 magnitude of its displacement multiplied by its frequency, which is
49 simply the mode frequency.
50
51 it is the velocity of the plasma element that we are able to measure.
52 the details of this measurement are beyond the scope of this discussion,
53 but in short we are actually measuring the frequency of spectral lines
54 of atoms at the solar surface. for plasma at rest, these spectral lines
55 have precise frequencies that can be measured in a laboratory. the
56 motion of the plasma causes these frequencies to be shifted slightly.
57 this is called the doppler effect, and it is exactly the same phenomenon
58 that you can observe with sound waves when a siren drives past: when the
59 siren is approaching it sounds higher, and when it is moving away it
60 sounds lower. we use the doppler effect to create velocity images of
61 the sun's surface, and for this reason we call them dopplergrams. each
62 pixel of the image corresponds to one of our plasma elements, and the
63 pixel value tells us how fast that element is moving toward or away from
64 tplarson 1.1 the observer.
65
66 here we run into a fundamental limitation of our method of observation.
67 we would like to know both the vertical and horizontal components of the
68 velocity. unfortunately, a dopplergram only gives us the velocity along
69 the observer's line of sight. at the center of the solar image, this is
70 exactly the vertical component, and at the very edge, it is the
71 horizontal component. everywhere else we measure a mix of the two, and
72 nowhere can we measure both of them. this is one of the reasons we have
73 to rely on a model.
74
75 so what quantity should we plot to represent a mode? it turns out that
76 for a large majority of modes, the horizontal component is very small
77 compared to the vertical component at the surface, and for many purposes
78 can be neglected entirely. hence, by default we will show the vertical
79 component of velocity only. the scripts, however, are also able to plot
80 the latitudinal and longitudinal components, as well as these two
81 combined into the total horizontal component. they can also plot the
82 magnitude of the total velocity, as well this magnitude squared.
83
84 we are also able to plot either a surface view or an interior view. in
85 tplarson 1.1 the latter case, we will show a plane containing the solar rotation
86 axis. also in this case we plot the velocity scaled by the square root
87 of the background density. this is for two reasons. firstly, since the
88 velocities drop off rapidly with depth while the density increases
89 rapidly, the quantity plotted shows visible variations throughout the
90 interior. secondly, the square of this scaled velocity is actually
91 equal to the energy density of the mode.
92
93 note that for every quantity plotted, the colormap will be scaled to use
94 the full range of the data. therefore, in a surface view, the velocity
95 squared is indistinguishable from the energy density, because the
96 background mass density is constant over the solar surface. similarly,
97 for plots of single modes, we do not scale the displacement amplitudes
98 by frequency to get a velocity, because this would also be invisible.
99 when we begin adding modes together, however, we must properly account
|
101 tplarson 1.1 necessary because the frequency will be different for different modes.
102
103
104 Section 2: Surface Views - drawharmonics.py and addharmonics.py
105
106 for instructions on running these scripts, see the quickstart guide. in
107 what follows i give a more detailed account of how they work. we start
108 with drawharmonics.py, which plots surface views of single modes.
109
110 but first we need a small amount of mathematical background. an
111 oscillation mode on the sun, also known as a harmonic, is given by the
112 product of a function of latitude and longitude and another function of
113 radius only. the first function is called a spherical harmonic, and for
114 these we have explicit analytical expressions. the function of radius
115 is called a radial eigenfunction, and these must be computed numerically
116 by solving a system of different equations.
117
118 for a surface view, we only need the spherical harmonic, which is given
119 by Ylm = Plm(cos(theta))*exp(im*phi), where the Plm's are associated
120 legendre polynomials. in the coordinate system that we use, theta is
121 the angle from the solar rotation axis, so it different from the
122 tplarson 1.1 latitude by 90 degrees. the angle phi is exactly the longitude. each
123 spherical harmonic is described by two integers, l and m. the spherical
124 harmonic degree l is >= 0, and the azimuthal order m ranges from -l to
125 l. the associated legendre polynomial, however, depends on l and |m|,
126 the absolute value of m. more details can be found in
127 instructions_audio.txt .
128
129 at the beginning of the script, the arrays containing theta and phi are
130 generated by the function image2sphere(). this function is exactly the
131 same as that used to generate artificial spherical harmonics for
132 scientific analysis. it takes as arguments the number of pixels in the
133 x and y directions (default 1000) as well as three parameters describing
134 the orientation of the solar image. the first of these is bangle, which
135 is defined as the latitude of the subobserver point. in other words, it
136 is the tilt of the solar rotation axis towards the observer. similarly,
137 the parameter pangle is the tilt in the plane of the image. these
138 angels default to 30 degrees and 0 respectively. a third parameter,
139 distobs, gives the distance of the observer in solar radii. it defaults
140 to 220 (the actual distance between the sun and the earth is about 215
141 solar radii), but this can be set lower or higher to zoom in or out of
142 the solar image.
143 tplarson 1.1
144 the script then enters an input loop wherein the user is prompted to
145 enter the l, m, and n of the desired mode. this final integer, the
146 radial order n, describes the behavior of the mode along the solar
147 radius. for a surface view, it serves only to determine the ratio of
148 the vertical and horizontal components.
149
150 next we must calculate the associated legendre polynomial. this is
151 ultimately done by the function setplm(), which again is exactly the
152 same routine used for scientific analysis. to find the Plm for a given
153 l and m, we start by using an explicit expression for the Plm with l=m,
154 and then use a recursion in l to find the Plm with the requested l.
155 hence, a total of l-m Plm's are actually calculated, along with their
156 corresponding derivatives. to save on computing, these are all saved
157 for potential use in subsequent iterations of the loop.
158
159 next we construction the spherical harmonic, Ylm, as well as its
160 derivative in both the theta and phi directions. the derivative in the
161 theta direction gives the theta component, and dividing the derivative
162 in the phi direction by sin(theta) gives the phi component. the radial
163 component is given directly by the Ylm. in all cases, we must take the
164 tplarson 1.1 real part to find the physically meaningful quantity.
165
166 next, we read the vertical to horizontal ratio from a table. this table
167 was extracted from a model using the function writesurfacemodel(). the
168 table contains both the magnitudes of the vertical and horizontal
169 components of the displacement at the solar surface as well as a
170 theoretical estimate of their ratio. more details of the model are
171 provided in the next section. the script drawharmonics.py arbitrarily
172 gives the vertical component (rc) a magnitude of 100 and sets the
173 magnitude of the horizontal component (hc) according to the theoretical
174 ratio.
175
176 finally, the python dictionary varlist is populated with the six
177 plotting options:
178 Vr - radial or vertical component
179 Vt - theta or latitudinal component
180 Vp - phi or longitudinal component
181 Vh - magnitude of horizontal component, =sqrt(Vt^2+Vp^2)
182 Vmag - magnitude of total velocity, =sqrt(Vsq)
183 Vsq - magnitude squared, =(rc*Vr)^2+(hc*Vh)^2
184 note that Vr, Vt and Vp are signed quantities, while Vh, Vmag and Vsq
185 tplarson 1.1 are strictly positive. once the script terminates, all six quantities
186 will be available in varlist.
187
188 unless the script was started with animate=0, the plot will be animated.
189 to do so, it uses the two functions defined at the top of the script.
190 the first, calcylmt(i), serves to evolve the mode in time. to do so, it
191 require the arrays ylm, dylmt and dylmp to already exist. the time
192 dependence is given by exp(-iwt), but rather than using the mode's
193 actual frequency, we simply set w=2pi/nframes so that the animation will
194 loop seamlessly. nframes is the total number of frames in the animation
195 and can be set when the script is called. note that calcylmt(i) only
196 calculates the requested quantity for each frame, rather than the full
197 six which were calculated for the first frame.
198
199 the second function needed for the animation, ylmanimate(i), simply
200 takes the array returned by calcylmt(i) and uses it to set the data for
201 a previously existing image instance. this is the function that must be
202 passed to the matplotlib animation utility, and it is also called to
203 save frames to disk.
204
205 we now move on to the script addharmonics.py, which as the name suggests
206 tplarson 1.1 plots sums of modes. it is a relatively straightforward extention of
207 the previous script, but the first notable difference is that the input
208 loop does not accept the list of modes to plot. rather, this list is
209 set once at the beginning of the script. that is, the user first enters
210 the number of modes, followed by the l, m and n of the desired modes.
211 the script then enters its input loop, wherein plotting parameters may
212 be changed. in other words, to plot a new combination of modes, the
213 script must be rerun.
214
215 for each mode, the Ylm and its derivatives are calculated as before and
216 the three resulting components are stored in lists. the same table is
217 read to get the model values, but now the vertical (rc) and horizontal
218 (hc) magnitudes are used directy. we will also scale the displacement
219 amplitudes by frequency (freq) to get the proper relative velocity
220 amplitudes. these three parameters are likewise stored in lists.
221
222 next the sums are performed and the six plotting variables are stored in
223 varlist as before. now, however, when we animate the plot, we give each
224 mode in the sum its proper relative frequency. in particular, we use
225 w=2pi*freq*freqscale/nframes, where freqscale is a parameter that can be
226 set when the script is called. because the frequencies have already
227 tplarson 1.1 been scaled to be in units of millihertz, the default value of freqscale
228 is 1/3. in general, the animation will no longer loop seamlessly. to
229 see the sum evolve further in time, one may increase freqscale. one may
230 also wish to simultaneously increase nframes. note that the animation
231 functions in this script are named calcsumt(i) and sumanimate(i).
232
233 here we must insert a short discussion of colormaps. for the signed
234 quantities Vr, Vt and Vp, we would typically like for the value of zero
235 to map to the center of the colormaps. for a single mode, this is
236 almost always accomplished automatically, because the data range will
237 usually be symmetric around zero. depending upon the orientation of the
238 solar image, however, certain modes (often those with m=0) will not
239 satisfy this requirement, perhaps leading to unexpected results. to
240 correct this situation, by default we force the color mapping to be
241 symmetric around zero for the signed quantities. this means that the
242 full range of the colormap will not be used. to revert to automatic
243 scaling, one may specify colorshift=0 when starting the script. the
244 default value is 1. in either case, the color mapping is determined by
245 the data in the first frame. this is adequate for single modes. for
246 sums, however, it may happen that a subsequent frame will exceed the
247 data range of the first frame, in which case the colormap will saturate.
248 tplarson 1.1 this is usually not a problem, but if one wishes to avoid this behavior,
249 one may specify colorshift=2. in this case, the script will step
250 through calcsumt(i) for i in range(nframes) and find the full range of
251 possible data values in advance. we may also add the capability to
252 change the range as the animation proceeds, but this is not yet fully
253 implemented.
254
255
256 Section 3: Interior Views - drawradial.py and addradial.py
257
258 these two scripts operate in close analogy to their surface
259 counterparts. the most important difference is that these require a
260 much more extensive model for input. in particular, the full radial
261 eigenfunctions must be read from a file, along the mesh on which they
262 were calculated. the full eigenfunction actually consists of two
263 functions, one giving the amplitude of the displacement in the vertical
264 direction, and one giving the amplitude of the displacement in the
265 horizontal direction. we call these two functions xir and xih
266 respectively, and they each depend on both spherical harmonic degree l
267 and radial order n. finally, we will also need the background mass
268 density used to construct the model.
269 tplarson 1.1
270 all of this data has been packaged into an hdf file available from the
|
271 tplarson 1.2 SoSH website: http://solar-center.stanford.edu/SoSH/#mods . at the
272 beginning of each of these scripts, we call the function loadmodel(),
273 which reads the mesh and mass density used for the model, as well as the
274 l, n and frequency of all the modes in the model, and puts them in
275 global variables. a second function, getradial(l, n), will return the
276 actual displacement eigenfunctions for the modes we wish to plot.
277
278 after loadmodel(), the scripts call image2rtheta(), which returns the
279 arrays containing the r and theta coordinates for each point in the
280 image. like its surface counterpart, this function can take as
281 arguments the number of pixels in the x and y directions (default 1000).
282 the only other paramter it can take, however, is distobs, the observer
283 distance in solar radii. in other words, the orientation of the r-theta
284 plane is fixed. distobs should in general be set to the same value used
285 for the surface views; in both cases it defaults to 220. note, however,
286 that these scripts can take another parameter, rsurf, which is the
287 fractional image radius at which to truncate the model (default 1.0).
288 in some cases you may wish to remove the surface values if you find they
289 dominate the plots. this will reduce the size of the area occupied by
290 the plot.
291
292 tplarson 1.2 as mentioned above, the radial eigenfunctions will be scaled by the mass
293 density. for the value of phi in drawradial.py, we have arbitrarily
294 chosen pi/4m for the right half of the image, so that the real and
295 imaginary parts of the spherical harmonic (recall the factor of
296 exp(imphi)) will be the same. this means that all three components will
297 be comparable. of course, the value of phi on the left half of the
298 image will be pi/4m + pi. for m=0, we set phi=0. in addradial.py, we
299 set phi according to the first mode entered.
|
Return to
instructions_visual.txt
CVS log
Up to [Development] / JSOC / proj / globalhs / sosh