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 
100              for their relative amplitudes.  so in this case the frequency scaling is 
101              necessary because the frequency will be different for different modes.
102              
103              
104              Section 2: Surface Views - drawharmonics.py and addharmonics.py
105              
106 tplarson 1.1 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              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 tplarson 1.1 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              
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 tplarson 1.1 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              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 tplarson 1.1 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              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 tplarson 1.1 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              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 tplarson 1.1 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              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 tplarson 1.1 
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              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 tplarson 1.1 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              
270              all of this data has been packaged into an hdf file available from the 
271              SoSH website: http://solar-center.stanford.edu/SoSH/#mods
272              
273              
274 tplarson 1.1 
275              
276