VIEWS: 14 PAGES: 109


            From left: Creighton Wilson (RET, Lovelady High School), Michelle McMillan (REU, Northern
            Arizona Univ.), Drew Medlin (Grad SRA, New Mexico Tech), Leah Simon (Grad SRA,
            Macalester College), Maria Kazachenko (Undergrad SRA, St. Petersburg Univ., Russia),
            Michael Sinclair (RET, Kalamazoo Math and Science Center), Brian Robinson (ATST Fellow,
            Univ. of Alabama-Huntsville), Stuart Robbins (REU, Case Western Reserve Univ.), Brian
            Harker-Lundberg (Grad SRA, Utah State Univ.), Frances Edelman (REU, Yale), Statia Luszcz
            (in front of Frances–REU, Cornell), Joel Lamb (REU, Univ. of Iowa), Mark Calhoun (RET,
            Sabino High School), and Heidi Gerhardt (in front–REU, Towson Univ.). Missing from photo are
            Matt Dawson (RET, Brockton High School), Cheryl-Annette Kincaid (Undergrad SRA, Univ. of
            North Texas), Kimberly Moody (Undergrad SRA, Univ. of Arizona), Anna Malanushenko (Under-
            grad SRA, St. Petersburg Univ. Russia).

2004 Annual Program Report
28 January 2005

             The NSO REU site program is co-funded by the Department of Defense
                           in partnership with the NSF REU Program
          NSO is operated by the Association of Universities for Research in Astronomy
                under cooperative agreement with the National Science Foundation
                                                              TABLE OF CONTENTS

Project Summary – 2004 NSO REU Site Program................................................................................................ 1
Participants ..............................................................................................................................................................1
2004 Program Activities ............................................................................................................................................2
Research Projects ........................................................................................................................................................2
Scientific Lectures and Field Trips ...........................................................................................................................4
APPENDIX I — NSO 2004 Research Experience for Teachers (RET) Program................................ 8
      •      Mark A. Calhoun “Imaging and Identification of Geostationary Satellites”........................................ 8
      •      Matthew Dawson “Nowcasting Air Temperature at the Advanced Technology
             Solar Telescope (ATST)”..........................................................................................................................16
      •      Michael Sinclair “Deconvolving the Point Spread Function from Solar Mass Ejection Imagery” ....33
      •      Creighton G. Wilson “The Neutral Line Method of Flare Prediction from Magnetograms” ............ 44
APPENDIX II – REU Student Reports .................................................................................................. 51
      •      Frances S. Edelman “The Effect of Spectral Line Shape Changes on GONG Observations
             of Oscillations and Flares”......................................................................................................................52
      •     Heidi H. Gerhardt “Understanding the Mechanism of the Penumbra and the
            Curled Light Bridge ” ................................................................................................................................55
      •      Joel B. Lamb “A Study of Plasma Flows in Emerging Active Regions”................................................63
      •      Statia H. Luszcz “Magnetic Field Measurements and Moving Magnetic Features (MMFs)
             Using the Infrared Fe Line at 1564.8 nm”..............................................................................................71
      •     Michelle T. McMillan “On the Accuracy of Polarimetric Inversion Codes for Spectroscopic
            Data Analysis” .........................................................................................................................................77
      •      Stuart J. Robbins “Solar Wind Forecasting with Coronal Holes” .......................................................85

i NSO REU/RET Annual Program Report 2004 – Table of Contents
Project Summary
The National Science Foundation (NSF) Research Experiences for Undergraduates (REU) program is
designed to encourage US college and university undergraduate students to pursue careers in science and
engineering. REU site programs (funded through the NSF Division of Astronomy) make it possible for
undergraduates to take part in independent research activities with professional scientists at major research
institutions, usually during the summer months. The National Solar Observatory (NSO) has participated in
the NSF REU program since inception of the program in 1986, and our tradition of contributing to the
program’s continued success has been rewarding. Several of our past REU and summer research assistants
are leaders in astronomy and astrophysics organizations in the US and abroad.

The 2004 NSO REU Program provided an opportunity for six US undergraduate students to participate in
astronomical research activities at the NSO sites in Tucson, Arizona and Sacramento Peak in Sunspot, New
Mexico. Three students were in residence at Sacramento Peak and three students were in Tucson. The
undergraduates’ experiences were enhanced by the presence of three NSO-supported summer graduate
students at Sac Peak and four undergraduate research assistants supported by NSO grant-funded and the US
Air Force undergraduate research programs (two each in Tucson and at Sac Peak). Another significant
aspect of the educational component of the NSO summer program was the participation of four high school
teachers supported through the NSF Research Experiences for Teachers (RET) program, which is funded as a
supplement to the NSO REU grant. The annual report for the Summer 2004 NSO RET program is presented
in Appendix I.

The success of the REU/RET Program at both NSO sites can be attributed to dedicated mentors and an
atmosphere conducive to research education where students interact with each other, with researchers,
engineers and visiting scientists at NSO and NOAO, and with similar students drawn from other programs
such as the NSO Summer Research Assistant Program and PhD students drawn from other US and foreign

The six undergraduate students (four women and two men) were selected from a total of 54 applicants.
Students were recruited for the 2004 program through the distribution of posters and application forms to a
broad spectrum of colleges and universities and through the NSO REU Web site (
More than 700 announcements were sent to astronomy, physics, engineering, mathematics, and natural
science departments throughout the United States and Puerto Rico. In an effort to attract students from
underrepresented areas, our mailing list includes: colleges from the Historically Black College List generated
by the NSF and a list of American Indian Science and Engineering Society Affiliates. An announcement
about the NSO Summer 2004 RET Program was posted on the NSO RET Web site (
and widely distributed electronically via the National Astronomy List Serves, Arizona Physics Teachers List
Serves, and the Arizona Science and Math Education Center; announcements were also sent to school
districts throughout New Mexico and in Tucson. Scientific program coordinators for each NSO site were
designated to oversee the selection of REU and RET project advisors, students, and teachers. Staff members
interested in participating in the REU/RET program developed research projects, which were reviewed for
scientific merit and quality of student involvement. Proposals that involved two or more students, or two or
more teachers, collaborating on a project as part of a team were encouraged.

1 NSO REU/RET Annual Program Report 2004                                          Project Summary & Participants
The REU students, their college/university affiliations, and the NSO staff designated as their research
advisors for the program are listed here.

Student                      Advisor(s)                     College/University                    NSO Site

Frances S. Edelman           Frank Hill                     Yale University                       Tucson
Heidi H. Gerhardt            K. Sankarasubramanian          Towson University                     Sunspot
Joel B. Lamb                 Alexei Pevtsov                 University of Iowa                    Sunspot
Statia H. Luszcz             Matt Penn                      Cornell University                    Tucson
Michelle T. McMillan         Han Uitenbroek                 Northern Arizona University           Sunspot
Stuart J. Robbins            Carl Henney & Jack Harvey      Case Western Reserve University       Tucson

2004 Program Activities
Research Projects
The Summer 2004 REU students at NSO spent an average of 10 weeks working as full-time research
assistants. They were involved in a wide variety of research projects, the results of which are described in
their reports in Appendix II. In addition to work on their individual projects, the Tucson students were
given opportunities to observe at the McMath-Pierce solar telescope on Kitt Peak and collaborated in
designing observational programs carried out at the Kitt Peak National Observatory (KPNO) 2.1-meter
telescope. The Sac Peak students participated in observational programs at the Dunn Solar Telescope and the
Evans Solar Facility. Toward the end of the program, each student presented a seminar on his or her research
work and discussed the results with the staff. The seminar is in a format similar to an AAS meeting,
providing students with experience in presenting research results in a professional scientific setting. The
informal discussion following each presentation is designed to give the students positive feedback and
encouragement to continue pursuing research. A brief description of each student’s project follows, and
copies of the students’ final reports are provided in Appendix II.

Frances S. Edelman (Yale University) computed the response of the Global Oscillation Network Group
(GONG) instrument to magnetically induced spectral line changes. Her work demonstrated that the flare-
associated magnetic field changes recently observed by GONG are not due to contamination of the measure-
ments and are therefore actual solar changes. In addition, the study showed that a substantial fraction of the
well-known decrease of the 5-minute oscillation amplitude in active regions is due to magnetically-induced
changes in the spectral line shape. Frances' work was presented at a GONG/SOHO meeting at Yale in July
2004, and has been published in the conference proceedings (Edelman, F., Hill, F., Howe, R., and Komm, R.
2004, in Helio- and Asteroseismology: Towards a Golden Future, SOHO14-GONG2004, ESA SP-559, ed. D.
Danesy, 416-419: The effect of spectral line shape changes on GONG observations of oscillations and flares).
Frances’ advisor was Frank Hill.

Heidi H. Gerhardt (Towson University), with advisor K. Sankarasubramanian (Sankar), worked with high
spatial resolution spectro-polarimetric data taken with the diffraction-limited spectro-polarimeter (DLSP) at
the Dunn Solar Telescope. She first took the time to understand the basic data calibration procedure for the
DLSP, then reduced and calibrated the DLSP data for a sunspot. Heidi worked on different inversion
procedures to obtain the vector magnetic fields and other physical parameters of the sunspot. The simple

2 NSO REU/RET Annual Program Report 2004                                                       Program Activities
Milne-Eddington inversion code was performed on the sunspot data. This inversion code doesn't invert
abnormal profiles, which were found in the penumbral regions of the sunspot and also at the light bridges.
To invert those profiles, Heidi used the Stokes inversion based on the Response function (SIR) code. This
code uses the line-of-sight gradients of the physical parameters for the abnormal profiles. The main result
from her analysis is that there are large downflows present at the photospheric layers in the penumbral region
which produce the large asymmetry or the abnormal profiles in the sunspot. Analysis on the light bridge
revealed that convective flows were present in these areas.

Joel B. Lamb (University of Iowa), with advisor Alexei Pevtsov, studied plasma flows in and around
emerging solar active regions using observations with the SOHO Michelson Doppler Imager (MDI). The
goal was to search for flow patterns in the photosphere before and during emergence of a new active region.
Main results indicate that a) there is no persistent pattern of downflows in the photosphere prior to
emergence of an active region. Downflows develop only after the beginning of a new flux emergence; b) at
very early stages of emergence, some active regions show plasma flowing from one polarity to the other.
These flows were predicted by theory. However, the direction of flows is not always in agreement with the
theoretical prediction; and c) what was observed, in several cases, was a formation of sunspot penumbra prior
to Evershed flow development. Joel presented the results of this work as a poster paper, "A Study of
Plasma Flows in Emerging Active Regions," at the January 2005 meeting of the American Astronomical
Society in San Diego.

Statia H. Luszcz (Cornell University), with advisor Matt Penn, examined a time sequence of vector magneto-
grams of a large sunspot region taken at the McMath-Pierce Solar Telescope during the summer of 2002. In
particular, Statia looked for small moving magnetic features (MMFs), as seen in the line-of-sight magnetic
field and also in the magnetic field inclination angle (relative to the line of sight). She found that of the 30 or
so clearly identified MMF objects seen in this data set, about half began their motions well inside the sunspot
penumbra. This is a completely new result, different from all previous studies of MMF objects! It is thought
that because the observations used an infrared spectral line, they are seeing MMF objects much deeper in the
solar atmosphere than the studies which use visible spectral lines. These new, deep, penumbral MMF
features currently are not explained by theory. Statia will present her results at the May 2005 AAS Solar
Physics Division Meeting in New Orleans.

Michelle T. McMillan (Northern Arizona University), together with advisors K. Sankarasubramanian and
H. Uitenbroek, investigated the reliability of spectral inversion codes for determining the strength of
magnetic fields on the solar surface. Such inversion codes are commonly used to retrieve values of physical
parameters from observed spectra. Several codes in use by the community were tested with a forward-
backward approach. In this method a spectrum is generated from a known solar model including a magnetic
field. The calculated spectrum is then fed to the inversion codes. A comparison of the results of this
inversion with the known values from the original model subsequently yields a good impression of the
reliability and accuracy of the inversion codes. Findings of these comparisons are that the values produced
by inversion codes should be interpreted with utmost care, as there is often a large spread between retrieved
values and actual values in the simulations. Since these simulations are highly realistic, the same can be
expected from actual observations. One of the problems identified is the structure of the solar surface
magnetic field with many changes in gradients along a typical line-of-sight; this is more complex than the
assumptions allowed in the inversion codes. Michelle presented her results as a poster paper, “Can
Photospheric Stokes Profiles Be Reliably Inverted?” at the January 2005 meeting of the American
Astronomical Society in San Diego.

3 NSO REU/RET Annual Program Report 2004                                                          Program Activities
Stuart J. Robbins (Case Western Reserve University), with advisors Carl Henney and Jack Harvey, developed
an empirical model for forecasting solar wind driven geomagnetic events. The forecasting model was derived
by cross-correlating Kitt Peak Vacuum Telescope coronal hole source maps with the solar wind velocity time
series provided by OMNIWeb. Using coronal hole size and location information, the model was found to
predict solar wind velocity up to 8.5 days in advance with a statistically significant correlation of 0.38 for the
11-year period investigated, from May 1992 through September 2003. In addition, Robbins developed a
prototype Web-based display for future solar wind forecasting using SOLIS vector spectromagnetograph
data. Robbins is preparing a paper to be submitted to Solar Physics summarizing his findings. He will also
present his results at the 205th Meeting of the AAS in San Diego in January 2005 and at May 2005 AAS Solar
Physics Division Meeting in New Orleans.

Scientific Lectures and Field Trips
In addition to their work on projects, students at both sites were exposed to a broad range of research topics
through a series of informal talks by scientific, engineering and technical staff, seminars by visiting
astronomers (both solar and non-solar), and through reports by the site staff on their ongoing projects. The
students also participated in field trips to nearby observatories and non-NSO facilities. Along with the
students in the concurrent REU program at the Kitt Peak National Observatory (KPNO), the Tucson
students traveled to Sunspot on a four-day field trip to visit the National Radio Astronomy Observatory’s
Very Large Array (NRAO/VLA) and the NSO facilities at Sac Peak. Likewise, the Sac Peak students traveled
to Tucson to visit the NSO and NOAO facilities at Kitt Peak and the NRAO/VLA in Socorro.

The topics of the 2004 REU informal talks are listed below by site.

a. NSO/Tucson Talks

                                 Table 1. Informal Staff Talks at Tucson
June 7            Dr. Richard Green (National Optical              Introduction to Facilities and Activities at the
                  Astronomy Observatory                            National Optical Astronomy Observatory
June 14           Dr. Frank Hill (National Solar Observatory)      The Advanced Technology Solar Telescope
June 21           Dr. Rachel Howe (National Solar                  An Introduction to Helioseismology
June 28           Dr. Michael Merrill (National Optical            An Introduction to Infrared Astrophysics
                  Astronomy Observatory)
July 7            Dr. Steve Howell (Wisconsin-Indiana-Yale-        The Kepler Mission
                  NOAO (WIYN) Observatory)
July 7            Dr. Patricia Knezek (Wisconsin-Indiana-Yale-     Evolution in Astronomy: Galaxies Do It, So
                  NOAO (WIYN) Observatory)                         Why Shouldn’t We?
July 12           Dr. Chuck Claver (National Optical               The Large Synoptic Survey Telescope (LSST)
                  Astronomy Observatory)

4 NSO REU/RET Annual Program Report 2004                                                          Program Activities
b. NSO/Sacramento Peak

                          Table 2. Informal Staff Talks at Sacramento Peak
June 10        Dr. Thomas Rimmele (National Solar                  High-Resolution Solar Observations with
               Observatory)                                        Adaptive Optics
June 18        Dr. Han Uitenbroek (National Solar Observatory)     What Spectral Lines Tell Us About the Sun

July 2         Dr. Steve Keil (National Solar Observatory)         The Evolution of the National Solar
July 2         Dr. Scot Kleinman (Apache Point Observatory)        The Sloan Digital Sky Survey
July 9         Dr. Richard Altrock (Air Force Research             Long Term Variation of Solar Coronal
               Laboratory/Sac Peak)                                Temperature
July 14        Mr. Dave Dooling (National Solar Observatory)       The Sun at 20 Kilometers
July 16        Dr. Joel Mozer (Air Force Research Laboratory,      Space Weather
               Sac Peak)
July 20        Dr. K.S. Balasubramaniam (National Solar            High-Resolution Magnetography of
               Observatory)                                        Sunspots
July 23        Dr. K. Sankarasubramaniam (National Solar           Solar Polarimetry
August 4       Mr. Creighton Wilson (2004 NSO Summer RET,          The Neutral Line Method of Flare
               Lovelady High School, Lovelady, TX)                 Prediction from Magnetograms
August 4       Mr. Michael Sinclair (2004 NSO Summer RET,          Deconvolving the Point Spread Function of
               Kalamazoo Math & Science Center, Kalamazoo, MI)     SMEI (Solar Mass Ejection Imager) Imagery
August 6       Mr. Brian Harker-Lundberg (2004 NSO Graduate        Observational Signatures of Flux Emergence
               Summer Research Assistant, Utah State University)
August 6       Mr. Andrew Medlin (2004 NSO Graduate                Searching for Alfven Waves in the Solar
               Summer Research Assistant, New Mexico Institute     Atmosphere
               of Technology)
August 9       Ms. Leah Simon (2004 NSO Graduate Summer            Acoustic Events and K2V Bright Points (or,
               Research Assistant, Macalester College)             Wrestling with Large IBIS Data Sets)
August 9       Dr. Raymond Smartt (Emeritus, National Solar        Point Diffraction Interferometry
August 10      Dr. Alexei Pevtsov & Dr. Raymond Smartt             A Little Information about Australia
               (National Solar Observatory)

August 11      Ms. Cheryl-Annette Kincaid (2004 Air Force          Solar Radio Burst Locator
               Research Laboratory/Sac Peak Undergraduate
               Space Scholar, North Texas State University)
August 12      Dr. Alessandro Cacciani (University La Sapienza,    Recent Progress in Analysis of Doppler
               Rome, Italy)                                        Magnetic Images at Two Levels in the Solar
August 27      Dr. Steve Tomczyk (High Altitude Observatory)       Initial Measurements from the Coronal
                                                                   Multi-Channel Polarimeter
September 13   Ms. Maria Kazachenko (2004 Undergraduate            Solar Coronal Heating By Emerging Active
               Summer Research Assistant, St. Petersburg State     Regions
               University, St. Petersburg, Russia)

5 NSO REU/RET Annual Program Report 2004                                                       Program Activities
6 NSO REU/RET Annual Program Report 2004   APPENDIX III – Program Evaluation
                                               APPENDIX I
                  NSO Research Experiences for Teachers (RET) Program
During the summer of 2004, four teachers—two at NSO/Sacramento Peak and two at NSO/Tucson—
participated in the RET program, which is funded as a supplement to the NSO REU grant.

Mark A. Calhoun (Sabino High School, Tucson, Arizona), a physics teacher, devised classroom topical
lessons concerning geostationary satellites (geosats) as his primary summer project. He also established a
program that would enable high school students in the Tucson Unified School District to receive college
preparation credit for astronomy courses. Geosats, the source of most television broadcasts and a vital part
of telephone communication systems, play a key role in society today. One might say that, except for the sun
and moon, geosats are the most important celestial objects in daily life; yet few people know anything about
them. Mark Calhoun learned about the celestial mechanics of these devices, their construction, and how to
map and identify them. These activities were converted to lesson plans that will enable students to photo-
graph and likewise identify these approximately 12th-magnitude objects. Mark’s mentor was Bill Livingston.

Matthew Dawson (Brockton High School, Brockton, Massachusetts) worked with Rob Hubbard on issues
relating to weather and temperature prediction algorithms. The Advanced Technology Solar Telescope
(ATST) project needs to predict the outside air temperature several hours in advance at the observing site.
This is necessary so that the primary mirror temperature can be maintained at or near the air temperature in
spite of the mirror's long thermal time constant. Matt investigated a variety of methods currently in use for
this purpose, and tried them out using historical temperature data obtained at the three candidate ATST
sites. He found that by applying relatively simple algorithms he was able to predict the temperature two
hours in advance with a root mean square (RMS) accuracy of about one degree centigrade. Matt plans to
challenge his science students at home to try their own prediction schemes to make educated guesses of
future temperatures.

Michael Sinclair (Kalamazoo Science and Math Center, Kalamazoo, Michigan) performed a photometric
analysis of stellar data from the Air Force's Solar Mass Ejection Imager (SMEI) instrument aboard the
Coriolis spacecraft. His primary contribution was to develop a deconvolution technique to remove the
significant coma aberration that is inherent to raw SMEI imagery. Using this algorithm, Mr. Sinclair was able
to reconstruct sharp, photometrically accurate whole-sky maps that will be used in the future for various
stellar variability and related studies. Joel Mozer was Michael’s mentor.

Creighton Wilson (Lovelady High School, Lovelady, Texas) investigated magnetic properties of emerging
active regions, and how these properties can be used in forecasting flare activity of solar active regions. He
worked with Alexei Pevtsov on developing the IDL software to identify and track the areas of emerging active
regions. The software was applied to full-disk SoHO/MDI magnetograms to compute various magnetic
parameters including magnetic flux, area, polarity separation, active region's tilt, and length of magnetic
neutral line. These parameters and their time-evolution were compared with the flare activity of mature
active regions. The results suggest that this approach may be useful in forecasting the flare activity of active
regions, but further studies are necessary to quantify this possible relationship. Creighton plans to develop a
classroom activity based on this project and implement it in his high school.

The project reports by Mark Calhoun, Matt Dawson, Michael Sinclair and Creighton Wilson follow.

7 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
                                  National Science Foundation
                                Research Experience for Teachers
                                      Summer 2004 Report

  Imaging and Identification of Geostationary Satellites
                                             Mark A. Calhoun
                                    Sabino High School, Tucson, Arizona

                                     National Solar Observatory, Tucson
                                      Advisor: Dr. William Livingston

NSO Research Experience for Teachers (RET)
         The National Solar Observatory RET Program, sponsored by the National Science Foundation,
provides summer opportunities to pre-college teachers to participate in contemporary scientific astronomy
research at their NSO facility on the University of Arizona campus. Enthusiasm is high for the opportunity to
pursue professional development opportunities to practice what they teach, unencumbered by the rigors of
classroom teaching. Participants spend summer months working directly with science research mentors,
faculty, and post-docs followed by the development of instruments for incorporating outcomes and for
assessing the impacts of the RET on K-12 classroom activities. The successful migration of the astronomy
research experience to the teacher’s classroom is a major RET program goal for the NSO. In general, it is
hoped that the results include a long-term renewal in the joy of science and the continual transfer of that joy
and interest to the teacher’s students. Increased student interest and motivation in science will occur only if
the excitement of cutting edge research and the skills gained by the teacher can be transferred to the
classroom successfully.

         The NSO RET projects vary significantly in approach to meet the program’s goals, from structured
mentor-driven research experiences to more open teacher-driven ones. Primarily, the focus of the program
must be on the design of the entire participation experience. This planning involves the nature of the
research work, the teacher’s readiness and commitment, the mentor’s resources, the development of
classroom activities, and the professional environment. Ideally, the teacher’s project is motivating,
challenging, and “do-able” within the timeframe of the brief appointment. The NSO RET project taken on by
the teacher will be determined in part by the mentor’s research agenda, resources, and summer activities.
(Here, the “mentor” is assumed to be the lead scientist on the research project.) Additionally, team members
can discuss and plan with the teacher the migration of research results to the teacher’s classroom, identifying
target courses, and monitoring the formulation of research-related classroom activities for the next academic
year. The optimum RET experience provides challenging, motivating, and meaningful research work while
establishing a new fertile instructional resource for the teacher. Such activities may, for example, require
development of visual, graphical, and textual materials (conventional or e-based). When the teacher returns
to the classroom, and transfers the excitement of science and scientific discovery to students, the research
experience is successful in achieving programmatic goals, whatever the outcomes of the teacher’s specific
experimental research project efforts.

8 NSO REU/RET Annual Program Report 2004                                           APPENDIX I – RET Program
Imaging and Identification of Geostationary Satellites
1. Background on Geostationary Satellites
         Arguably the greatest impact that the Space Age has had on our everyday lives is the invention of the
communications satellite. In just a few decades, these devices have brought together the most far-flung
reaches of the globe in ways that not that long ago were barely imaginable. Today it is possible to talk
directly to anyone anywhere on and above the Earth, or communicate via the Internet with virtually any
computer system on the face of the planet—all with the help of communications satellites. While
communications satellites perform their missions in many types of orbits, from near-earth orbits like Global
Star to the highly-eccentric Molniya orbits used by the Russian Federation, one of the more important classes
of orbits for these satellites is the geostationary orbit. I'd like to examine the unique aspects of this class of
orbits (geostationary) which make it suitable for not only satellite communications, but for early warning and
weather observations, too.

         The satellites that orbit high above Earth's equator at a mean height of 35,786 kilometers (22,220
miles) have orbital speeds that match the turning of Earth far below, allowing them to appear as a stationary
target for television signal reception dishes and ground-based communications antennae. As a result, they are
called "geostationary" satellites, or geo-sats. It should also be clear that it is not possible to orbit a satellite
which is stationary over a single point on Earth which is not on the equator. This limitation is not serious
since most of the earth's surface is visible from geostationary orbit. In fact, a single geostationary satellite
can see 42 percent of the earth's surface and a constellation of three geostationary satellites can see the
earth’s entire surface between 81° S and 81° N.

         However, satellites in this orbit do experience a slight drift from the ideal position, which arises due
to anomalies in the Earth's gravitational field. The noncircular shape of the earth's equator causes these
satellites to be slowly drawn to one of two stable equilibrium points along the equator, resulting in an east-
west librations (drifting back and forth) about these points. The gravitational influences of the Sun and Moon
provides an out-of-plane force too, which gradually increases the orbital inclination towards that of the Moon
about the Earth. The geostationary satellites now tend to describe a figure-of-eight ground track; ground
controllers aim to restrict this to some predetermined box of error given that enough orbit-keeping fuel
remains. To counteract these perturbations, sufficient fuel is loaded into all geostationary satellites to
periodically correct any changes over the planned lifetime of the satellite. These periodic corrections are
known as station-keeping. North-south station-keeping corrects the slowly increasing inclination back to zero
and east-west station-keeping keeps the satellite at its assigned position within the geostationary belt. These
maneuvers are planned to maintain the geostationary satellite within a small distance of its ideal location.

        Due to the popularity of this orbit (geostationary slots over many regions are highly crowded and
greatly valued) some agencies de-orbit their satellites into a graveyard orbit some 500 to 1000 km above or
below their operational altitude. This would be the ideal situation but the final graveyard orbit is really
dictated by the remaining fuel on board. Failure to put these geo-sats in a graveyard orbit would put some of
the valuable active satellites at risk once station-keeping fuel was exhausted and the dead satellites would be
at the mercy of gravitational forces.

2. Research Introduction
          Unlike objects in low Earth, or highly elliptical orbits, geostationary satellites are visible throughout
every night of the year around a couple of weeks either side of each equinox. During the same period the
satellite tends to brighten over several days, twice a year, when the satellites orientation favors the 'beaming'

9 NSO REU/RET Annual Program Report 2004                                               APPENDIX I – RET Program
of the Sun in the direction of the observer. The source of the light detected is reflected Sunlight from the
satellite’s solar cell-covered wings that can reach more than 25 meters across. Typically the satellite’s
reflection will be in the mag. +11 to +14 range (or dimmer), but brightening by several magnitudes when the
geometry is favorable (around mag. +5 to +6 is not untypical). Surprisingly, given dark enough skies, it is
possible, armed with a stationary off-the-shelf camera to capture most of the satellites nestling in the
geostationary ring. Stars and gaseous nebulae leave linear streaks in the images, due to the rotation of Earth
during the exposure. The photograph below includes pointers at the Orion Nebula and the bright star Rigel,
as examples. These natural features can be identified with some basic astronomical knowledge and a star

        My NSO research project during the summer of 2004 was to successfully identify all detectable
geostationary satellites along the visible portion of the celestial equator from Kitt Peak observatory just
outside Tucson, Arizona. The method used to collect and analyze the data started in the form of
photographs, then using image processing software and some simple mathematics the detectable satellites
could theoretically then be identified using a published index of two line elements (TLE) for satellites. The
matter of identifying all of the detected satellites becomes a repeatable process of measuring the spacing
between the satellites through the image processing software and based on a scale derived from the
magnification of the camera lens. Once converted to degrees of longitude, the satellite coordinates and
spacing is then compared to the actual values as listed in the published TLE index or similar records
available on the internet.

3. Data and Analysis
         Minimal equipment is required to collect, analyze, and identify the data in this research project: a
plain old camera with Ektachrome 100G film, a simple compass, any brand of image processing software,
image scanner, and an updated catalog of satellite locations or internet connection. The steps involved in the
research were acquiring, processing, and analyzing the data, then finally identifying the geostationary

         The first step in the process of acquiring the data was to determine the best time and location for
photographing the geo-sats. As previously pointed out it is more rewarding to photograph these satellites
around the spring and fall equinoxes as well as around the new Moon phase so that the satellite will be more
visible. Data collection (photographs) over successive nights before and after this time will allow you to view
an increased brightening of the object. Of course this will also avoid the disappointment of poor data
collections whilst they may be in eclipse or washed out by the light of the Moon! In addition, it would be
best if ambient light were kept to a minimum. For this reason the location chosen was the Kitt Peak
Observatory about 50 miles outside Tucson, Arizona.

       The next step in capturing the images effectively was to determine the specific direction to aim the
camera to capture the geostationary satellites. By virtue of these satellites being stationary relative to the

10 NSO REU/RET Annual Program Report 2004                                           APPENDIX I – RET Program
Earth, they would have to travel along the celestial equator. The camera would need to be pointed towards
the azimuth and altitude corresponding to the visible celestial equator for the observer’s latitude. Except at
the equator, one must realize that an angle compensation for the latitude-dependent parallax must be
incorporated since the latitude (90o - ω), where ω is the local latitude that is directed at the celestial equator
that extends to infinity yet the geostationary orbit is at a distance of only 35,786 kilometers above the surface
of the Earth. A correction was necessary in order to compensate for this parallax effect that exists when using
celestial coordinates that actually extend to infinity. The trigonometry involved in determining this corrective
angle for the location at the Kitt Peak Observatory is as follows:

Equatorial plane                                            32 degree latitude plane (zenith)

True direction towards geo-sats                             32 degree down from Kitt Peak zenith

Geo-sat location

                                                                    ω     A


        Kitt Peak is located on the circumference of the circle (Earth) at the intersection point of the three
lines and is approximately at +32 o latitude, which extended along the zenith. If the camera angle were
directed 32 o below this intersection point, then the line of site would be parallel to the celestial equator and
not towards the geo-sats that are to be photographed. As a result a corrective angle must be determined to
ensure they are photographically captured. As mentioned before, the line of sight (AD) along (90 o - ω) would
be parallel to the celestial equator line BC, but the camera should be angled along the line AC to capture the
geo-stats. With the law of cosines this corrective angle at the vertex “C” of triangle “ABC” can be

        Radius of Earth          Geo-sat Orbital Distance           Kitt Peak Latitude
        AB = 6372 km             BC = 42164 km                      Angle ω = 31.95 o

         First side AC must be calculated knowing the angle opposite that side and the other two sides of
triangle “ABC”. Thereafter, the three sides can be used to calculate angle “C” again using the law of cosines.
The corrective value of 5.24 degrees was found for the Kitt Peak location and thus the camera needed to be
directed at [(90 o - ω) - 5.2 o ] or 52.8 o above the horizon.

11 NSO REU/RET Annual Program Report 2004                                             APPENDIX I – RET Program
           Finally the photos can be acquired. This process requires the use of long-exposure photography in
order to collect enough light reflected from the faint satellites that only have apparent brightness factors from
+10 to +18. As it is these brightness factors are about 100 times fainter than the faintest stars viewed from
Earth. Using an exposure duration of six to eight hours, the camera captured the geo-sats in the field of view
whilst the stars drift in and out courtesy of the Earth's rotation. It should be noted here that a tripod or other
means of maintaining a stationary camera was necessary for this exposure time. Else the images will become
blurred as the position of the satellites would appear to change relative to the shifting camera.

         Now that the images of the geo-sats have been acquired and after they are developed they need to be
processed for data reduction that would be used to finally identify them individually by name. This process
begins by scanning the photos into any image processing software for optimizing and then analysis. These
specific steps will vary slightly with the brand of scanner and image processing (IP) software used. Once
these photos were imported to the IP software they were magnified, “sharpened”, and the contrast was
adjusted to optimize image quality before image analysis could be effectively completed. At this point the
vertical and horizontal measurement scales were imposed on the image in order to identify and record
relative positions and spacing of each of the detected satellites.

         The relative positioning data of the satellites from the IP software scale values was then
mathematically manipulated to determine the experimental latitude and longitude positions of the geo-sats.
This was achieved by first calculating the field of view angle for the camera used. This angle was used to
calculating the degrees per centimeter value for the final optimized image that spanned 437 cm. The camera
used here had an 80 mm focal length and a 55.5 mm imaging width. This was used to calculate a field of
view angle of 38.26 o for the photographs taken of the geo-sats in orbit.

                                     υ = Tan-1 [(55.5/2)/80]
                                     υ = 19.13 o

                                     Total camera field angle is 2υ
                                     or 38.26 o

                                              The degrees of longitude per centimeter of IP
                                              image would then be:

                                              38.26 o / 437 cm = .0876 deg./cm


Given that the longitude at Kitt Peak is – 111.6 o and the field of view of the camera is 38.26 o, it was
established that the far eastern point of the IP geo-sat images would be at 92.48 o of longitude and the most
western point would be at 130.74 o of longitude. This being established, the degrees of longitude for each
satellite position previously recorded in centimeters was used to determine their corresponding longitudinal
position. Once the longitudinal positions of all detectable objects along the equatorial plan are logged, then
they can be used in conjunction with a published index of current satellites and their TLE to identify the
satellites by name. I used these internet locations for referencing the location and name of each satellite
photographically       detected     in      this    research,          and .

12 NSO REU/RET Annual Program Report 2004                                             APPENDIX I – RET Program
        Results:                                   4. Results

                                     Telstar 6               I had determined that there were a
                                                   total of twenty-six detectable objects along
                                   Galaxy 8-i      just over thirty-eight degrees of the
                                                   geostationary orbit above Kitt Peak. Based
                                                   on the calculated longitudinal positions I was
                                                   able to identify all of them by name with an
                                     Telstar 5     acceptable measure of uncertainty. My main
                                                   objective was to develop a unique context
                                                   that could effectively be used to instruct
                                   Galaxy 4R       concepts related to orbits through Kepler’s
                                                   and Newton’s Laws. Through this research I
                                                   was to determine if the collection,
                                                   mathematical reduction, and analysis of the
                             Direct TV-2           data in the form of real images taken of the
                             AMSC-1                geostationary satellites above Tucson,
                         }   GE 4
                             Direct TV-1R
                                                   Arizona was possible within the realm of
                             Direct TV-4S          students’ interest and abilities. Having gone
                                                   through the entire process four times and
                                         GE 1      working out any foreseeable pitfalls, I
                                                   believe that this is not only a unique
                                     GOES 11       approach to teaching orbital mechanics, but
                                                   one that students will be able to relate as
                                                   these satellites are a large part of their world
                                       MSAT        in the form of television, cell phones, and the
                                                   fact that these geostationary satellites are
                                                   visually accessible in the night sky outside
                                       Anik F1     their bedroom window. Students will be able
                             Direct TV-1           to experience the skills involved in scientific
                         }   Echostar 8            process that not only require hands-on
                             Echostar 5            learning, but minds-on as well. My approach
                                        Anik E2    in getting students interested in a particular
                                   Solidaridad 2   concept in physics or astronomy can be best
                                                   summed up by this quote by Antoine de
                                                   Saint-Exupery, “If you want to build a ship,
                                                   don’t herd people together to collect wood
                                                   and don’t assign them tasks and work, but
                                                   rather teach them to long for the endless
                                       Satmex 5    immensity of the sea.” In today’s words,
                                                   teaching a love of learning will help my
                               Anik E1
                                                   students motivate themselves. I truly believe
                                                   that it is the essence of teaching.
                          }    Tempo 2
                               Echostar 7
                               Echostar 6
                               Echostar 4
                                     Galaxy 10R

13 NSO REU/RET Annual Program Report 2004                                      APPENDIX I – RET Program
5. Developed Educational Activities
         As a result of my research project, geostationary satellite identification, I have developed several
related lessons that can be used to study satellites and the mechanics involved in their orbits. An emphasis is
given to geosynchronous satellites, but general orbital mechanics is also present in several of the activities.
The order that they are listed is only a suggested completion sequence. It is hoped that these activities will
provide a unique and interesting means for studying orbital mechanics, while also coming to understand the
technologies involved in providing television, cell phone, and other amenities that we have taken for granted.

         Satellite Technology & Purpose: This activity involves a video (Eye in the Sky) about the history,
nature, type, and purpose of satellites. A video question sheet is included to aid in the student comprehension
of the video. Also included is a teacher guide and grading rubric to facilitate a student-centered research and
discussion activity.

          Satellite Orbital Activity: Students will explore various Internet websites that simulate the orbit of
a planet around the sun and provide current orbital data about earth-orbiting satellites. Using a simulation
they will investigate characteristic of satellite orbits. They will then use this information along with the
centripetal force equation and Newton’s law of Universal gravitation to calculate other features of the
satellites’ orbit (velocity, altitude, period of revolution).

        Sinusoidal Orbital Mapped Paths: This is a hands-on construction activity for students to discover
why a sine curve orbital ground plot is created when an orbit is portrayed on a flat 2-D map.

        Geostationary & Geosynchronous Satellites: A colorful teacher slide presentation on
geostationary vs. geosynchronous orbits and their corresponding ground tracks. A follow-up student activity
on geostationary and geosynchronous satellite ground tracks and location visibility of satellites in these

        Slider Map Activity: In this activity students become familiar with the positioning of a “Slider”
Map and marking the longitude of day and the longitude of night (i.e. identifying the portion of the ground
track that is in daylight. Several of the sample orbits represent special circumstances which require students
to reposition the “Slider” Map in order to view all their targeting opportunities.

        Satellite Ground Tracks: This lesson included a teacher colorful slide presentation of orbital paths
for satellites and their ground tracks called “footprints”. Using satellite ground tracks students are to
mathematically determine the orbital, horizon, and ground track distances based on satellites orbit.

         Identifying Geostationary Satellites: This is a student research-type activity where students use
geo-sat images, image processing software, simple trig-based mathematics, and the scientific process to
identify by name the satellites in the given image. Students will require an overview of geostationary
satellites which should be provided by the “Observing Geostationary Satellites” and other prerequisite
activities listed herein. These prerequisites will act as a means of providing interest, motivation, and
additional background for this student-centered group activity pertaining to geostationary satellites.

         Satellite Tracking via the Internet: Students will use Newton's Universal Law of Gravity and the
equation for centripetal force for calculating the orbital velocities, location, and periods for the ISS. They
will view and comprehend real-time flat maps showing satellite orbital paths, location, and data on the
Internet. They will then use this Internet resource to confirm all results.

14 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
        Changing Satellite Orbits: A colorful teachers’ slide presentation on the logistics and mathematics
involved in changing the orbit of satellites. A follow-up activity pertaining to the vector mathematics
involved in changing the orbits of satellites.

6. References

Hernandez, C., Jehn, R. Classification of Geosynchronous Objects, Issue 6: Jan. 2004.
NASA, ISS EarthKAM website:
Canadian Space Agency web page:
Geostationary Satellites web site:
Douglas Isbell, “Satellites Spanning the Sky”, NOAO News, Feb. 28, 2001

7. Acknowledgements
         I’d like to graciously thank my RET advisor Dr. Bill Livingston for his time, guidance, and
motivation while working under him for the summer. His patience and support was of great value throughout
the process of my research project. The RET program wouldn’t be possible without such scientist as Dr.
Livingston who put a part of their lives aside for ten weeks to assist a teacher who they never knew or will
probable ever see again.
         I’d also like to thank Priscilla Piano for making me feel as though I was a part of the NSO team by
providing a great workspace, mailbox, ID card, e-mail account, invitations to NOAO and NSO lectures, and
entrusting me to chaperone the REU students to Sacramento Peak. Her logistical contributions really helped
make my experiences here at NSO a productive, rewarding, and fun endeavor. I really enjoyed the trips to the
VLA, Sunspot, and White Sands.
         I’d like to also thank the Sunspot staff for a great stay at Sacramento Peak observatory. The July 4
picnic and activities was such a great time.
         Finally, I’d like to express my gratitude to the National Science Foundation for seeing a value in
making available opportunities for teachers to experience scientific research first hand. I truly feel
reenergized and hyper-motivated to return to my classroom with a renewed interest in getting students
excited about science and the scientific process.

15 NSO REU/RET Annual Program Report 2004                                          APPENDIX I – RET Program
                                   National Science Foundation
                                 Research Experience for Teachers
                                       Summer 2004 Report

                  Nowcasting Air Temperature at the
              Advanced Technology Solar Telescope (ATST)
                                              Matthew Dawson
                               Brockton High School, Brockton, Massachusetts

                                      National Solar Observatory, Tucson
                                          Advisor: Robert Hubbard

             “Those who have knowledge, don't predict. Those who predict, don't have knowledge.”
                                        Lao Tzu, Father of Taoism

              “No weather predictions or forecasts are available for the summit of the mountain.”

                        “The weather at Haleakala is unpredictable and ever-changing.”

“Since the weather at both the summit and Kipahulu areas of the park is unpredictable and can change quickly, it
                                 is important to be prepared for any condition.”
                                           The National Park Service

  “Sunshine is delicious, rain is refreshing, wind braces up, snow is exhilarating; there is no such thing as bad
                                   weather just different kinds of good weather.”
                                           John Ruskin, Victorian thinker

                             “Climate is what you expect, weather is what you get.”
                                        Ancient Meteorological Proverb

               “Climate tells you what clothes to buy, but weather tells you what clothes to wear."
                                          Some Middle School Student


The primary mirror of the Advanced Technology Solar Telescope (ATST) has a diameter of 4 meters. The
purpose of this mirror is to collect light from the Sun over extended periods of time. Like anything
constantly bathed in sunlight (think of your steering wheel on a summer day in Tucson), the ATST’s primary
mirror will get very, very hot. The mirror will absorb so much of the Sun’s heat that it will cause the air
immediately around it to heat up as well. When the temperature of air increases, its density decreases.
Therefore, the warm air near the mirror will be less dense than the cool air above it. Fluids with lower

16 NSO REU/RET Annual Program Report 2004                                               APPENDIX I – RET Program
densities prefer to “float” on fluids with higher densities, so the warm air will rise and be replaced by the
cool air that sinks into its place. In this cyclic process, known as convection, warm air and cold air are
always on the move as they constantly trade places. This motion is an obstacle for telescope observations.
Light travels at different speeds through air with different densities, so circulating pockets of cold and warm
air cause light waves to refract in complex patterns. This causes the image recorded by the telescope to be
shifted, distorted, blurred, and reduced in both contrast and resolution. These visual anomalies caused by the
telescope’s mirror are called “mirror seeing.” To reduce mirror seeing to an acceptable level, the convection
of warm air rising from the primary mirror must be minimized by maintaining a small temperature difference
(∆T) between the mirror and the air around it.

        How small does the temperature difference (∆T) between the mirror and air have to be? According
to an analysis by Dalrymple (2002):

        • “Still air over a heated mirror results in bad seeing.
        • Specifically, the mirror DT must be kept below about 0.5 K in still air to stay within the budgeted
          error of 0.15 arcsec mirror seeing.
        • Adding even weak flushing or wind over the mirror reduces the mirror seeing dramatically. At a
          flushing velocity of 1 m/s, the allowable DT increases to 1.7 K and at 2 m/s, DT can be several
          degrees.”(Dalrymple, 2002)

        In a later report, Dalrymple notes that “Estimates indicate that a surface–air temperature difference
∆T ≤ 1 K or so will adequately control mirror seeing.” (Dalrymple, 2003) To accomplish this, the mirror’s
surface will be cooled with a combination of air jets and liquid coolant systems.

         Due to the mirror’s size, up to one hour may be required to adjust its temperature by 1 K (or 1°C).
This time lag requires the ambient air temperature to be anticipated a few hours in advance. Therefore, a
method of accurately forecasting the ambient air temperature around the ATST must be developed. The goal
of this project is to develop a technique to “nowcast” the air temperature at the ATST within °C, within a 3
hour horizon.

         Nowcasts of the ambient air temperature around the ATST serve a few functions. Air temperature
nowcasts can be used to preset the cooling system to minimize the temperature difference between the
primary mirror and the air surrounding it, thereby reducing mirror seeing and improving observation quality.
Air temperature nowcasts might help to predict when thermal conditions are favorable (or unfavorable) for
certain types of observations, thus serving as a scheduling or prioritization tool. Nowcasts of air temperature
can also be used as a guide for optimization the set-up of the ATST. In general, the ability to accurately
nowcast the ambient air temperature around the ATST will enhance the telescope’s overall efficiency and

Categories of Forecasts
         In the field of meteorology, there are three general categories of forecasts: long-term forecasts, short-
term forecasts, and “nowcasts.” The main difference between these forecast categories is the “horizon,” the
period of time covered by a forecast’s predictions. Long-term forecasts have horizons of 24 hours or greater
and are typically generated by complex physical models developed by meteorological institutes. Long-term
forecasts usually take on a global perspective with relatively low spatial (tens of km) and temporal (6 to 12
hours) resolution. Short-term forecasts have horizons between 1 to 12 hours. Like long-term forecasts,
short-term forecasts may be generated by physical models, but their resolution may still be too low for some
purposes. Nowcasts, the category of forecast dealt with in this project, are made when the desired horizon is
less than 6 hours, oftentimes only 1 or 2 hours. (Murtag, Sarazin)

17 NSO REU/RET Annual Program Report 2004                                              APPENDIX I – RET Program
Approaches to Forecasting
        There are three basic approaches to producing meteorological forecasts. Forecasts may be generated
through the use of physical models, statistical analyses, or a combination of both. Each approach to
forecasting has unique advantages and disadvantages.

         Physical models simulate the actual conditions and processes in Earth’s atmosphere. Like Earth’s
atmosphere, physical models can be very complex; they often require large quantities of input data, extended
computing time, and a specialized knowledge of meteorology. Because the equations used in physical
models are solved numerically by computer, these models are commonly known as Numerical Weather
Predictions (NWPs). NWPs are mostly used to generate long-term forecasts that cover large regions, but at
relatively low resolutions. Sophisticated NWPs are used by the National Oceanic and Atmospheric
Administration (NOAA), National Weather Service (NWS), National Centers for Environmental Prediction
(NCEP), National Center for Atmospheric Research (NCAR), European Centre for Medium-Range Weather
Forecasts (ECMWF), and other major institutes of meteorology to generate both long-term and short-term
forecasts. Because of the complexity and expertise involved, production of a physical model to forecast
ambient air temperature at the ATST is beyond the scope of this project. However, the Hawaii
Weather/Climate Modeling Ohana (HWCMO) uses two different NWP systems to forecast a number of
meteorological variables, like temperature, precipitation, wind, and cloud cover, at high resolution over the
major Hawaiian Islands, including Maui, home of Haleakala. One system, a Regional Spectral Model
(RSM), generates forecasts over a 48 hour horizon, at 6 hour intervals, with a 3 km spatial resolution over all
of Maui. The other system, called Mesoscale Model 5 (MM5), generates forecasts over a 48 hour horizon, at
1 hour intervals, with a 3 km spatial resolution over all of Maui, and with a 1 km resolution over the summit
of Haleakala. Although these forecasts are considered experimental, the motivation for the 1 km resolution
forecasts at Haleakala is to assist astronomical observations made at the observatories located on its summit.
(HWCMO website) If the ATST is built at Haleakala, MM5 forecasts would certainly be a useful aid in
telescope operations. If the ATST is built elsewhere, the NSO might consider using MM5 itself, since the
model can be adapted for use anywhere, as long as sufficient data is available. The MM5 software package
was developed by UCAR and Pennsylvania State University, and it is available for free download at the
MM5 website. Using MM5 to produce weather forecasts requires meteorological knowledge, numerical
modeling experience, and good computing facilities, so interdisciplinary collaboration with a local
meteorological institute, computing center, or university may be helpful (or even necessary). (MM5 website)

        ECMWF temperature forecasts are considered not precise enough for astronomical purposes;
standard error is “well higher than 2 degrees,” horizontal resolution is 60 kilometers, and although forecasts
have a horizon of 10 days, they are only updated every 12 hours. A more focused Limited Area Model
(LAM, also Local Area Model) could have a horizontal resolution of 30-15 km and might achieve standard
errors below 2°C, possibly as low as 1°C for a “stable region” like La Palma. The limitations are said to be
“on the computing side.” Several European countries are working on a High Resolution Limited Area Model
(HIRLAM). (Buffa & Porceddu 1997)

         Forecasts can also be made by conducting statistical analyses of historical time series data. Patterns
in past meteorological variables are used to predict future values of the variables. Unlike the physical
models of NWPs, little or no knowledge of the actual physical processes that govern the weather is used in
the statistical approach to forecasting. Statistical forecasts are advantageous because they perform well for
short- and intermediate-term forecasting, can be customized to a specific site, can easily be updated as new
data is collected, require much less input data than NWPs, use data that is easy to obtain, and need much less
computing time than NWPs. Although more overall input data improves the accuracy of statistical forecasts,
some statistical forecasts work best with few input variables, sometimes as few as one. The main
disadvantage of statistical forecasts is that they are not very good at predicting extreme changes in weather
conditions. Statistical forecasts have been used, with varying degrees of success, to predict a diverse array of

18 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
variables, including electrical load, wind power generation, avalanches, air temperature, stock prices, sunspot
numbers, and even astronomical seeing. To date, statistical forecasts have not been used in the context
presented here: predicting daytime air temperatures on very short horizons with an accuracy better than 1°C.

         To generate accurate forecasts that have a high resolution and the ability to anticipate extreme
weather changes, a combination of statistical analyses and physical models may be desirable. One possibility
is to use statistical methods to relate broad forecasts generated by NWPs to local conditions. For example,
NWS temperature forecasts for the Ranger Station at Haleakala National Park might be statistically “filtered”
to generate forecasts for Haleakala’s summit. Or perhaps statistical methods could be used to predict small-
scale variations from or errors in forecasts generated by physical modeling. Another possibility is to rely on
statistical forecasts for optimizing general telescope operations while using physical models to anticipate
extreme events that may provide opportunities for telescope downtime and maintenance.

Assessing Forecasts
        The errors in a forecast can be measured in a few different ways. The two error measures most
widely used by forecasters are Root Mean Square Error (RMSE) and Mean Absolute Error (MAE). The
RMSE is given by

                                      Σ(X -X )
                                             p   o


where Xp is the predicted value, Xo is the observed value, and N is the total number of predictions. The
MAE is just the mean of the absolute values of all of the errors, where (Xp-Xo) is the error. Forecast errors
can also be displayed as a histogram showing the frequency distribution of error. Error frequencies are
typically normalized to 1 by dividing them by the total number of predictions or measurements. Another
way to measure forecast error is by determining the correlation function, R, and R2 values that relate the
predictions to the observations. (Giebel 2003) In this project, RMSE, MAE, and histograms are used to
show the errors in air temperature nowcasts.

         The quality of a weather forecast can be assessed by comparing its errors to the errors produced by a
simple “benchmark” forecast. Forecasts based on “persistence” and “climate” are the simplest and most
traditional benchmark forecasts. A persistence forecast is made by assuming that conditions at time t will be
the same as the conditions at time t-dt. For example, to forecast the temperature at noon on June 21, 2005,
one may assume it will be the same as the temperature at noon on June 20, 2005. The persistence method of
forecasting is also called the “carbon copy” method or the “naïve predictor.” In many parts of the globe,
time scales in the atmosphere are in the order of days, so significant changes in weather usually take place
over a few days rather than a few hours. As a result, persistence forecasts work well for short horizons,
making them an appropriate benchmark for ATST ambient air temperature forecasts. (Giebel 2003)

         A climate based forecast is made by assuming that conditions at time t within a given year (or other
period or season) will be the same as the conditions at time t averaged over all previous years (or as many
years as possible). For example, to forecast the temperature at noon on June 21, 2005, one may assume it
will be equal to the average of all temperatures recorded at noon on June 21, through as many previous years
as possible. Climate forecasts seem easy to make, but a thorough model of a given area’s climate requires
many years of historical data, which may not be available. Also, climate forecasts are hampered by both
random climate fluctuations and general climate trends (like global warming). Therefore, while climate
forecasts are good at capturing broad seasonal trends, they often miss variations from those trends that might
otherwise be captured by persistence forecasts. Hence the saying, “climate is what you expect, weather is
what you get.”

19 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
Literature Review
        Murtagh et al. (1992) describe a series of meteorological forecast experiments carried out for Cerro
Paranal, Chile, during the austral summer of 1989-1990. Murtagh et al. attempted to predict the nighttime air
temperature 24 hours in advance and astronomical seeing a few hours in advance. The forecast methods
used were purely statistical; a 1-nearest neighbor method and a “multilayer perceptron” (neural network)
were used for both temperature and seeing forecasts. The 1-nearest neighbor method “performed well” for
both temperature forecasts and seeing forecasts. On the other hand, the neural network worked “acceptably
well” for seeing forecasts but performed “poorly” for temperature forecasts.

         The candidates and references in the 1-nearest neighbor temperature forecasts were 5-tuples,
consisting of {Tt-24, Tt-48, Tt-72, Pt-24, Pt-48}, where Tt-24 is the temperature at time t-24 hours and Pt-24 is the
pressure at time t-24 hours. “The sequence of temperatures is used in order to capture contemporary trend.
Pressure is used for its information on quickly changing weather conditions such as incoming weather fronts.
Longer sequences of temperatures, alone, were also assessed, but not further reported on.” 5260 sequential
tuples were used in all, with the first 4000 tuples serving as candidates and the last 1260 tuples serving as the
references for an equal number of predictions. All pressure and temperature values were range-normalized
to the interval [0,1], and nearest neighbors were found using Euclidean distance. No mention was made of
taking the first difference of the temperature or pressure time series. 1-nearest neighbor temperature
forecasts had errors < 0.5°C 62.3% of the time, < 1.0°C 74.3% of the time, and < 2°C 85.1% of the time. In
comparison, forecasts made using persistence from 24 hours before the forecast time had errors < 0.5°C
20.5% of the time. It was also found that “using temperature values alone and not pressure had indicated
improved results for 2° prediction, but poorer results for 0.5° prediction.” 1-nearest neighbor temperature
forecasts were also evaluated with respect to “predictability of large temperature changes (defined as > 3°
difference).” “The hit rate (a big change is correctly predicted, given that a large change actually takes
place) was 67.6%. The false alarm rate (a large change predicted, but does not ensue in reality) was 6.5%.”

        Alternative forecasting methods were suggested, including “nonlinear time series analysis,
clusterwise regression, and parametric, Bayesian approaches.” The use of a larger candidate set was also
proposed. (Murtagh et al. 1992, Murtagh & Sarazin 1995)

         Aussem et al. (1994) is an extension of the work presented by Murtagh et al. (1992). Aussem et al.
states “the use of nonlinear models is strongly motivated by the work carried out by Murtagh et al. (1992)
which shows that a 1-nearest neighbor method outperforms the fitting of linear autoregressive models.”
First, Aussem et al. forecasted the air temperature at sunset, 24 hours in advance, at Cerro La Silla and Cerro
Paranal, Chile, sites of the European South Observatory (ESO). The table below summarizes the errors in
three different forecasting methods used for La Silla. (Aussem et al. 1994)

Figure 1. Forecast errors for La Silla
 method                                error < 2.0°C             error < 1.0°C               error < 0.5°C
 neural network                           72.5%                     42.5%                       20.0%
 1-nearest neighbor                       59.8%                     37.8%                       18.3%
 persistence from 24 hours                56.7%                     31.3%                       18.8%
                                                                                              (Aussem et al. 1994)

The neural network forecasts performed better, with Normalized Mean Square Errors (NMSE) of 0.70°C for
La Silla and 1.20°C for Paranal, while the 1-nearest neighbor forecasts had a NMSE of 1.69°C at La Silla
(NMSE was not reported for Paranal). (Aussem et al. 1994)

        Aussem et al. also forecasted the air temperature at each hour of the day, 24 hours in advance. The
explored forecast methods include persistence from 24 hours before the forecast time, a multidimensional

20 NSO REU/RET Annual Program Report 2004                                                APPENDIX I – RET Program
autoregressive model, two neural network variations, and three 1-nearest neighbor variations. As in Murtagh
et al (1992), the candidates and references for the nearest neighbor forecasts were 5-tuples consisting of {Tt-
24, Tt-48, Tt-72, Pt-24, Pt-48}, and closeness was measured using “Euclidean distance between range normalized
variables.” The three 1-nearest neighbor forecasts were made as follows:

        • 1-nearest neighbor (1)- 1260 predictions were made using 1260 reference tuples compared to
          chronologically prior tuples.
        • 1-nearest neighbor (2)- 1260 predictions were made using 1260 reference tuples compared to a set
          of 4000 prototype tuples.
        • 1-nearest neighbor (3)- like 1-nearest neighbor (1), except pressure variables are omitted.

Errors in the forecasts of air temperature at each hour of the day, 24 hours in advance are displayed in the
table below. (Aussem et al. 1994)

Figure 2. Forecast errors
method                                   error < 2.0°C             error < 1.0°C            error < 0.5°C
1-nearest neighbor (1)                      83.6%                     72.7%                    62.3%
1-nearest neighbor (2)                      66.7%                     48.6%                    32.3%
1-nearest neighbor (3)                      68.1%                     50.2%                    37.1%
neural network (1)                          71.1%                     45.1%                    22.0%
neural network (2)                          69.7%                     41.8%                    21.2%
persistence from 24 hours                    63.2%                    36.3%                    18.8%
multidimensional autoregressive             63.7%                     35.8%                    16.3%
                                                                                          (Aussem et al. 1994)

        Of all the methods used to forecast air temperature at each hour of the day, 24 hours in advance, the
1-nearest neighbor (1) approach is clearly the best. Aussem et al. suggest that the 1-nearest neighbor (1)
forecasts perform better than the 1-nearest neighbor (2) forecasts because “more historical information is
being used.” They also suggest that the 1-nearest neighbor (1) forecasts perform better than the 1-nearest
neighbor (3) forecasts because of the information contained in the pressure variables. Interestingly,
persistence seems to work as well or better than the multidimensional autoregressive approach. These results
must be taken with caution, as they only apply to Paranal and La Silla, the sites investigated in this study.
Furthermore, while the nearest neighbor method produced the best forecasts of air temperature at each hour
of the day, a neural network produced the best forecasts of the air temperature at sunset, leading the authors
to declare that “neither approach is uniquely preferable.” Therefore, the best conclusions that can be drawn
from this study are:

        • Short-term air temperature forecasts for the purpose of aiding astronomical observation are
        • The best forecast method depends on the specific characteristics of the site and the type of forecast
          desired. (Aussem et al. 1994)

         Buffa and Porceddu (1997) carried out “a feasibility study of temperature forecast in support to
active air conditioning of a telescope’s dome as part of a research activity which is done within the
meteorological support to the Italian national telescope Galileo (TNG), located at the Observatorio Roque de
Los Muchachos (ORM), in the Canary Islands. A neural network modeling is presented as compared to
classic linear filter algorithm.” (Buffa & Porceddu 1997)

        Short-range (1-12 hour) temperature forecasts were made using linear (L) and non-linear (NL)
autoregressive (AR) models. Non-linear autoregressive (NLAR) models were carried out using neural

21 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
networks. In addition to purely autoregressive models, autoregressive moving average (ARMA) models and
autoregressive moving average models with exogenous inputs (ARMAX) were also explored. The use of
non-linear models was motivated by the idea that a “non-linear approach allows modeling of complex
dynamics in climatic variables.” (Buffa & Porceddu 1997)

         A time series of temperature and pressure data was collected from the automatic weather station at
the Carlsberg Meridian Telescope (CMT), part of the Observatorio Del Roque de los Muchachos (ORM) on
La Palma island. The time series “corresponds to the period April-May 1993 for training, May-June 1993 for
validations.” Data were standardized to µ=0 and σ=1. Pressure values constituted the “exogenous inputs” in
the ARMAX models and were thus not used in any other model. The inclusion of pressure was suggested by
“an analysis of the cross-correlation structure between temperature and all the other meteorological
variables” collected by the CMT weather station. However, the results show that the use of pressure does not
improve forecasts and that the inclusion of other meteorological variables actually worsens forecasts. The
results, which are summarized in the table below, also indicate that “on site output prediction of a neural
network are competitive with respect to a linear prediction approach.” Buffa and Porceddu suggest exploring
additional neural network models that use the locally observed time series of temperature data and ECMWF
forecasts as inputs. (Buffa & Porceddu 1997)

Figure 3. Forecast errors for the CMT, at the ORM, La Palma
model                               # samples                horizon (hours)                RMSE (°C)
LAR                                    1320                         1                            0.7
LAR                                     220                         6                            2.2
LARMAX                                  220                         6                            3.1
NLAR                                   1320                         1                            0.6
NLAR                                    220                          6                           1.3
NLARMA                                  220                          6                           1.4
NLARMAX                                 220                         6                            1.4
NLARMAX(p)                             220                          6                            1.5
NLAR                                    110                         12                           2.1
ECMWF                                    ---               updated every 12 hrs        “ well higher than 2° ”
                                                                                     (Buffa & Porceddu 1997)

        Erasmus (2004) outlines a meteorological forecast program for the Southern African Large
Telescope (SALT). The SALT is to be located at the Southern African Astronomical Observatory (SAAO)
near Sutherland, in the Northern Cape (a four hour drive from Cape Town). The SALT forecasts combine a
physical model (or NWP) with some statistical analysis. According to the SALT weather website, “weather
forecasts for SALT are made using numerical meteorological forecast (Eta) model output data from the
South African Weather Service (SAWS) and meteorological observations made in situ at Sutherland
Observatory. These data are ingested by a forecast programme containing specially developed predictive
algorithms and a forecast of meteorological parameters relevant to SALT operations is made.” (SALT
weather website) These meteorological parameters include temperature, pressure, relative humidity, wind
speed, and wind direction. (Erasmus 2004)

        Temperature forecasts made using persistence from 24 hours before the forecast time had RMSE
between 4-6°C and MAE between 3-5°C. Temperature forecasts based on SAWS Eta output had RMSE
between 2-2.5°C during nighttime hours and 2.5-5.5°C during daytime hours. The MAE was between 1.5-
2°C during nighttime hours and between 2-5°C during daytime hours. It was expected that the use of a
Kalman filter would improve the temperature forecasts based on SAWS Eta output, particularly during the
daytime. The use of a Kalman filter reduced errors in the short-range forecasts (up to 9 hours), but actually
increased errors in the medium to long-range forecasts, where the errors sometimes exceeded those of the

22 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
persistence forecasts. Thus, the Kalman filter approach was determined to have “little value” for this
particular model and site. Further research will examine whether a Kalman filter “can be trained or used
differently in order to accomplish the desired result.” Other possibilities include using “additional Eta model
data” and using “the observed data in the days prior to the forecast period in some manner other than the
Kalman filter approach.” (Erasmus 2004)

        Raible et al. (2001) describes the use of a linear empirical forecast model to produce 24 hour ahead
forecasts of maximum and minimum daily temperatures at Hamburg, Germany. Linear empirical forecast
models are described as playing “the minor role of supporting the highly advanced numerical weather
prediction (NWP) systems.” The linear empirical forecast model used here is a “variant of the generalized
equivalent Markov system.” (Raible et al. 2001)

        Forecasts of maximum daily temperatures had an RMSE of 2.95°C. Forecasts of minimum daily
temperatures had an RMSE of 2.76°C. In comparison, the German Weather Service Europa Model NWP
forecasts had an RMSE of 3.34°C for maximum daily temperatures and 2.69°C for minimum daily
temperatures. When the linear empirical forecast model was combined with the NWP forecasts, the RMSE
was 2.74°C for maximum daily temperatures and 2.38°C for minimum daily temperatures. These results
support the practice of combining physical models (NWPs) with statistically based forecasts when possible.
(Raible et al. 2001)

The Sites
        Air temperature nowcasts at three locations: Haleakala, Tucson, and La Palma Island. Haleakala and
La Palma Island were obviously chosen because of their status as candidate sites for the ATST. Tucson was
chosen for two reasons: its stable weather might provide a good benchmark for assessing forecasting
methods, and good quality observations and forecasts are available for download at the NWS website.

         Meteorological data for the summit of Haleakala, elevation 10,023 feet, were collected by the
Advanced Electro-Optical System (AEOS) Haleakala Atmospheric Characterization (AHAC) group. The
data were downloaded from the AEOS AHAC ftp site by Frank Hill and then shared for use in this project.
The Haleakala data include air temperature, air pressure, wind speed, wind direction, precipitation, and
pyranometer readings. The Haleakala data covered the period from 1994 to 2002, and measurements were
recorded every ten minutes. To make the data sets more manageable, only measurements recorded on the
hour were used. Also, due to the amount of bad or missing data, lots of “cleaning” was necessary, so
forecasts were produced only for the year 2000. Time series data for the year 2000 were divided into two
continuous sub-series with no missing data, one sub-series running from 01:00:00 on 01/07/00 to 14:00:00
on 05/27/00 (~142 days) and the other running from 11:00:00 on 06/19/00 to 14:00:00 on 12/11/00 (~186
days). Altogether, fewer than a dozen “bad” values were “fixed.” The two sub-series of data were appended
together and treated as one period of time, lasting approximately 328 days, nearly a full year. By looking at
nearly an entire year of data, forecasts can be made and evaluated through almost all seasons. This provides
a more realistic assessment of forecast quality, since the ATST will observe the sun through all seasons, not
just the seasons favorable to air temperature predictions.

         Meteorological data for Tucson were collected by the National Weather Service (NWS) at Tucson
International Airport and then downloaded from the NWS website by Rob Hubbard. These observations
were recorded hourly and include air temperature, air pressure, wind speed, wind direction, visibility,
dewpoint, relative humidity, and cloud cover. Observations were downloaded from 00:00:00 on 03/29/04
through 10:00:00 on 06/18/04. For comparison purposes, forecasts for Tucson generated by the NWS were
also downloaded. NWS forecasts include hourly air temperature, daily maximum air temperature, daily
minimum air temperature, wind speed, wind direction, dewpoint, relative humidity, cloud cover, and chance
of precipitation. Predicted values are listed in three hour intervals, so predictions were approximated for the

23 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
hours in between by using linear interpolation. Forecasts for Tucson were posted to the NWS website once
about every 12 hours, but only the forecasts available in the early morning were downloaded. Forecasts were
downloaded from 06:00:00 on 03/29/04 through 18:00:00 on 06/20/04. There were a few small gaps, “bad”
records, and partial records in both the observation and forecast data. These anomalies were deleted in a
manner that allowed the time series to remain as continuous as possible. A few days are missing here and
there, but the sequence of hours is consistent; in other words, data recorded at midnight are always
immediately preceded by data recorded at 11pm, but possibly on a different day.

          Meteorological data for La Palma island were collected by the weather station at the Nordic Optical
Telescope (NOT), which is part of the Observatorio Roque de los Muchachos (ORM), located at an elevation
of 7815 feet (2382 meters) on La Palma island. For use in this project, the data were downloaded from the
NOT weather archive, which is publicly available on the NOT website. The NOT weather archive contains
observations from 02/03/97 to present. Observations are recorded every 5 minutes and include air
temperature, air pressure, wind speed, wind direction, and relative humidity. To make the data sets more
manageable, only observations recorded on the hour were used. Two separate time series, one from 02/11/04
to 08/11/04 (~183 days) and another from 02/03/97 to 12/31/97 (~332 days), were downloaded and
“cleaned” in a manner similar to the Tucson data. For the first NOT time series, cleaning involved replacing
a few “bad” values with “good” ones. For the second NOT time series, cleaning involved the deletion of a
few days with incomplete data, and inserting about 30 hours worth of missing data (using a rough linear
interpolation). Considering the relatively large size of these two time series, the cleaning is thought to have
little effect on forecast results.

         The table below summarizes some of the basic statistical characteristics of the air temperature time
series data for the three sites examined in this project.

Figure 4. Site statistics
                     # of data     mean T       median T      min. T       max. T       st. dev.     variance
site (year)
                      points        (°C)          (°C)         (°C)         (°C)          (°C)         (°C)
Tucson (2004)           1540        23.86         23.89        6.11         40.00         7.39        54.56
Haleakala (2000)        7600        10.59         10.78        -1.38        17.99         2.74        7.51
La Palma (1997)         7400        7.94            8.0         -3.9        20.9          4.79        22.91
La Palma (2004)         4300        7.59           7.34        -6.94        23.02         5.95        35.40

k-Nearest Neighbor Method
         The k-nearest neighbor method is a technique for interpolating and extrapolating values in time
series and other types of data sets. When using this method with time series data, one compares a sub-series
of interest, the “reference,” to all previous sub-series, the “candidates,” to find the closest matches, or
“nearest neighbors.” The conditions in and around the nearest neighbors are used as estimates or predictions
of the conditions in and around the reference. The coefficient k denotes the number of nearest neighbors that
are sought; the values they yield are averaged together to produce a single estimate. The size of the sub-
series, called the “window size,” may be varied to achieve better results.

        Consider the example illustrated in figure 1.0, where the k-nearest neighbor method is used to
forecast a value one time-step ahead in a univariate time series. In this example, k equals 2 and the window
size equals 4. The reference is the sub-series of four consecutive values leading up to the forecast time. The
candidates are all the other possible sub-series of four consecutive values. In this example, there are 12
candidates. Each candidate is compared to the reference, and the difference between them, or error, is
calculated. A common way to calculate the error is by using Euclidean difference. Since k equals 2, the two
candidates with the smallest errors are considered the nearest neighbors. To forecast the value one time step

24 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
beyond the reference, the values that appear one time step beyond the two nearest neighbors are averaged

               Figure 5. K-nearest-neighbor forecasting with window size of four and k=2.
                                               (Plummer 2000)

        An IDL procedure ( was written to generate temperature forecasts. For comparison
purposes, forecasts were generated by a variety of methods (described below). Note that x represents the
forecast horizon, measured in hours, t represents the time at which a forecast is made, measured in hours, and
d represents a number of days.

Figure 6. Forecast methods used in this project
abb.                                            description                                            sites
P24      “Persistence” from 24 hours before t. T at time t will be the same as T at time t-24.          all
P24a     Adjusted “persistence” from 24 hours before t. Assume the error at time t will be the          all
         same as the error at time t-x.
Px       “Persistence” from x hours before t. T at time t will be the same as T at time t-x.            all
Pd       “Persistence” over d days before t. T at time t will equal the average of the Ts at time t-    all
         24, t-48, … , t-24d.
C        “Climate” averaged through the years 1994-2002, excluding the year 2000.                       Hal
Ca       Adjusted “climate.” Assume the error at time t will be the same as the error at time t-x.      Hal
NWS      National Weather Service forecasts.                                                            Tuc
NWSa     Adjusted NWS forecasts. Assume the error at time t will be the same as the error at time       Tuc
NN       K-nearest neighbor forecasts.                                                                  all
NNCa     K-nearest neighbor forecasts (NN) averaged with adjusted “climate” forecasts (Ca).             all

        Nearest neighbor forecasts can be made using “raw,” un-transformed temperature data. However, a
few simple transformations improve results by reducing errors. For example, the use of transformed data
reduces both the RMSE and MAE of forecasts made for Haleakala by ~0.2°C.

25 NSO REU/RET Annual Program Report 2004                                           APPENDIX I – RET Program
         One reason transformations are needed is that air temperature time series data are generally non-
stationary. Stationary data “has the property that the mean, variance, and autocorrelation structure do not
change over time.” A stationary series is a “flat looking series, without trend, constant variance over time, a
constant autocorrelation structure over time and no periodic fluctuations (seasonality).” (NIST website) Put
another way, a time series is considered to be stationary (or to exhibit stationarity) if its statistical properties
do not vary with time. Nearest neighbor forecasts do not easily capture the broad trends or seasonal
fluctuations exhibited by non-stationary data, so the data must be made stationary by removing the trends and
seasonalities. This is accomplished by taking the “first difference” of the data. For the temperature time
series used in this project, T(t) represents the temperature (°C) at time t (hours), and its first difference, T’(t),
is found by T’(t)=T(t)-T(t-x), where x is the forecast horizon (hours). The nearest neighbor process is then
performed on the differenced temperatures, T’(t). This means that the forecasts produced by the nearest
neighbor method are actually forecasts of T’(t). Forecasts of T’(t) are easily transformed into forecasts of
T(t) by T(t)=T’(t)+T(t-x). Although calculating the first difference using x, the forecast horizon, rather than a
value of 1 may seem unnatural, it is motivated by the fact that when forecasting a noontime temperature,
T(12), 3 hours in advance, T(12-1) is not yet known. However, T(12-3) is known; it is simply the
temperature at the time the forecast is being made. Like temperature, the other meteorological time series
used in this project are non-stationary, so they are also differenced in the same manner.

         Another reason to transform the temperature data is to allow the inclusion of other meteorological
variables in the nearest neighbor process. If variables like air pressure or wind speed are to be included in
the references and candidates used by the nearest neighbor method, they must be normalized, along with
temperature, to a common interval. Otherwise, the variable with the largest magnitude or range would
dominate, perhaps creating undesirable results. All meteorological variables in this project are normalized to
the interval [0,1] before differencing. Differencing before normalization was explored, but it provided nearly
identical results.

       For the nearest neighbor forecasts, tuples of various structures were tested. The size and dimension
were varied to produce the best possible forecast. The general form of the tuples was
                               {Tt, Tt-1, … , Tt-n, Pt, Pt-1, … , Pt-n, St, St-1, … , St-n}
where T is temperature, P is pressure, S is a pyranometer reading (or other meteorological variable, like wind
speed), n is the number of hourly measurements used, and t is the time. For a reference tuple, t is the time of
the desired forecast less the horizon time (x), and for a candidate tuple, t may have any value between n and
the time of the desired forecast less 2x. Based on the literature, it was believed that no more than one or two
variables other than temperature would improve nearest neighbor forecasts. In fact, tuples consisting solely
of 24 consecutive hourly temperature measurements produced the best forecasts for all three sites. Using
fewer than 24 temperature measurements increased forecast error, and using more than 24 temperature
measurements reduced forecast error, but by insignificant amounts. For Haleakala, the use of air pressure
and pyranometer readings was explored. The inclusion of 24 consecutive hourly air pressure measurements
reduced both RMSE and MAE by only ~0.03°C. The inclusion of 24 consecutive hourly pyranometer
readings actually increased RMSE and MAE by ~0.25°C. For La Palma, the use of air pressure and wind
speed was explored, but neither significantly reduced forecast error. For Tucson, no variables other than air
temperature were used in the nearest neighbor forecasts.

        The value of k, the number of nearest neighbors used to make a single prediction, was varied
between 1 and 100. Higher k values provided lower forecast errors. Generally, for k values greater than 40,
the reduction in error became negligible, so unless otherwise noted, all reported results were achieved with
        Candidate tuples are compared to reference tuples using Euclidian distance, which is given

                                         (r1-c1)2 + (r2-c2)2 + … + (rn-cn)2

26 NSO REU/RET Annual Program Report 2004                                                     APPENDIX I – RET Program
where rn is the nth element in the reference tuple and cn is the nth element in the candidate tuple. The k
candidate tuples with the least distances from the reference tuple are considered the k nearest neighbors. The
k predictions yielded by the k nearest neighbors are then averaged together to provide a single temperature
prediction at a given time. Weighted averages are sometimes used in the nearest neighbor process, but they
were not attempted in this project.

Ret Project, Summer 2005?
Once the NSO has perfected a temperature forecasting scheme, perhaps some effort should be geared toward
developing methods to forecast earthquakes (Big Bear Lake), volcanic eruptions (Haleakala, La Palma), and
mega-tsunami generating landslides (La Palma).

A. Aussem, F. Murtagh, M. Sarazin, Dynamical Recurrent Neural Networks: Towards Environmental Time
Series Prediction, International Journal of Neural Systems 6(2):pp. 145-170, 1995,

F. Buffa and I. Porceddu, Temperature forecast and dome seeing minimization I. A case study using a neural
network model, Astron. Astrophys. Suppl. Ser. 126, 547-553 (1997),

Nathan E. Dalrymple, Mirror Seeing, Project Documentation Report #0003 Revision #A, 9/26/02,

Nathan E. Dalrymple, Primary Mirror Cooling Concepts, Project Documentation Technical Note #0020
Revision #B, 3/20/03,,

D. Andre Erasmus, Meteorological Forecasts for the Southern African Large Telescope (SALT):
Phase 1; Setup and Maintenance of an Automated System for Operational Forecasts, Interim Report, 31st
January, 2004,

Gregor Giebel, The State-Of-The-Art in Short-Term Prediction of Wind Power, A Literature Overview,
Version 1.1,

NIST website, Introduction to Time Series Analysis,

27 NSO REU/RET Annual Program Report 2004                                           APPENDIX I – RET Program
NOT weather website,

Eric A. Plummer, Time Series Forecasting With Feed-Forward Neural Networks, Guidelines and Limitations
(Master’s Thesis),

Christoph C. Raible, Georg Bischof, Klaus Fraedrich, and Edlbert Kirk, Meteorologisches Institut,
Universitat Hamburg, Hamburg, Germany, “Reply,” Weather and Forecasting, Aug. 2, 2001, Volume 16, p.

SALT weather website,

Other Relevant Information Sources
Selected astroclimatology papers at the ESO,

M. Sarazin bibliography, [site-surveys and astroclimatology],

M. Sarazin, ESO Internal Workshop on Forecasting Astronomical Observing Conditions, 29-30 May, 1997
Workshop Executive Summary, Predicting Observing Conditions at ESO Observatories, reality and
perspectives, October 7, 1997,

Lorenzo Zago, The Effect of the Local Atmospheric Environment on Astronomical Observations, PhD
thesis, Ecole Polytechnique Fédérale de Lausanne, February 1995,

Kevin Roe, Duane Stevens, One Kilometer Numerical Weather Forecasting to assist Telescope Operations,

D. Andre Erasmus website,

D. Andre Erasmus, C. A. van Staden, Meteorological Site Evaluation and Forecasting needs for the Southern
African Large Telescope (SALT),

D. Andre Erasmus, Relating Seeing Quality To Meteorological Conditions: Development of Seeing Quality
Forecasts and Improvement of Site Selection Procedures, 1988 and 1996,

Ignazio Porceddu, Cagliari Astronomical Observatory, Centro Meteorologico Territorial de Canarias
Occidental, “Astronomy And Meteorology II, The TNG interdisciplinary activity”,

Wojciech Kowalczyk, Averaging and data enrichment: two approaches to electricity load forecasting,

Time-Critical Decision Making for Economics and Finance, [general forecasting info],

28 NSO REU/RET Annual Program Report 2004                                        APPENDIX I – RET Program
An Introduction to Neural Networks, Prof. Leslie Smith, Centre for Cognitive and Computational
Neuroscience, Department of Computing and Mathematics, University of Stirling, Scotland, April, 2003,

Neural Networks, by Christos Stergiou and Dimitrios Siganos, Imperial College, London, England,

Neural Networks on,

Kalman Filters on,

How is a weather forecast made? (NWS) -
Predicting highs, lows, and temperature trends -
Weather Computer Models -
Some operational forecast models -
forecast verification

Chatfield, C., The Analysis of Time Series. Fifth Edition, Chapman and Hall,1996.
Efron, B., and R. Tibshirani, An introduction to the bootstrap. Chapman and Hall, 1993.
Rousseeuw, P.J., and A.M. Leroy, Robust regression and outlier detection. John Wiley, 1987.

Meteorological Data Sources

Hawaii Weather/Climate Modeling Ohana (featuring MM5 forecasts) -
Haleakala Site Survey Home Page (weather info & data) -
Haleakala Ranger Station climate data -
Haleakala Ranger Station “point forecast matrices” -
NWS forecasts for HI (click “Experimental Text” under “Forecasts”) -
Haleakala National Park weather page (by the NPS) -

Big Bear:

Big Bear Lake WRCC climate data -
Big Bear Lake NWS climate data (click on map) -
Big Bear City, CA, “point forecast matrices” -
Fawnskin, CA, weather observations – (click Fawnskin on map), weather - -
Ben Brissey’s Weather, Big Bear, CA -
Bill Marquette -

La Palma:

Instituto de Astrofísica de Canarias (IAC) weather page -
Instituto Nacional de Meteorología (INM) forecasts -
INM forecasts for Canary Islands -

29 NSO REU/RET Annual Program Report 2004                                         APPENDIX I – RET Program
IAC Sky Quality Group (site survey) -
Nordic Optical Telescope (NOT) weather data (& archive) -


Point Forecast Matrices (predictions) -
Digital Forecast Products (prototype predictions) -
Hourly Observations (aviation format) -
Climate data -


ECMWF homepage -
National Center for Environmental Prediction (NCEP) homepage -
National Center for Atmospheric Research (NCAR/UCAR) homepage -
National Oceanic and Atmospheric Administration (NOAA) homepage -
National Weather Service (NWS) homepage -
Mesoscale Model 5 (MM5) homepage -
Western Regional Climate Center (WRCC) (see “Historical Climate Information”) -
Climate Data & Weather Observations -
NNDC Climate Data Online -

30 NSO REU/RET Annual Program Report 2004                                        APPENDIX I – RET Program
                                                                  Haleakala (2000): RMSE vs. Horizon

                                                                                                                     P24 RMSE
    RMSE (degrees C)

                                                                                                                     P24a RMSE
                                                                                                                     Px RMSE
                                                                                                                     C RMSE
                                                                                                                     Ca RMSE
                                                                                                                     NN RMSE
                              0   1   2   3   4   5   6   7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
                                                                      Horizon (hours)

                                                                  La Palma (1997): RMSE vs. Horizon

    RMSE (degrees C)

                                                                                                                     P24 RMSE
                                                                                                                     P24a RMSE
                                                                                                                     Px RMSE
                                                                                                                     NN RMSE
                              0   1   2   3   4   5   6   7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
                                                                      Horizon (hours)

31 NSO REU/RET Annual Program Report 2004                                                                  APPENDIX I – RET Program
                                                                La Palma (2004): RMSE vs. Horizon

  RMSE (degrees C)

                                                                                                                     P24 RMSE
                                                                                                                     P24a RMSE
                                                                                                                     Px RMSE
                                                                                                                     NN RMSE
                            0   1   2   3   4   5   6   7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
                                                                   Horizon (hours)

                                                                Tucson (2004): RMSE vs. Horizon

                                                                                                                    P24 RMSE
  RMSE (degrees C)

                                                                                                                    P24a RMSE
                                                                                                                    Px RMSE
                                                                                                                    NWS RMSE
                                                                                                                    NWSa RMSE
                                                                                                                    NN RMSE
                            0   1   2   3   4   5   6   7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
                                                                   Horizon (hours)

32 NSO REU/RET Annual Program Report 2004                                                                   APPENDIX I – RET Program
                                 National Science Foundation
                              Research Experiences for Teachers
                                    Summer 2004 Report

                                            Michael Sinclair
                             Kalamazoo Area Mathematics & Science Center
                             National Solar Observatory/Sacramento Peak

                                        Advisor: Dr. Joel Mozer

                 Deconvolving the Point Spread Function
                from Solar Mass Ejection Imager Imagery
       The removal of optical artifacts from the United States Air Force Solar Mass Ejection
       Imager’s data is required to correct a “smearing” of stellar images in order to complete
       precise photometric analysis of stars. Employing a multistage approach, first by mapping
       each pixel within a given camera array using polar coordinates onto an oversized grid and
       then translating that data into a corrected rectangular array, we are able to create a more
       accurate map of the sky. Further image processing by deconvolution of the point spread
       function from the rectangularized array produced corrected stellar images that show point-
       like characteristics of actual stars. With this rectangularization and deconvolution technique,
       one can then use the SMEI images for additional astronomical studies of stellar variability
       and periodicity as well as imaging and studying novæ and supernovæ.

1.   Introduction.

     The Solar Mass Ejection Imager is a satellite-based experiment designed to image the entire
sky in white light with respect to the Earth in order to identify and track coronal mass ejections
propagating through space. By doing this, it is able to provide one- to three-day forecasts of
approaching geomagnetic storms, allowing initiation of preventative measures in order to mitigate
damage to satellites and power grids on and around the Earth. Launched on January 6, 2003 from
                              Vandenberg Air Force Base in California, the Coriolis Mission carries
                              two camera experiments: the U.S. Navy Windsat and the U.S. Air
                              Force SMEI. The satellite circles the Earth in a Sun-synchronous orbit
                              at an altitude of approximately 840 km (or 504 miles) above the
                              surface, roughly along the terminator. In addition to its primary
                              mission of imaging coronal mass ejections, the Solar Mass Ejection
                              Imager can also be used to study heliospheric shock waves, solar
                              phenomenon which co-rotate with the Sun, the Zodiacal light, near-
                              Earth objects, asteroids and comets, transits of bright stars by their
                              planets, stellar variability, and discoveries of novæ and supernovæ.

33 NSO REU/RET Annual Program Report 2004                                          APPENDIX I – RET Program
2.   Satellite Configuration and Instrumentation.

     The Coriolis satellite pictured at right, which contains the U.S.
Navy Windsat microwave polarimetric radiometer and the Air Force
Solar Mass Ejection Imager, was assembled and launched by the
U.S.A.F. Space and Missile Systems Center’s Detachment 12 Space
Test Program Office at Kirtland Air Force Base in Albuquerque, New
Mexico using a refurbished Lockheed-Martin Titan 2 booster and a
Spectrum Astro spacecraft bus. Spectrum Astro integrated the
spacecraft bus and payloads with the Titan booster prior to its launch
from Vandenberg Air Force Base. The total mass of the Solar Mass
Ejection Imager is 36 kg and uses an average of 38 W of power.
     The two major components of the Solar Mass Ejection Imager are
an electronic optical system and a data-handling unit. The optical array
is shielded by a baffle system designed to minimize stray light. The
three SMEI baffle systems are attached at the bottom of the spacecraft
in the section which was despun in orbit. The baffles can be seen
projecting from this section in the photograph at right.
     The camera baffles or “sunshades” were designed by the
University of California at San Diego and the University of
Birmingham. The baffles reduce the amount of scattered light reaching the CCD cameras by
approximately ten orders of magnitude when viewing near the Sun. For imaging faint heliospheric
objects (such as coronal mass ejections), a 15-fold rejection of scattered light is required and the
remaining magnitude reduction is accomplished within the optical system as described on the
following page.

     The SMEI imaging array consists of three CCD cameras as shown in the diagram above. Each
camera covers a different 3˚ 60˚ field of view such that, together, they view a thin 180˚-wide slice
of the zenith-facing hemisphere of the sky, sweeping over the entire sky through the course of each
101-minute orbit. A shutter/sensor system protects the detectors from direct solar illumination.

34 NSO REU/RET Annual Program Report 2004                                     APPENDIX I – RET Program
     The camera and optics were designed and constructed by
the University of California at San Diego. As shown in the
exploded view of the SMEI optics at right, the entrance
aperture from the baffle is Z0. Two mirrors, labeled M1 and
M2, reflect the incoming light to the CCD detector plane. The
vane V restricts the amount of light entering the image plane
and serves as a Lyot stop of a coronagraph, thus reducing
additional scattering. The primary mirror, M1, has a spherical
surface which reflects the light to the angled flat (M2) which,
in turn, reflects and focuses the image on the CCD.
     The camera electronics and data handling unit (DHU)
were designed by the University of Birmingham and
subcontracted to Rutherford Appleton Laboratories. The
CCD is an EEV CCD05-30-231A sensor and is a frame
transfer device with 1242 576 pixels (22.5 m).
     The DHU controls the three on-board cameras and
interfaces between the Coriolis spacecraft and SMEI. The
power supply printed wiring board is structured around Texas
Instrument’s SMJ320C5AHFGM-40 DSP microprocessor
operating at 19 MHz. The flight software controls all data
transferred between the cameras, data compression, and data transfer from SMEI to the spacecraft
telemetry system, which also handles operational command to SMEI from the ground.

3.   Converting Image Arrays.

     As described earlier, each camera frame covers a roughly 3˚ 60˚ section of the sky. The figure
below shows three simultaneous slices from each camera. Camera 1 shows the night sky in the anti-
Solar direction while Camera 2 is
centered along the terminator and
Camera 3 points closest to the Sun. Each
region is arrayed as an arc rather than as
a rectangle. Together, the three camera
fields-of-view, when “stitched” end to
end, correspond to an approximate great
circle in the sky. When these frames are
mapped into ecliptic coordinates, a
complete sky map can be built up over
an orbit (except for small areas centered
on the Sun and in the anti-Solar
direction). One such map, produced in
March 2003, is shown on the following
page. This is an equal-area Hammer-Aitoff projection centered on the Sun with the north ecliptic
pole at the top and the south ecliptic pole at the bottom. The dark circle at the center of the picture is
a zone of exclusion 20˚ in radius centered on the Sun. Zodiacal light can be seen on either side of
the solar disk, Venus is to the right just off the limb of the Sun, the Milky Way galaxy appears to be
“wrapped” around the Sun, and the Magellanic Clouds are at the bottom of the image.

35 NSO REU/RET Annual Program Report 2004                                        APPENDIX I – RET Program
     One of the key difficulties in working with a Solar Mass Ejection Imager array is the non-linear
structure of the image. In order to perform photometric analysis of stars, the primary focus of this
project was to redefine the curved 3˚ 60˚ array into a corrected rectangular array. This
rectangularization effectively removes a geometric rotation in the coma distortion that is inherent to
SMEI imagery. From there, one can perform a deconvolution routine which will “clean up” and
sharpen stellar images by removing the coma-dominated point spread function.
     The first step of this project was to convert the raw curved array (as shown below) into a simple
rectangular array using a brute-force approach. This required a shift of all data from a given xarc,yarc
position in the arc to an assigned xrectangular,yrectangular location.

                                           Raw curved array.
    As can be seen in the following array, however, it leaves a number of problems in the final
image. Besides added noise and image displacement, the most important difficulty is the extended
nature of a given star.

                                        Raw rectangular arrayed

36 NSO REU/RET Annual Program Report 2004                                      APPENDIX I – RET Program
     This can be seen in the figure at right, where a single star image
has been pulled out of the rectangular array and enlarged. Notice the
extended “fish-like” and tilted shape of the star as well as the
considerable pixelation of the image. Of course, stars do not look like
this. Therefore, another approach was required in order to properly
orient and smooth the image prior to deconvolving it.
     The next approach in converting the curved arrays into a
rectangular format was to attack the problem by calculating the
position of each pixel in the SMEI CCD image plane with respect to
the center of the curved array. In this case, a translation using the
polar coordinate position with the relative radius and angle from the center of each camera created a
mapped array over the SMEI-generated arc. Applying a bilinear function for each mapped pixel,
one could then interpolate the pixel position from the curved array to a corrected rectangular array.
This, in turn, would produce an accurate xi,yi location for every pixel with respect to its radius and
angle from any one of the three given satellite cameras.
     Thus, the two-step reformulated rectangular array would look like this:

                                    Corrected raw rectangular array.
This corrected array is clearly cleaner and shows a more uniform shape for each star. In particular,
the axes of each “fish” are all in the same direction.
     Finally, the rectangularized array should be deconvolved with the point spread function in order
to sharpen up all of the individual stellar images. Deconvolution is an
image processing technique that removes features in an image that are
caused by the telescope or the arrangement of the optics in the telescope
rather than from actual light coming from the sky. The IDL program
max_entropy, originally developed by Frank Valosi at NASA’s Goddard
Space Flight Center, was used to produce a final corrected image.
Deconvolution requires two steps in order to improve image quality.
First, a reasonably bright star is selected to serve as a baseline for the
point-spread function. An example of such a star image is shown at right.
Notice how much smoother and more uniform the image is as compared
to the previous example of an individual star shown above. There is, however, a slight
triangularization of the image but this can be removed by applying a rotation of 180˚ to the image
and running it through the deconvolution process. This image will serve as a kernel for the
deconvolution procedure.
     Mapping the kernel with respect to any section of the rectangular array will then produce a
clean, smooth point image of a star with very little noise. Therefore, in the second stage of the
deconvolution procedure, the kernel cycles through every pixel position in the rectangular array and
“focuses” the image in order to improve it. An example of the deconvolved SMEI array is shown

                                   Deconvolved rectangularized array.
     Notice that the stars are now clearer and crisper in shape. The relative brightness of each star is
also substantially more uniform and there are no saturated stellar images.

37 NSO REU/RET Annual Program Report 2004                                       APPENDIX I – RET Program
     A careful examination of a single star pulled out from the deconvolved array shows a distinct
and obviously more star-like structure. A dramatic decrease in pixelation
also improves the quality of the image. An example of a star taken from
the deconvolved array on the previous page is shown at right. Notice the
definite point-like structure of the star, as would be expected. The
pixelation has also been minimized and the image quality has been clearly
improved as well.
     In addition, a contour plot of the pixel brightness of the star and the
surrounding region would more accurately show whether or not the stellar
image is, in fact, properly circularized. Upon examination of such a plot,
only a slight ellipticity of the image is observed.
     This can be readily seen in the illustration shown at right. It is obvious
that any further correction of the array would only slightly improve the
image, although that step may be necessary in order to allow for extremely
precise photometic analysis of individual stars. One possible extension to
this project would be to continue to improve the circularization of the final
stellar image.
     The IDL program was used to rectangularize the Solar
Mass Ejection Imager curved array and can be found as an Attachment to
this paper.

4.   Extensions.

     A future application using this rectangularization technique would be to extract raw data from
the Solar Mass Ejection Imager archives and photometricly analyze a series of “guide” stars. This
would allow for calibration of the SMEI optics and evaluate the “seeing” ability of SMEI in order to
compare actual magnitude seeing with that of its projected capabilities (which were expected to be
on the order of 12th magnitude). A re-drawing of the all-sky map would also be a useful and
interesting extension as it would provide astronomers with a more accurate representation of the
night sky and would allow for easy identification of near-Earth objects such as asteroids or comets,
since those objects would appear as smeared stars rather than points of light. Finally, a careful and
long-term examination of active SMEI images may allow astronomers to discover and study short-
and long-term stellar phenomenon, such as variable stars, novæ, and supernovæ.

•    Conclusions.

     A careful mathematical mapping of the raw curved Solar Mass Ejection Imager data using a
translation from polar coordinates (relative to each CCD camera) to rectangular coordinates
combined with a bilinear function can produce a more accurate rectangularized array. With a
carefully selected baseline guide star to serve as a measure of the point spread function, this
function can be used to deconvolve the entire rectangularized array to successfully sharpen and
improve the SMEI data. This, in turn, will allow astronomers to apply more precise photometric
analysis to stellar objects, opening the door to using the SMEI as a robotic telescope studying stellar
periodicity or variability, discovering novæ and supernovæ, and searching for near-Earth objects
such as asteroids and comets.

38 NSO REU/RET Annual Program Report 2004                                     APPENDIX I – RET Program

     I would like to extend my sincere thanks and appreciation to Dr. Joel Mozer. I would also like
to thank Timothy “Tex” Henry and Cheryl-Annette Kincaid for their able assistance throughout my
apprenticeship with IDL programming. Finally, I wish to express my deep appreciation to the entire
Sunspot community at Sacramento Peak for a wonderful summer.
     This research was funded by the National Science Foundation (NSF) Research Experience for
Teachers (RET) program at the National Solar Observatory at Sacramento Peak, New Mexico.


Beck, Rainer,, Solar Astronomy Handbook; Richmond, VA: Willmann-Bell, 1995.
Fanning, David W., IDL Programming Techniques; Boulder, CO: Fanning Software, 1997.
Mizuno, D., “Some Results for SMEI Brightness Calibration”, Boston College, 6 July 2004.
No author given, Building IDL Applications; Boulder, CO: Research Systems, 1991.
No author given, The IDL Reference Guide; Boulder, CO: Research Systems, 1991.
No author given, Using IDL; Boulder, CO: Research Systems, 1991.
Wentzel, Donat G., The Restless Sun; Washington, DC: Smithsonian Institution Press, 1989.

39 NSO REU/RET Annual Program Report 2004                                   APPENDIX I – RET Program

     The program to examine a SMEI .buf file and convert it to a corrected rectangular array and deconvolve
the stellar images is given below:

; Name:
; Purpose:       To convert "arc" array into rectangular array using polar
;               coordinates of each pixel from the frame data, pull out an
;               individual star image for a baseline threshold and point-
;               function, and deconvolve the array so that stars appear as
;               in order to be analyzed for satellite degradation and stellar
;               periodicity or variability (if any).
; Input:        1.   SMEI file (L1A...)
;               2.   camera # (n)
;               3.   lower-left x-value         (raw)
;               4.   lower-left y-value         (raw)
;               5.   lower-left x-value         (dcv)
;               6.   lower-left y-value         (dcv)
;               7.   pixel threshold
; Output:       IMAGES.
;               1. raw rectangularized array
;               2. convolved rectangular array
;               3. deconvolved rectangular array
; By:                             M. Sinclair, NSO/SP 2004
;                              F. Varosi, NASA/GSFC 1992
;                  with IDL 5.0 conversion                 W. Landsman, AIT/EKU 1997
;               max_entropy link                             J. Moser, AFRL/SP 2004
;               read data & location                    C. Kincaid, NSO/SP 2003
;               pgm design & assistance                 T. Henry, NSO/SP 2004

pro sinclair,infil,n,frame_data=frame_data
    ; Set the basic parameters for converting the raw SMEI data

   infil= ''
   read,'Input <L1A...filename>:',infil
   read,'Camera #?',n
           for i=0,317 do begin
                 w=where(jnput[i,*] LE 8250)
                 if w[0] GE 0 then input[i,w]=jnput[i,w]
    ; Create simplified rectangular array with raw data by "shifting"


40 NSO REU/RET Annual Program Report 2004                                         APPENDIX I – RET Program
            for i=0,nx-1 do begin
                  w=where(frame_data(i,*) LT 10000)
                  if w[0] ge 0 then begin

    ; Constants

      if n EQ 1 then begin
      if n EQ 2 then begin
      if n EQ 3 then begin
    dtheta=0.10                                   ;arc-sec/pix
    dy=1.0                                   ;delta-r/pix

    ; Set radius and angle
    theta=(findgen(751)*dtheta-35.0)*rpd     ;theta
    r=reverse(findgen(101))*dy+1049./4.           ;r

    ; Map array over SMEI arc and convert to rectangular array

   for i=0,750 do begin
      for j=0,100 do begin
           p[i,j]=camnx+r[j]*stheta[i]            ;xpos in data
           q[i,j]=camny-r[j]*ctheta[i]      ;ypos in data
      if (pij GT 0) and (pij LT 317) and (qij GT 0) and (qij LT 63) then

41 NSO REU/RET Annual Program Report 2004                APPENDIX I – RET Program

    ; Pull out rectangular array and enlarge

    ; Pull out star image from rectangular array

       print,'click on center of selected star'
            ; a=''
            ; b=''
            ; read,'select star x-value:',a
            ; read,'select star y-value:',b
            ; rdpix,star_image

    ; Smooth selected star image

            ; tv,star_pix

    ; Convolve rectangular array with star_data

    ; Deconvolve the array using max_entropy procedure

       if N_elements(niter) EQ 0 then niter=10
                  for i=1,niter do begin

    ; Examine deconvolved star image


42 NSO REU/RET Annual Program Report 2004                APPENDIX I – RET Program
      read,'select star x-value:',c
      read,'select star y-value:',d
      read,'threshold (pixel value):',THr
        tvscl,final_star GE THr


43 NSO REU/RET Annual Program Report 2004                APPENDIX I – RET Program
                           NSF Research Experiences for Teachers
                                   Summer 2004 Report

                                        Creighton G. Wilson, Sr.
                                  Lovelady High School, Lovelady, Texas

                                National Solar Observatory/Sacramento Peak
                                        Advisor: Dr. Alexei Pevtsov

               The Neutral Line Method of Flare Prediction
                          from Magnetograms

        The purpose of this research was to determine if there is a relationship between the length of
        neutral lines and flare production from active regions. Active regions were chosen from the
        Solar-Geophysical Data prompt reports based on the region originating on the face of the
        Sun. From these same reports the active region was track through the magnetograms to
        determine the number of days the region could be observed before it reached the far limb of
        the Sun. SOHO Sun images were then requested through Stanford University from their SOI
        Data Sets. An IDL program, image2_mod, was written and used to mark the neutral lines
        for all images available for the chosen duration of each active region. A ratio of the length
        of the neutral line and the diameter of the active region determined the normal_length of the
        active region. The number of flares and the flare index was calculated using another
        newly written program, extract, which obtain data from tables of the National
        Geophysical Data Center. The last program to be written, superextractor4, removed
        duplicate data, and flawed data, before plotting a graph of average normal_length versus
        flare index to determine if there was a relationship between the two.

•   Sequential Historical Progression of Sunspots Analysis

         Galileo first introduced the concept of sunspots in 1609. His observations, though simple, were
revolutionary. Nearly 200 years passed before anything much more than the fact of their existence was
known. Based on the discovery of narrow dark lines in the solar spectrum by William Wollaston in 1802,
Joseph von Fraunhofer continued the study using high magnification of the Sun and discovered over 300
gaps in the solar spectrum by 1814. Laboratory experiments in the 1800's revealed that every element has its
own characteristic spectrum. This came to be known as the “fingerprint” of the element. Taking these
findings and applying them to Fraunhofer's spectrum of the Sun resulted in the conclusion that vaporized iron
existed on the Sun. The next major development in the study of sunspots came with the work of Charles
Young at Halstead Observatory. Young applied spectroscopy to the study of sunspots. Young discovered
double lines in the spectra of the spots, and hypothesized that they were the results of layers of hot gases.
This explanation was shown to be incorrect by the experiments of Pieter Zeeman in 1886. Zeeman
discovered that the emission lines of the spectrum of a sodium flame were split when the flame was placed
between two magnets. More tests revealed that the amount of splitting was related to the strength of the

44 NSO REU/RET Annual Program Report 2004                                           APPENDIX I – RET Program
magnetic field through which the light passed. By 1905 George Ellery Hale had applied Zeeman splitting to
the study of sunspots using the Mt. Wilson Observatory. Based on the principles of Zeeman splitting he
proved that there are magnetic fields associated with sunspots. See Zeeman Splitting image below. Zeeman's
work and Hale's application of it to sunspots provides the basis for present study and analysis of sunspots.
From these experiments tremendous amounts of information are able to be derived.

                                              Zeeman Splitting

•   Magnetic Fields and Magnetographs

         Magnetic fields in sunspots may be studied using either a longitudinal magnetograph or a vector
magnetograph. The study of magnetographs is much like the analysis of motion. In motion speed is
concerned with magnitude only. Velocity is concerned with magnitude and direction of motion. Magnitude
of motion is expressed by the length of the arrow and direction is expressed by the way in which the vector
arrow point is pointing. This is much the way it is in vector magnetographs analysis. Magnetograph
measure both size and direction. The longitudinal or line-of-sight magnetograph measures only the size of
the magnetic field coming toward a person or going directly away from a person. The vector magnetograph
measures both the size of the component of the field in the line-of-sight, the size of the component directed
across the line-of-sight, the “transverse” component, and the direction of the transverse component. When
vector magnetic field measurements are available, a map can be reproduced that indicates the magnitude and
direction of the actual magnetic field under consideration. The dipole field mapping (first image below) is a
result of analysis of light passing through the magnetic field shown in the image below the map.

                                           Dipole Field Mapping

45 NSO REU/RET Annual Program Report 2004                                          APPENDIX I – RET Program
                 Magnetic Field through Which Light Was Passed for Vector Field Analysis

•   Magnetographs, Sunspots, and Flare Prediction

          It has been discovered that polarized light at the wavelength of 5250.2 A is the best wavelength for
the determination of magnetic fields in sunspots and the construction of magnetographs. This wavelength is
very sensitive to magnetic fields, Zeeman splitting, and it originates from the photosphere where sunspots are
present. This polarized light from the sunspot is made up of a circular part and a linear part. The circular
part produces the line-of-sight component of the magnetic field. The linear polarized part produces the
vector component of the magnetic field, providing the magnitude and direction of the transverse to the line-
of-sight. Therefore, based on the light analysis of a sunspot, a map of the magnetic field of the sunspot can
be reproduced.
          This mapping has shown that magnetic fields normally rise up out of the photosphere of Sun at one
point of the field, continues into the corona, and then return to the Sun at the another point on the
photosphere. Sometimes the magnetic field becomes twisted or sheared by fluid motion of the Sun as it
rotates. Continual analysis of the magnetic fields by magnetographs quantifies the amount of twisting and
shearing. From that data various correlations can be made to contribute to the science of flare prediction.
          Flare predictions science has shown a relationship between the complexity of active regions and the
number of flares produced. The question that has not been answered is whether the complexity is caused by
process below the photosphere which rises to the surface or is the complexity caused on the surface over time
as the active region moves? The more complex the active region is, the longer and more complex the neutral
line is that runs from one side of the region to the other. If early in the existence of the active region
complex lines are found, then this evidence would contribute to the hypothesis that the complex is form by
forces and process occurring below the surface of the photosphere. If the complexity occurs over time then
the hypothesis of movement on the surface is of greatest probability for their formation. The related question
that this experiment addresses is the relationship of this neutral line length to flare production. Whether the
line first appears as long and complex or grows in length over time is not the question, but rather does the
length itself have in correlation with the number of flares produced and the flare index of an active region?

46 NSO REU/RET Annual Program Report 2004                                           APPENDIX I – RET Program
•   Neutral Lines, Flare Index, IDL Programs

        Magnetographs are used to form magnetograms, a picture-like image of the Sun that shows the
strength and location of magnetic fields. See image below.

Through the development of two IDL programs, image2_mod and MDI_EFR by Dr. Alexei Pevtsov of the
National Solar Observatory at Sunspot, New Mexico, analysis of magnetograms to determine the length of
neutral lines as well the size of the active region and various other information was possible. A description
of MDI_EFR which was called by image2 mod follows below.


;   for a given MDI magnetogram, extract area of emerging flux and; calculate various parameters

;   mdi_efr,filename

;    filename - name of fits file with MDI full disk magnetogram
;    outfile - name of output file (def: 'efr_parameters.dat'
;    arh - 2-elements array (Carrington longitude and latitude),
;        coordinates of active region. If arh is absent, program
;        will ask to mark position of active region on image.
;    noaa - NOA

outfile is created or appended.
;     It contains following parameters:
;    1. date and time of magnetogram
;    2. working directory, filename
;    3. xy- coordinates of center or AR, and lower-left/upper-right; corners of selected box.
;    4. NOAA AR number
;    5. heliographic, and Carrington coordinates of central area (emerging region) for magnetogram.
;    6. total and net magnetic flux of the region, area, footpoint
;     separation, title, leading polarity
;    7. length and footpoint separation of the main neutral line,
;    xy-coordinates of two ends of the neutral line
;    8. n-elements in neutral line, xy-coordinates of the neutral
;     line, and magnetic flux gradient across the neutral line.

47 NSO REU/RET Annual Program Report 2004                                            APPENDIX I – RET Program
Using image2_mod data was collected for active regions 8103, 8107, 8115, 8116, 8118, 8119, 8226, 8236,
8760, 8824, 8828, and 8829. The next step was to write the program extractor. This IDL program retrieved
flare data on the active regions being studied from the National Geophysical Data Center, and calculated the
flare index based on the duration of the flare and the magnitude assigned to the type of flare produced by the
active region. Another program, readdata, was also written for the purpose of retrieving the data that
image2_mod had obtained and saved. There were approximately twelve images of each active region per
day and as many as ten days of viewing per active region. The final program written, superextractor4,
pulled together data from extract and readdata, calculated the average normal_length and plotted a graph of
the average normal_length to the flare index. The normal_length is the length of the neutral line divide by
the diameter of the active region. It also eliminate false data.

•   Graph, Correlation

       Superextractor4 plotted the graph above and showed a correlation of 0.221169 for the relationship of
normal_length to flare index. That is a very low correlation. Of the active regions that were observed, three
produced no flares, and only three produced flares having a flare index greater than 0.5. Nearly all flares
observed fell into the 0.5 range.

•    Flare Prediction and Conclusions

         The purpose of this experiment was to determine if normal_length of neutral lines in active region
could be used to predict flare production. The data presented here shows such a low correlation that one of
two conclusion must be drawn. Either neutral lines do not predict flare production or the data is not
sufficient to draw a conclusion either way. Due to the small number of active regions observed that
produced flares with flare indexes greater than 0.5, it is suggested that the experiment be repeated with a
much larger number of active regions being measured before a firm conclusion is drawn.

48 NSO REU/RET Annual Program Report 2004                                           APPENDIX I – RET Program

        I am so grateful to Dr. Alexei Pevtsov for his participation in the RET Program and for having
accepted me to tutor. It is truly an honor to be his student. Dr. K.S. Balasubramaniam has also provided
programming assistance that has been very valuable to the completion of my research. On a number of
occasions Mr. Tim Henry has gone beyond the call of duty to spend hours explaining UNIX and IDL
programs. Jeneen Sommers, Science Data Coordinator, SOI/MDI, Stanford University, gave of her own time
to make my data requests a priority when the system was on overload, just because she wanted to. I would
be amiss if I also did not thank student SRA Cheryl Annette Kincaid and Joel Lamb for providing feedback,
advice, and opinions on my programming and research.
        My research was funded by the National Science Foundation. Jackie Diehl and Rebecca Coleman
have taken care of all the administrative needs to make the research experience possible and enjoyable.


Adams, Mitzi, Vector Magnetographs, NASA/Marshall Space Flight Center, Huntsville, AL (2002).

Anderson, Gail & Paul, The UNIX C Shell Field Guide, Englewood Cliffs, N.J., (1986).

Building IDL Applications, Research Systems, Inc., Version 5.1, Boulder, (1998).

Chaisson, Eric, and McMillan Steve, Astronomy Today, The Sun, p.353-376, Upper Saddle River, New
Jersey, (1999).

Falconer, D.A. (2001): Journal of Geophysical Research, 106, p.25,185-25,190.

IDL Quick Reference Guide, Research Systems, Inc., Boulder, (1999).

IDL Reference Guide, Research Systems, Inc., Version 5.1, Vol. 1 & 2, Boulder, (1998).

Interactive Data Language Reference Guide, Research Systems, Inc., Version 4, Vol. 1

Kaler, James B., Astronomy!, The Sun, p.310-331, Reading, Massachusetts, (1997).

Loughridge, Michael, Solar-Geophysical Data prompt reports, 629-676, (1997-2000).

National Geophysical Data Center, National Environmental Satellite, Data, and Information Service, (1997,

Solar Oscillations Investigation, W.W. Hansen Experimental Physics Laboratory, Stanford University and
the Solar and Astrophysics Laboratory of the Lockheed- Martin Advanced Technology Center.

SOHO, SOHO (Solar & Heliospheric Observatory) project is being carried out by the European Space
Agency (ESA) and the US National Aeronautics and Space Administration (NASA), (1995).

Scherrer, Deborah, The Sun – A Magnetic Star, Stanford University, (1997)

Using IDL, Research Systems, Inc., Version 5.3, Boulder, (1999).

Wentzel, Donat G., The Restless Sun, Sunspots and Magnetism, p.81-98, Washington, (1989).

49 NSO REU/RET Annual Program Report 2004                                          APPENDIX I – RET Program
50 NSO REU/RET Annual Program Report 2004   APPENDIX I – RET Program
                                        APPENDIX II

                                 REU Student Reports

51 NSO REU/RET Annual Program Report 2004              APPENDIX II – REU Student Reports
                        NFS Research Experiences for Undergraduates
                                   Summer 2004 Report

                                            Frances S. Edelman
                                              Yale University

                                    National Solar Observatory/Tucson
                                          Advisor: Dr. Frank Hill


         Oscillation amplitudes are observed to be reduced in regions of strong magnetic field, and abrupt
changes in the magnetic field are observed during the course of strong flares. However, it is not clear
whether both of these effects are due to actual changes in the solar phenomenon, or to variations in the
diagnostic used to measure velocity and magnetic fields. Most instruments do not record the full spectral
line profile, so changes in the shape of the spectral line used for the observations can produce spurious
results. Ultimately, studying this issue will help improve observations.

       Most of the work involved writing and running computer programs in order to produce models that
would help to answer these questions.

          All the data used in this project were obtained in 256 x 256 pixel images from the GONG
instruments. We make use of the Doppler velocity images of the full solar disk, magnetic field strength
images, total intensity and line depth/modulation images for a 1-Å band centered on the Ni I (6768 Å) line,
all at a 60 second cadence.

         During white-light flares, we observe changes in the localized region's intensity, modulation,
velocity, and magnetic field. Of particular interest are the magnetic fields, which are thought to twist up
tightly until they break and release energy in the form of a flare. While we expect the spectral line profile to
change during a flare, there is uncertainty as to whether the observed variations represent actual disturbances
in the velocity and magnetic fields or whether they are due in part to reading inaccuracies on the part of the
instrument. To investigate this question, we look at the GONG instrument readings of the X10 WLF flare of
Oct 29, 2003. Through a series of models, we generate an input spectrum that is 'observed' by a numerical
model of the Gong instrument. In order to determine the accuracy of the actual observations, we compare
the velocities and magnetic fields produced through the GONG model to the actual velocity and magnetic
field observed to asses the line shape effects. We will see how much v and b are affected by the spectral

52 NSO REU/RET Annual Program Report 2004                                    APPENDIX II – REU Student Reports
          We first produce a model of the Ni I 6768 spectral line that includes three parameters that affect the
shape of the line. The modeled spectral line we generate is based on photospheric emission during a strong
white-light flare, and hence the emission is not as great as it would be if it were based in the chromosphere.
Over the course of such a flare, many photospheric lines undergo reversals in their core—the lines go from
absorption to partial emission. In the time series image of the modeled spectrum, we see a clear shift towards
emission where the line becomes white. We prescribe two parameters for the emission core produced during
the flare thus accounting for its strength; one controls the emission core width and the other controls the
emission core amplitude. These values range from 0.1 to 0.55 and 0.2 to1.0, respectively, and increase in
increments of 0.05. The final parameter that contributes to the spectrum accounts for the motion of the entire
continuum by providing a shift around 1.0 in a range from 0.7 to 1.3. The changes are applied to the actual
line profile observed by the FTS at NSO/KP. Finally, we Doppler shift the entire spectrum to account for
velocity and generate left and right circularly polarized Zeeman shifted components.

         This modeled input spectral line with all the parameter variations is then passed through the model of
the GONG instrument. The Gong program provides modeled output intensities, modulations, and spectra that
are then compared to the actual data taken by GONG over the course of the WLF of Oct. 20. The
appropriate set of parameters and hence spectral lines are selected to match the observed changes in total
intensity and modulation. We use a model based on intensity and modulation values because these values are
better diagnostics for the spectral line than velocity and magnetic field values—they are more readily
observed with certainty. From there, we proceed to run the GONG model again with those corresponding
modeled spectra along with the velocity and magnetic field values given by the corresponding actual data
(since GONG runs on v and b). All this is processed as the input that is normally seen only by the
instrument. From this run-through, we are provided with output velocity and magnetic field values, which is
what we, as observers, would see as the readout data from GONG. Since our model provides us with insight
to the normally unknown input velocity and magnetic field realities in addition to the velocities and
magnetic fields traditionally 'observed' by GONG, we can compare the “before” and “after” values to
determine the accuracy of our instrument/whether what we observe is true or is mostly due to the instrument.

         From our results, we examine the active region before (frame 570), at the peak (583), and after the
flare (598). By comparing the actual observations to our modeled observations for each of these times, we
see that the velocity and magnetic fields appear to change little from image to image. In other words, the
measurements of these observables do not appear to be much affected by changes in the spectral line shape.
Magnetic field and velocity during a flare, then, are not dependent on the strength of photospheric emission.

         Examination of the histograms of the differences and ratios gives a different view of the results.
Here, we see that the velocity is actually altered by the line shape changes, but that the magnetic field is not
as much affected. The mean velocity difference is -86 m/s, and the mean magnetic field difference is -0.57 G.
The ratios of the intensity and modulation indicate that our grid of models needs to be improved—neither
distribution is peaked at a value of 1. The inadequate model grid may also be the cause of the additional
peaks in the ∆B distribution.

        We are faced with similar instrument uncertainty in the observations of surface oscillations. It has
been observed that the surface of the sun moves up and down, and that its velocity oscillates for periods of 5
minutes. Of particular interest here are the observed decreases in wave amplitudes in the presence of
magnetic fields. In our model, we generate a sin wave that undergoes a Fourier transform. We pass this
input through a series of similar programs and look for changes in the spectrum.

53 NSO REU/RET Annual Program Report 2004                                   APPENDIX II – REU Student Reports
         First, we model the effect of a non-flaring magnetic field region on the spectral line shape by
changing the “C-shape” line bisector. This simulates the suppression of the granulation and the convective
blue shift that is observed in active regions. Observational studies have shown that the bisector moves
towards longer wavelengths (redward) in the presence of magnetic fields, and that the curvature of the
bisector is reduced. To model this, we first obtain the actual bisector of the line H(λ) from the FTS spectrum.
We then alter the bisector according to the formula

                             H’(λ) = (1 - 0.9 (|B|/4000) ) H(λ) + 0.002 (|B|/4000)

where B is the magnetic field strength in G.

         To model the line shape effects on the oscillations, we generate a sine wave of velocity values with a
5-min period and a time span of 1 hour. We then generate a line profile for a given B, and send the model
profile, V time series, and input B through the GONG instrument model. We compute the FFT and power
spectrum of the output time series and compare it with the input spectrum.

        Ultimately we look to compare the input and output velocity time series to look for amplitude
changes. For a magnetic field of 1000G, while the original input time series has an amplitude of ~1 m/s and
a 0 m/s offset, the output sine wave has an amplitude of ~0.5 m/s and an offset of -65 m/s. The results thus
demonstrate a decrease in the oscillation amplitude, and an overall red shift which is in qualitative agreement
with oscillations. By plotting the modeled and observed decrease of oscillation power as a function of
magnetic field strength, we can conclude that at least some of the observed decrease is due to spectral line
shape changes (since they both demonstrate similar trends).

        We have modeled the effect of spectral line shape changes on the velocity and magnetic field
observed in flares, and on oscillation amplitudes in active regions. We find that the measured B is only
weakly affected by profile changes. Since there is little change in the observations of magnetic field as the
shape of the spectral line changes, the changes in the magnetic field that we observe over the course of a flare
are indeed due to changes in the Sun. However the observed velocity changes quite a bit with the spectral
line. The instrument readings of velocity are clearly more sensitive to the line shape. These results indicate
that we may need to change certain interpretations of observations. It is widely believed that oscillation
amplitudes in active regions decrease because of magnetic fields, but our findings imply that at least part of
that observed decrease in oscillation amplitude is due to profile changes. Thus, we must re-examine how we
interpret magnetic field effects on the oscillations. These results, however, are still in the early stages and
require further testing. We need to improve our models and may want to model other instrument in addition
to the GONG instrument. From this point, we plan to write a full paper and go further in depth into our

54 NSO REU/RET Annual Program Report 2004                                   APPENDIX II – REU Student Reports
                       NFS Research Experiences for Undergraduates
                                  Summer 2004 Report

                                           Heidi H. Gerhardt
                                           Towson University

                             National Solar Observatory/Sacramento Peak
                                  Advisor: Dr. K. Sankarasubramanian

                 Understanding the Mechanism of the Penumbra
                         and the Curled Light Bridge
          We analyzed the mechanisms of the penumbra and the curled light bridge by utilizing the
          SIR (Stokes Inversion based on Response functions) inversion program. The magnetic
          field strength, LOS (line of sight) velocity, and inclination were the main focus of both
          stratifications. The research on the penumbra was to compare the data with the siphon
          flow model and interchange convection model. From analyzing the data, pertaining to
          the penumbra, we concluded both models gave a logical representation of the penumbral
          mechanism. The light bridge system was under investigation for comparison to the
          cluster model. There is no definite answer, however with our data collection we
          concluded the cluster model gave a good theory on the mechanism of the light bridge.

1. Introduction
         The sunspot, NOAA 0484, was obtained at the Dunn Solar Telescope in Sunspot, New Mexico with
the DLSP (Diffraction-Limited Spectro-Polarimeter) on October 24, 2003. DLSP is a collaborative project
between NSO and HAO. The aim of this instrument is to produce high resolution spectro-polarimetric data.
The spatial resolution of this new instrument will be 4-times better than the well-known Advanced Stokes-
Polarimeter (ASP) (Sankarasubramanian et al., 2004).
         Two different aspects of NOAA 0484 will be discussed; the abnormal profiles in the penumbral
region and the curled light bridge. The main interest in understanding the mechanism in the penumbra was
the siphon flow and interchange convection models. We inverted the abnormal profiles using SIR to better
understand the magnetic field strength, LOS velocity, and inclination. The definition of abnormal profiles is
a deviation of the simple Zeeman Effect. These profiles are predominately located around the outer edge of
the penumbra (see Figure 1).
         The cluster model was utilized for understanding the mechanism of the light bridge. The observed
light bridge data was also run through the SIR inversion program. For this project we studied the magnetic
field strength, LOS velocity, and inclination. Now we will discuss the mechanism of the penumbra and how
the data correlates to the models.

55 NSO REU/RET Annual Program Report 2004                                 APPENDIX II – REU Student Reports
                                                                        Fig. 1: The white spots represent the
                                                                        abnormal profiles on the continuum.

2. Penumbra Study

2.1. Interchange Convection
       Schlichenmaier and his colleagues believe the interchange convection model is the mechanism of the
penumbra (2000). The theory is the surrounding quiet Sun heats the magnetic flux tube thereby making it
expand and less dense causing the tube to become buoyant. In order for the buoyancy forces to work, the
magnetopause of the background must exceed a critical inclination angle.

2.2. Siphon Flow
         The siphon flow model, introduced by Thomas et al. (2002) is similar to how pressure drops from
each end of a straw. Assuming the penumbra as a collection of flux-tubes (or straws of magnetic flux tubes),
the total pressure at the two end points has to be equal for the flux tube to remain stable. Observations show
that the flows in these flux tubes are from the inner penumbra to the outer penumbra. In order to obtain such
a flow, P2 must be less than P1. From the steady state condition, the above equality transform to B2 must be
greater than B1. Note that the assumption of constant radius of the tube is necessary and that the effect of
gravitation is neglected. In this case the magnetic pressure drops at the outer footpoint of the arch. In order
for this model to work the arch of the flux tube must have the following characteristics seen in figure 2.

                B12 + P1 = B22 + P2
                B2 > B1 and P2 < P1
Fig. 2: The left hand model is a representation of how the siphon flow model works in a magnetic flux tube.
        The direction of the flow is from left to right (image on the right, Thomas et al. 2002).

2.3. Procedure
         We studied the Stokes parameters (I, Q, U, and V) in two Fe I (6301.5 Å and 6302.5 Å) lines for the
SIR inversion program (see figure 3). Figure 3 provides an example of an abnormal profile, showing that
there are more than two lobes in the V parameter (one in the positive and one in the negative).

56 NSO REU/RET Annual Program Report 2004                                  APPENDIX II – REU Student Reports
         The principle of this procedure was to compare the magnetic field strength, LOS velocity, and
inclination against the mean height and compare the two models to our results.

Fig. 3: These plots are from the SIR inversion program, there were 1135 abnormal profiles. This is an
        example of abnormal profile 78. The orange line represents the input model and the black lines are
        the output models of the SIR inversion program.

So far we have found that the magnetic field is stronger at some points around the outer penumbral
region rather than in the inner penumbral region, as seen here in figures 4a, b.

                                         Field Strength in the inner penumbral region

   2450.00                                                                                             Fig. 4a: These points did
   2350.00                                                                                             not go through the SIR
   2300.00                                                                                             inversion; instead, they
                                                                                                       were from the M–E
   2150.00                                                                                             (Milne–Eddington)
                                                                                                       inversion. Each number
   2000.00                                                                                             in the x-axis corresponds
   1950.00                                                                                             to one pixel on the inner
             1   2   3   4   5   6   7   8   9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

57 NSO REU/RET Annual Program Report 2004                                                   APPENDIX II – REU Student Reports
Fig. 4b: These plots are the results from the SIR inversion of the abnormal profiles located in the
        upper right hand region of the penumbra. The top plot is the mean field strength in G over
        the height of formation, the middle plot is the mean LOS velocity in km/s over the height of
        formation, and the bottom plot is the mean inclination in degrees over the height of
        formation. Each number in the x-axis corresponds to one pixel on the outer penumbra.

These plots show the magnetic field strength in the inner penumbral region (see figure 4a) as well as in the
outer penumbral region (see figure 4b). The magnetic field strength is about 2,300 G in the inner penumbral
region and varies between 2,000 – 5,000 G in the outer penumbral region showing stronger downflows in the
outer footpoints. It means that some points at the outer penumbra show higher field strengths compared to
the inner penumbra. Those places (marked with arrows on figure 4b) also show larger downflows and those
field lines are more vertical (see figure 4b). These results may be taken as evidence for the siphon flow
mechanism. However, previous analysis by Sankarasubramanian, Rimmele, and B. Lites (2004) of the inner
penumbra of the same sunspot showed evidence for interchange convection mechanism [smaller field
strength, upflows, and critical angle of inclination at the inner penumbral footpoints].

2.4. Penumbra Conclusion
        Our results present both models to be sufficient in their explanation of the penumbral mechanism. It
is not definite whether one model is more predominate. It may be possible that both mechanisms operate
together. A detailed study of the effect of interchange convection mechanism for the outer penumbral region
is necessary, but has not been done to the knowledge of the author.

3. Curled Light Bridge Study

3.1. Cluster Model
         The cluster model consists of three characteristics: strong field strength at higher layers and a weaker
(or field free) field strength at lower layers, in the deeper layers the inclination is more horizontal than in the
higher layers, and the deeper layers should show convective flows. This theory was predominately for the
mechanism of umbral dots, however it has recently been considered for the mechanism of light bridges (see
figure 5a/b).

58 NSO REU/RET Annual Program Report 2004                                     APPENDIX II – REU Student Reports
Fig. 5a: The cluster model (Solanki, 2003).

3.2. Procedure
         After inverting the 96 profiles of the light bridge, we studied the mean magnetic field strength, LOS
velocity, and inclination at different optical depths to understand how the cluster model fits in with our data
(see figure 5a/b). Evidence of these characteristics is obvious.

Fig. 5b: These plots show field strength, inclination, and LOS velocity over the optical depth.

The field strength is stronger in the higher layers and weaker in the lower layers (around 2000 G), inclination
is more horizontal in the deeper layers, and the LOS velocity has upflows (2-4 km/s) in the deeper layers and
downflows (3 km/s) in the higher layers. We also looked at the individual points on the light bridge by
taking an average around each. In doing so, we found the field strength as low as 500 – 600 G seen in the
middle of the light bridge at the 48th profile (see figure 7). This is lower than seen in the previous procedure.

59 NSO REU/RET Annual Program Report 2004                                    APPENDIX II – REU Student Reports
Fig. 7: These are the graphs from the averaged points around 96 data points of the light bridge data.

In studying the physical properties more closely we found a few correlations between the magnetic field
strength and the continuum intensity and the inclination (see figure 6).

Fig. 6: The top plot shows an inverse correlation between the magnetic field strength (B) and the continuum
intensity (Ic). The equation defines when B decreases by 100 G, Ic increases by 3%. The bottom plot shows
when B decreases the inclination (ψ) becomes more horizontal, therefore the equation defines when there is
a change in about 100 G, the inclination changes about 3o. The symbol ‘V’ represents velocity.

60 NSO REU/RET Annual Program Report 2004                                 APPENDIX II – REU Student Reports
We then plotted graphs showing the relationship between field strength, inclination, continuum intensity, and
LOS velocity (see figure 8). The graphs show how the field strength decreases near the middle of the light
bridge and the LOS velocity has stronger downflows and the inclination becomes more horizontal.

Fig. 8: The solid line is field strength (top graph), solid line is LOS velocity (bottom graph), dashed line is
inclination, and the dotted line is continuum intensity.

3.3. Light Bridge Conclusion
         Our data confirms the earlier observational findings of Leka (1997) and Rimmele (1997). The
studied light bridge data follows the cluster models’ three main characteristic of having a strong field
strength at higher layers and a weaker (or field free) field strength at lower layers, in the deeper layers the
inclination is more horizontal than in the higher layers, and the deeper layers show convective flows.

        I would like to thank the NSF REU program at the National Solar Observatory for giving me this
wonderful learning experience and a special thank you to my advisor, Dr. K. Sankarasubramanian, for all his
help and intellectual influence.

       This work was supported by the National Science Foundation (NSF) and the Association of
Universities for Research in Astronomy. This work is carried out through the National Solar Observatory
Research Experiences for Undergraduate (REU) site program, which is co-funded by the Department of
Defense in partnership with the National Science Foundation REU Program.

61 NSO REU/RET Annual Program Report 2004                                    APPENDIX II – REU Student Reports

   1. Leka, K.D., The Vector Magnetic Fields and Thermodynamics of Sunspot Light Bridges: The Case
      for Field-Free Disruptions in Sunspots, The Astrophysical Journal, 484: 900-919, 1 August 1997.

   2. Rimmele, T.R., Evidence for Megnetoconvection in a Sunspot Light Bridge. The Astrophysical
      Journal, 490: 485-469, 20 November 1997.

   3. Sankarasubramanian, Rimmele, Lites, B., Diffraction Limited Spectro-Polarimetry at the Dunn Solar
      Telescope, American Astronomical Society 2004, 204.2006S, May 2004.

   4. Sankarasubramanian et al., The Diffraction Limited Spectro-Polarimeter: a new instrument for high-
      resolution solar polarimetry, SPIE 2004, 5171.207S.

   5. Schlichenmaier, R. et al., The Evolution of a Thin Magnetic Flux Tube Embedded in a Sunspot, 15
      February 2000,

   6. Solanki, S. K., Sunspots: An overview, The Astronomy and Astrophysical Review, 11:153 – 286, 3
      March 2003.

   7. Thomas, J. H., et al., Downward pumping of magnetic flux as the cause of filamentary structures in
      sunspot penumbra. Nature, 420: 390-393, 28 November 2002.

62 NSO REU/RET Annual Program Report 2004                               APPENDIX II – REU Student Reports
                        NSF Research Experiences for Undergraduates
                                   Summer 2004 Report

                                               Joel B. Lamb
                                             University of Iowa

                              National Solar Observatory/Sacramento Peak
                                        Advisor: Dr. Alexei Pevtsov

             A Study of Plasma Flows in Emerging Active Regions

                 We analyze plasma flows with respect to magnetic fields and visual phenomena in
          emerging active regions at the photosphere. We initially focused our study on two
          aspects of flows in emerging active regions; (1) flows that appear before any magnetic
          flux or white light indications of an active region and (2) the emergence of Evershed flow
          with respect to the emergence of the penumbra in white light. We found no consistent
          flows before emergence; however, after emergence downflows were observed in most
          regions. No conclusive results were obtained on the emergence of Evershed flow and
          penumbra. Additionally, we found in our data evidence of asymmetric flows which we
          will also cover in this paper.

         Plasma flows in emerging active regions are not well understood. In the case of flows before
emergence, very little research has been done to the knowledge of this author. However, one theoretical
model by Lorrain & Koutchmy (1998) proposes that downflowing plasma is the mechanism by which
sunspots are self-excited magnetic-flux-tube dynamos. According to this model, one could expect to find
downflows immediately preceding the appearance of a sunspot in an otherwise quiet region of the sun. This
model is our impetus for looking at flows prior to active region emergence. After emergence, many authors
have reported downflows in the area around developing pores and sunspots, see, for example, Keil,
Balasubramaniam, Smaldone, & Reger (1999).
         Evershed flows are a well-studied yet poorly understood aspect of developed sunspots. Since its
discovery, (Evershed, 1909) the literature has been filled with observations and studies of the phenomenon.
It is widely accepted that the onset of Evershed flow occurs simultaneously with the appearance of the
sunspot penumbra. For examples, see Rimmele (1995) and Leka & Skumanich (1998). In our comparison
between the emergence of the Evershed effect and penumbra, we attempt to determine a causal relationship
between the appearance of the penumbra and the beginning of Evershed flow.
         Little research has been done on asymmetric flows in emerging active regions. An asymmetric flow
in emerging active regions was first predicted in a series of three-dimensional numerical simulations of the
dynamical evolution of emerging flux loops done by Fan, Fisher, & DeLuca (1993). This simulation was
successful in describing a number of properties observed in active regions. The asymmetric flow predicted
by this simulation is a flow from the preceding polarity or sunspot to the following polarity. In an attempt to
confirm these predicted asymmetric flows, Cauzzi, Canfield, & Fisher (1996) embark on a search for these
flows in emerging active regions. While they were successful in finding asymmetric flows, all of the flows
were in the opposite direction (following to preceding) of that predicted by the simulation. Although not an
intended outcome of our research, the presence of a few asymmetric flows in our data are worthy of
consideration in this study.

63 NSO REU/RET Annual Program Report 2004                                   APPENDIX II – REU Student Reports
Data Processing
         The data for this study came from the Michelson Doppler Imager (MDI) instrument aboard the Solar
and Heliospheric Observatory (SOHO) satellite. Dopplergrams, magnetograms, and full continuum images
are used. The images are 1024 by 1024 pixel images centered on the Ni I photospheric absorption line at
6767.8 Angstroms. The Dopplergrams are created using wavelength shifts in the Ni I absorption line to
obtain line of sight velocities. The Magnetograms are taken at a cadence of 96 minutes (15 per day) while
both the Dopplergrams and the full continuum images are taken at a cadence of 6 hours (4 per day). The
images have a resolution of 4 arcsec and a plate scale of 2 arcsec per pixel.
         We first used the magnetograms to calculate the center of gravity of magnetic flux for the emerging
active region over 3 to 6 days worth of data. Graphing the coordinates of the center of gravity of magnetic
flux vs. the time that the image was taken and applying a least squares linear fit to the data gave us a linear
function whose slope is the rotation rate of the active region. Extrapolating this rotation rate backwards in
time allowed us to calculate the coordinates of the active region on the sun at any time before it emerged.
Using this rotation rate, we studied the emerging region in the dopplergrams from up to two days before
emergence until up to four days after emergence. In order to gather uniform data from the dopplergrams,
two corrections had to be made to the images. The first was a correction for the line of sight rotational
velocity of the sun for each pixel in the image. The second correction was an overall calibration for each
image. A horizontal line of pixels, taken approximately at the equator, was created and the line of sight
velocity values for each pixel along this line were plotted in sequence. Applying a fitted line to this linear
function provided a velocity value at a central meridian distance (cmd) of -90 degrees and 90 degrees. Since
these two values should be of equal magnitude and opposite sign, simply adding the two values and dividing
by two gives the correction.
         With the corrected and calibrated dopplergrams, the rotational rate linear fit was used to determine
the center of the active region on each image. Around this point, an 11 by 11 pixel box was constructed in
which the average velocity and standard deviation were calculated from these points. This 11 by 11 pixel
box was also created in the magnetograms where the sum of the magnetic flux, the average magnetic flux
and the standard deviation were all recorded. Finally, the full continuum images also had the 11 by 11 pixel
box constructed in them where the minimum intensity, average intensity, and standard deviation were all
found. All of these values were plotted vs. time to check for temporal trends in these values for emerging
active regions.
         For the research concerning the Evershed effect, the same rotational rate linear fit was used to
determine the center of the emerging active region. Using this point as a center, a box of variable size
corresponding to the size of the active region was created. This box was split vertically into two evenly sized
smaller boxes so each polarity sunspot could be studied separately. Inside of each smaller box, the minimum
intensity value was found and the heliographic coordinates of that point were noted. Using these
heliographic coordinates as a center point, a line 51 pixels in length was constructed that, if extended, would
pass through the center of the solar disk. Both the velocity values and the full continuum intensity values
were saved for every pixel along this line. These values were plotted and compared to determine the
existence of both the Evershed flow and penumbra.
         To determine the existence of asymmetric flows, a combination of viewing the corrected/calibrated
dopplergrams plus examination of the line profiles created for detecting the Evershed effect were used.
Since these methods are a sort of ad hoc approach, only regions with a very noticeable and distinct
asymmetric flow are discussed in this report.

         The correlation or predictable changes through time in the 11 by 11 pixel box parameters were either
weak or nonexistent. Obviously, the minimum continuum intensity value decreases with the formation of
pores and sunspots. By the same token, standard deviation of continuum intensity increases with the
appearance of pores and sunspots. The average intensity value decreases with the emergence of pores and
sunspots. However, this effect is overshadowed by the limb-darkening effect which will increase the average
intensity as the active region moves toward disk center or decrease the average intensity as the active region

64 NSO REU/RET Annual Program Report 2004                                  APPENDIX II – REU Student Reports
moves away from disk center. The average values for magnetic flux show very little correlation temporally
from one active region to the next. Only the standard deviation of magnetic flux shows a consistent increase
in time for each active region, which of course is expected. The dopplergrams show no steady flows before
emergence of the active region. Only right at emergence do most active regions exhibit a strong downflow.
The change in average velocity is within the standard deviation of the velocities; however, the number of
regions that exhibit this downflow at emergence are enough to validate a real difference in plasma flow
before and after emergence. A general downflow trend persists in the majority of the active regions as they
continue to emerge. Refer to figure 1 for average and minimum intensity vs. time plot of a typical average
region. Figure 2 shows the average magnetic flux vs. time and figure 3 shows the average velocity vs. time.
All error bars are one standard deviation for the pixel values inside the box.

         Fig 1: Intensity Parameters vs. Time                  Fig 2: Average Magnetic Flux vs. Time

Line: Minimum Intensity
Points: Average Intensity

                                                                          Fig. 3: Average Velocity vs. Time

65 NSO REU/RET Annual Program Report 2004                                 APPENDIX II – REU Student Reports
        The line profiles of dopplergrams and white light images provide some very interesting results
concerning the formation of Evershed flows and the penumbra. A typical line profile of the Evershed effect
on continuum and dopplergram images is shown in figures 7 and 8 respectively with the continuum image of
this sunspot shown in figure 4. Notice the distinguished change in the line profile between the umbra and
penumbra in figure 7. In figure 8, the distinctive Evershed effect is seen on the line profile as first an upflow
(blueshift) which is traditionally a negative value followed by a downflow (redshift) represented as a positive
value. One emerging active region in particular, AR 8861, shows very trace elements of a penumbra yet
displays a rather robust Evershed flow. Figure 5 show the continuum image of AR 8861 while figures 9 and
10 show the line profiles of continuum and dopplergram images respectively. Finally, one active region, AR
8737, displays a full formed sunspot with a large penumbra while showing no sign of any Evershed flow.
Figure 6 shows the continuum image of AR 8737 with figures 11 and 12 showing the continuum and
dopplergram line profiles.

          Fig. 4: AR 8200                   Fig. 5: AR 8861                   Fig. 6: AR 8737
          Sunspot on right                  Sunspot on right                  Sunspot on left

       Fig. 7: AR 8200 intensity line profile                       Fig. 8: AR 8200 velocity line profile

66 NSO REU/RET Annual Program Report 2004                                    APPENDIX II – REU Student Reports
         Fig. 9: AR 8861 intensity line profile                     Fig. 10: AR 8861 velocity line profile

        Fig. 11: AR 8737 intensity line profile                     Fig. 12: AR 8737 velocity line profile

         Three emerging active regions in this study showed asymmetric plasma flows during emergence.
Each of these flows appear somewhat short in duration, the longest lasting for about one day. The first active
region, AR 8699, shows a striking asymmetric flow from the preceding sunspot to the following sunspot. To
this author's knowledge, this is the first observation of an asymmetric flow from the preceding to following
spot in an emerging active region. This data provides the first observational confirmation of this type of
asymmetric flow predicted by the three dimensional numerical simulations of Fan, Fisher & DeLuca (1993).

67 NSO REU/RET Annual Program Report 2004                                  APPENDIX II – REU Student Reports
The continuum and dopplergram images of AR 8699 are shown in figures 13 and 14 respectively. The
dopplergram shows black as upflow and white as downflow. The second active region, AR 8617, shows an
asymmetric flow in the opposite direction, from the following sunspot to the leading sunspot. This flow is
similar to the flows found in Cauzzi, Canfield, & Fisher (1996). Figures 15 and 16 show the continuum and
dopplergram images of this flow. Even more striking perhaps are the line profiles of the dopplergrams for
each polarity of sunspot, which show very large upflows and downflows in the following and leading spots
respectively. The following polarity line profile is shown in figure 17 while the leading polarity line profile
is shown in figure 18. The final active region with asymmetric flows, AR 9244, again shows the following
to leading flow path. The asymmetric flow in this active region is particularly interesting because of its
unique evolution in time. At first, the asymmetric flow is concentrated on the opposite polarities of magnetic
flux. The next images taken 6 hours later show that the asymmetric flow is now concentrated at the sunspots
of each polarity instead. The first set of continuum, magnetogram, and dopplergram images for AR 9244 are
shown in figure 19 while the second set are shown in figure 20.

     Fig. 13: AR 8699            Fig. 14: AR 8699           Fig. 15: AR 8617       Fig. 16: AR 8617
     Continuum image             Dopplergram image          Continuum image        Dopplergram image

       Fig. 17: AR 8617 following polarity                           Fig. 18: AR 8617 preceding polarity
             Velocity line profile                                            Velocity line profile

68 NSO REU/RET Annual Program Report 2004                                   APPENDIX II – REU Student Reports
          Fig. 19: AR 9244 continuum, magnetogram, and dopplergram
          Notice strong correlation between magnetogram and dopplergram

          Fig. 20: AR 9244 continuum, magnetogram, and dopplergram
          6 hours later. Notice correlation between continuum and dopplergram.

Conclusions and Discussion
         We have studied plasma flows in emerging active regions, both before and after emergence. Our
data analysis has focused on consistent flows before emergence, Evershed/penumbral appearance, and
asymmetry in plasma flows after emergence.
         We have not found any consistent flows before active region emergence. This appears to be in
conflict with the idea of a sunspot as a self-excited dynamo proposed by Lorrain & Koutchmy (1998). After
emergence, an overall downflow occurs in the area of the active region. This downflow is sustained as the
active region continues to emerge and mature.
         Our data concerning the Evershed effect is mixed as we have found examples of Evershed with little
or no penumbra while also finding quite clearly a penumbra without any sign of Evershed flow. The data is
wholly insufficient to determine a causal relationship between the Evershed effect and the appearance of a
penumbra. Clearly, more research at higher resolution, with a larger dataset, and with greater temporal
resolution needs to be undertaken for any conclusions to be drawn on this subject.
         We have found conflicting examples of asymmetric flows in emerging active regions. One active
region provides observational confirmation of the predicted asymmetric flow from the preceding polarity to
the following polarity described by Fan, Fisher, & DeLuca (1993). Two other active regions display an
opposite asymmetric flow in conflict with Fan, Fisher, & DeLuca (1993). This following to preceding flow
agrees with the asymmetric flows found by Cauzzi, Canfield, & Fisher (1996). It appears more research
needs to be done on asymmetric flows both theoretically and experimentally before definitive conclusions
can be drawn.

        This work was supported by the National Science Foundation (NSF) and the Association of
Universities for Research in Astronomy. This work is carried out through the National Solar Observatory
Research Experiences for Undergraduate (REU) site program, which is co-funded by the Department of
Defense in partnership with the National Science Foundation REU Program.
        I am immensely grateful for the time and guidance so graciously given by my advisor Dr. Alexei
Pevtsov. His mentoring helped bring this project to fruition. I would like to thank Tim “Tex” Henry for
sharing his programming expertise whenever a program hit a snag. I also would like to thank the other
summer students for their help and support throughout the project. A big thank you goes out to the entire
Sunspot community for a fabulous summer.

69 NSO REU/RET Annual Program Report 2004                                 APPENDIX II – REU Student Reports
Cauzzi, G., Canfield, R. C., & Fisher, G. H. 1996, ApJ, 456, 850
Evershed, J. 1909, MNRAS, 69, 454
Fan, Y., Fisher, G. H., & DeLuca, E. E. 1993, ApJ, 405, 390
Keil, S. L., Balasubramaniam, K. S., Smaldone, L. A., & Reger, B 1999, ApJ, 510, 422
Leka, K. D., & Skumanich, A. 1998, ApJ, 507, 454
Lorrain, P., & Koutchmy, O. 1998, A&A, 339, 610
Rimmele, T. R. 1995, A&A, 298, 260t

70 NSO REU/RET Annual Program Report 2004                               APPENDIX II – REU Student Reports
                      NSF Research Experiences for Undergraduates
                                 Summer 2004 Report

                                         Statia H. Luszcz
                                         Cornell University
                                National Solar Observatory/Tucson
                                        Advisor: Matt Penn

            Magnetic Field Measurements and
          Moving Magnetic Features (MMFs) Using
             the Infrared Fe Line at 1564.8 nm

        Over 30 years of observations and theories comprise the study of moving magnetic features
(MMFs). MMFs, whichoften occur in opposite polarity pairs, migrate radially outwards through
the sunspot moat at speeds of about 1 km/s. Magnetogram images and maps of the magnetic field
vector using the IR Fe line at 1564.8 nm have penumbral magnetic features that appear to move
radially outwards at similar velocities and azimuth angles as those of MMFs in the moat, some-
thing not predicted by preceding work. We use polar time slice images to measure the radial
velocities of both types of features in an attempt to determine the validity of this relationship, as
well as to evaluate its implications for our understanding of sunspot magnetic fields.

MMFs and BPs: Previous work
        When Hale discovered magnetic fields in spots, he effectively introduced magneto-
hydrodynamics into solar physics. The study of solar magnetic fields has since become an
important component of solar physics. Sunspots are the best observable features in the
intermittent solar magnetic field. A typical sunspot consists of an umbra that is dark in
continuum, cooled by its strong vertical magnetic field. The umbra is surrounded by a brighter
penumbra, where the magnetic field flares out, reaching a near-horizontal inclination. During its
decay, the penumbra is often surrounded by a moat, an annular zone extending to about 20 Mm
and lacking a stationary magnetic field. (Zwaan)
        Magnetic flux in a sunspot is contained in narrow filaments called flux tubes. These flux
tubes, while below the current limit of telescope resolution, are the key to the fine structure
dynamics of the sunspot magnetic field. Therefore, our ideas of small magnetic features, such as
MMFs and bright points, are as much theory as observation, and continue to be studied with
increasing interest.
        In 1971 Vrabec noted the systematic outward streaming of small (<1’) magnetic knots of
both polarities within annular areas surrounding several sunspots. Harvey and Harvey (1973)
termed these knots ‘moving magnetic features’ (MMFs) in what was the first real, systematic study
of these phenomena. MMFs aren't visible in continuum movies, so H73 studied the MMFs using
time series of magnetograms and H-alpha filtergrams. They found that MMFs appear only around

71 NSO REU/RET Annual Program Report 2004                            APPENDIX II – REU Student Reports
sunspots which possess a moat. Like their predecessors, Sheeley and Vrabec, H73 found that
MMF's first appear at the outer penumbral boundary and travel in near-radial directions at speeds
of ~1 km/s, until they vanish or reach the network. They estimate the flux of individual features at
6e17 to 8e19 Mx. While MMFs tend to appear in opposite-polarity pairs, and may carry 2-10 times
the flux of the parent sunspot, the rate net flux carried away from the sunspot in MMFs is of the
same magnitude of the decay rate of sunspot flux.
         While moat MMFs have repeatedly observed, reports on temporal variations in the sunspot
penumbra are far less common. Bovelet and Weihr looked at flow in both the moat and penumbra
by measuring the motions of G-band bright points. In 1969 Sheeley noted in his spectroheliograms
BP’s moving outwards from sunspots at speeds of ~1km/s. Some of these BPs were cospatial with
the magnetic field concentrations. Bovelet and Weihr studied the moat flow by looking at a large
population of BPs that were cospatial with magnetic field concentrations. To track the BPs over
their lifetime, they used an algorithm called multiple level tracking. MLT was capable of tracking
1426 isolated BPs by separating and extending each recognized feature into a realistic shape, and
using that shape to determine flux and exact position from the moment each feature reached an
intensity threshold, until the moment it sank back below that threshold. They were then able to
determine the magnitude and angular offset from radial of the velocity of each BP. Their
measurements had a large random component, but showed an overall systematic outward radial
drift of .3 km/s. They found these speeds to be at a maximum at the outer penumbral boundary,
decreasing to nearly zero at ~20 Mm, which they took to be the extent of the moat.
         Bovelet and Weihr also studied ~275 bright structures within the penumbra. Unlike the
BPs in the moat, the Penumbral Bright Structures possess both inward and outward velocities,
with the outward velocities concentrated in the outer penumbra. They also studied Penumbral
Dark Structures, which are located primarily in the outer penumbra and move predominantly
outwards. They believe that some of the moat ‘drift’ might persist in the outer area of the
penumbra, causing this outward motion.
         Another study that observed inward and outward motions in the penumbra was done by
Lites et al in 1997. They found that at low frequencies in maps of fluctuations in continuum
intensity, magnetic field strength, and field inclination, distinct features moving radially inward
can be seen in the inner penumbra. In the outer penumbra and surrounding quiet sun, distinct
outward–moving features are visible. They believed these features to be signatures of the
convective interchange of magnetic flux tubes in the sunspot.

         The data I worked with are spectropolarimetry scans taken in June of 2002 at the McMath-
Pierce telescope on Kitt Peak. The camera used was the CSUN-NSO infrared camera, using a
256 x 256 pixel detector. The observers made polarized, east-west scans of the spectrograph over
the active region NOAA 10008, a sunspot located just southeast of disk center. Sixteen scans were
taken from 17:38 to 21:59 UT at a 15 minute cadence. Two of these scans were eventually thrown
out because they were corrupted. The magnetic field measurements were made using the
instrument -corrected Fe line 1564.8 nm. As far as we could find, nobody has yet studied sunspot
magnetic fields using this line: however, it is ideal for studying moving magnetic features for a
variety of reasons.
         The data were analyzed by measuring the Zeeman splitting of the line, which is directly
proportional to the strength of the field causing the splitting, allowing us to directly measure that
field strength . The pattern and amount of splitting is also a property of the magnetic the spectral
line used, namely the wavelength and a property of the transition known as the g-factor. Our IR

72 NSO REU/RET Annual Program Report 2004                             APPENDIX II – REU Student Reports
Fe I line at 1564.8 nm has a longer wavelength than any optical choice, as well as a very high g
factor, making it an excellent choice for magnetic field study.
        Besides our choice of line, our data has the benefit of a long data-collecting period of 4.5
hours of data. The slow .2-2 km/s speeds of the type of features we are looking at require a long
period of observation for measurable results. These features also have a lifetime of ~4 hours, so an
observational period of this length is preferable. The combined benefits of the line choice, full
vector measurements, and long observation period make this data uniquely good for studying the
sunspot magnetic fields.
        Despite these attributes, our data was difficult to analyze because of seeing. MMFs, on the
scale of only ~1’, require high resolution. However, in our data, atmospheric turbulence causes
image motion and degradation that severely limit resolution. Turbulence in the atmosphere is
compounded by turbulence that builds up inside the telescope tube as it is heated by the sun. This
problem can be minimized by evacuation of the tube or by adaptive optics. Other telescopes, like
the DOT on the Canary Islands, take advantage of the location's steady, cool breeze to remove
turbulence between the mirror and detection apparatus. Our data were taken without benefit of
any of these devices. Additionally, the slow speed of scanning means seeing can change
drastically over the course of a single frame.

         Carl Henney created maps of continuum intensity, magnetic field strength, mf inclination
to the line of sight, and azimuth angle of the transverse field by applying a least-squares fit to the
magnetic field components at every pixel. The line of sight magnetic field strength map shows a
double umbra. The northward umbra is darker in continuum, about .45 times quiet sun intensity,
and has a stronger magnetic field ~2800 G, than the southern umbra (I=.55 qsun, field=2300G).
Also visible is the penumbra, cooler than quiet sun, but with a weaker magnetic field (1200-1400G).
In inclination, the umbra is dark, as it's inclination to vertical is near zero. In contrast, the
penumbra shows an inclination of about 90 degrees.
         These inclination maps also reveal another interesting phenomenon: ‘wiggles’, or regions of
less inclination to vertical. In contrast to the near-horizontal inclination of the rest of the
penumbra, these regions possess field inclinations of 40-60 degrees to the line of sight. Time
sequences show that these features are dynamic on the period of observation, moving radially
outward through the penumbra. Our goal was to look at the velocities and azimuth locations of
these features, and explore the possible relationship between these features and the MMFs of the
         The line fitting method used to make the magnetic field maps does not work well on the
weaker field outside the spot. Therefore, a different and simpler method was used to study
motions in the moat area: we created a time series of magnetograms by differencing the line wings
of the Fe 1564.8 nm Stokes V spectra, effectively integrating over that line. After shifting the
images so that the spot center was aligned between frames, and removing two bad frames, we
were left with a time series of 14 magnetograms spaced over the 4.5 hour period.
         Our goal with the magnetogram images was to find the velocities of the radially moving
MMFs, with hopes of later comparing these velocities to those of the penumbral inclination
‘wiggles’. Measuring the radial velocities was made difficult by atmospheric seeing, which often
blurred some of the features in a frame. Therefore, automated methods of feature tracking, like
MLT developed by Bovelet and Weihr, were not applicable.
         My first attempt to identify and track MMFs was to visually and manually track moat
features over the course of the 14 good frames, and use the positional changes to calculate the
radial velocities of the MMFs. This revealed an average outward velocity of ~.25 km/s, with

73 NSO REU/RET Annual Program Report 2004                            APPENDIX II – REU Student Reports
directions varying only to about 20 degrees of radial. However, the paths of these features could
not be visualized in time slice diagrams produced along the direction of motion of each feature.
         Next we attempted to use a technique of local correlation tracking, similar to the method
November and Simon used for proper-motion measurements of solar granulation. Local
correlation tracking corrects for atmospheric geometric distortions in a time series by cross
correlating small segments of the images in time. We attempted a two-fold purpose with local
correlation tracking: we hoped to remove large-scale atmospheric shifts, as well as uncover the
organized radial shift due to moat flow. To do this, we adapted our software to output shift
values, and plotted them to a vector field. The result was noisy and lacked any organization: we
realized that the program, when lacking a real correlation, would output a false and often large
shift value. To improve our correlations, we tried adapted versions of several other correlation
programs, none of which provided a clear visualization of the moat flow. We attempted to again
adapt the software to use only well-correlated shifts, once again resulting in a disorganized shift
field. Our ultimate conclusion was that while local correlation works well for images of the quality
of November and Simon, our significantly larger noise levels made use of LCT impossible.
        Our final attempt at MMF identification and tracking was to revert to using time slice
diagrams. However, this time, we approached it by imaging the sunspot in polar coordinates. The
goal of time slices is to image in one spatial dimension and a time dimension. In polar coordinates,
radial motion becomes spatially one dimensional. By creating time slices over the range of
azimuthal angles, I identified nearly 40 MMFs in the magnetograms. We measured the motion of
each feature in space and time to find the radial velocity (dr/ dt) and then map these velocities to
the the spots. We also repeated the time slice procedure to look for moving features in the maps of
magnetic field inclination. This method revealed an additional 10 moving features.

        The MMFs in the magnetograms had a variety of velocities, from .18 to 2.51 km/s. These
results were consistent with previous work, including Ryutova’s observations that MMFs possess
a wide velocity range. The median velocity was about .4 km/s. As Zhang and Solanki found, we
also noticed that the MMFs appear to cluster at particular azimuths around the sunspot. The
majority of the features were members of opposite polarity pairs.
        The most unexpected result from the magnetogram studies was that several of the MMFs
identified moved within the penumbra. Nearly half of the 38 magnetic features from the time
slices were first discovered moving outwards through the outer penumbra, lying in contradiction
to the strict outer-penumbral starting boundary previously supposed for MMFs. These features
appear to move radially outwards in lanes between more vertical field concentrations. Their
velocities remain relatively constant as they pass through the outer penumbral boundary into the
moat, and there is no azimuthal or velocity distinction between the penumbral-start MMFs and the
moat-start MMFs. Additionally, the 10 features we tracked in the inclination maps of the magnetic
field move radially outwards at similar velocities to the MMFs, in the range of .18-1.7 km/s. These
features appeared to cluster at the same azimuth angles as the MMFs, suggesting that the types of
features are related. Our discovery of these two new features in the penumbra has an important
impact on our understanding of the structure and formation of these features.

        Harvey and Harvey first proposed that MMFs are the mechanism of the slow decay of
sunspots, and this idea has been adopted by a variety of later works. While most sunspots decay
fast by fragmentation, many larger spots—the ones surrounded by moat cells—decay gradually.

74 NSO REU/RET Annual Program Report 2004                           APPENDIX II – REU Student Reports
Across the moats of such spots stream MMFs of both polarities, with the total rate of flux removal
in MMFs indicative of the the decay rate of the sunspot.
        Because of the significant role of MMFs in sunspot dynamics, it is important to try to
understand the MMF in the context of the sunspot environment. Features on the scale of single
flux tubes are beyond our resolution capabilities, so we have been forced to settle for models of
MMFs, which take into account observations and what we know about the magnetic and
convective properties of a sunspot.
        Three types of MMFs have been observed around sunspots. According to Thomas et al,
Type I MMFs, the most commonly observed, are opposite-polarity pairs. Type II MMFs are single
magnetic elements with the same polarity of the sunspot, moving outward across the moat with
speeds similar to those of Type I MMFs. Type III MMFs are single magnetic elements with polarity
opposite that of the sunspot. These MMFs move rapidly outward through the moat. Types I and
III MMFs are only seen around sunspots with penumbrae, while a penumbra is not a requisite for
Type II MMFs. Most models tend to either omit talk of Type II and III MMFs, or attribute them to
opposite-polarity pairs in which one member is not visible. The main focus of MMF models is on
Type I MMFs. While most theories conclude that MMF pairs are a the footpoints of a small part of
a magnetic field line that loops through the solar surface, the cause and nature of these U loops
differs greatly between models.
        Harvey and Harvey first proposed that a series of MMF's is actually the repeated
emergence of a sub-photospheric sunspot field line. Ryutova adopted this intuitive picture, citing
the wide range of flow velocities observed in ‘near sunspot’ regions as the cause of flux tube
deformation and the generation of a stable kink. This ‘kink’ may then evolve into a traveling MMF
or series of MMFs. Each flux tube has its own range of critical flow velocities which may cause it
to deform. The exchange of energy and momentum between magnetic flux and outer motions
result in the onset of nonlinear shear stabilities which lead to the observed macroscopic effects on
the flux tubes.
        Ryutova also actually claims to see the flux tube emerging, its ends becoming a pair of
opposite-polarity MMFs. In accordance with these observations, Spruit et al present a second
model of rising loops. The properties of active regions produce U-loops: flux tubes having two
ends at the photosphere but otherwise still embedded in the convection zone. Once a U loop
develops, the bottom of the flux tube becomes magnetically buoyant and rises. As it rises, the flux
tube expands, and field strength drops, becoming controlled by convection on top of the
convection zone. At the surface, the large loop is fragmented into a series of closed loops, the
footpoints of which are the members of an opposite polarity MMF pair. Flow in the moat cell
grabs the top part of the rising loops and carries them radially outward.
        Unlike Spruit’s the rising U loop model, Zhang and Solanki propose that MMFs are formed
when the field lines in a small part of the magnetic canopy dip down to produce a ‘falling’ U-loop.
They begin their scenario with a packet of dense, outward-flowing Evershed gas in the penumbra.
The vertical field gradient of the penumbra provides a force that resists the gravitational force. At
the outer penumbral edge, this supporting force disappears, and sufficiently massive and dense
gas can no longer be supported by the flux tube field. Therefore, the gas sinks, taking the magnetic
field with it, causing a U loop to be created near the penumbral edge, with gravity and magnetic
tension working to submerge and push up the flux tube, respectively.
        Thomas and Weiss explain the production of MMFs through the interaction of the magnetic
field with turbulent, compressible convection in the solar granules and supergranules. In their
model, granules grab flux tubes, with their associated outflows of gas, from the sunspot penumbra
and pull them under the surface of the sun. They are then carried outwards by a large

75 NSO REU/RET Annual Program Report 2004                            APPENDIX II – REU Student Reports
supergranule, ~30 Mm in diameter that is centered on and stabilized by the sunspot. Occasionally,
a submerged flux tube surfaces, producing outward-moving magnetic loops.
       Unlike other theorists, Thomas et al do not dismiss single polarity MMFs. Type I MMFs,
along with Type III MMFs rely on the near- horizontal flux tubes of the penumbra, and correspond
to deeply submerged flux tubes. Type II MMFs correspond to vertical flux tubes stripped off of the
sunspot umbra and carried outwards.

        These three models all predict MMFs, predominantly in opposite–polarity pairs.
However, none of the theories or observations thus far has observed an MMF manifestation in the
moat. As Shang and Solanki said, the MMFs tend to first appear at a distance of 1000 – 5000 km
from the outer boundary of the sunspot. However, we note MMFs in the penumbra, as well as
magnetic inclination features in the penumbra that seem to be correlated with MMFs.
       If we look at our penumbral features, it is most likely they are not related to the PBS or PDS
observed by Bovelet and Weihr, as the speeds and predominant directions of these features differ
from our observations. The velocities and azimuth angles of our penumbral features are more
comparable to the velocities and clustering azimuth angles of the MMFs we measured in the moat.
       Most of the theories thus far either don’t address or completely disregard a possible
penumbral MMF manifestation, probably because no such manifestation has yet been reported.
These theories also depend, to varying degrees, on the environmental change between the
penumbra and moat, and this dependence allows us to eliminate some theories based on our
observations. For example, the justification of Zhang and Solanki that the lack of magnetic field in
the moat is what causes the MMFs to form, lies in direct contradiction to our observations of the
moving penumbral magnetic features. More analysis must be done to relate these features to the
other models. However, it is clear that the discovery of penumbral MMFs forces us to reevaluate
how and where MMFs are formed. Further study is necessary to confirm our results and to
provide more evidence towards the relationship between the penumbral and moat features.
Accordingly, we have repeated the observations of June 2002 this summer, using a better camera
and the McMath’s adaptive optics bench. Hopefully, analysis of this new data will help solve the
mystery of the dynamics of MMFs.

References/ Summer Reading:

Harvey, J., & Harvey, K. 1973, Sol. Phys., 28, 61
Vrabec, 1971
Sheeley, 1969
Simon & November, 1988, Precise proper motion measurements of solar granulation, 1988 ApJ 333
Zhang, J. et al 2003, A&A 399, 755-761
Spruit, H.C., et al 1987, Sol. Phys., 110, 115
Lites, B.W. et al 1998 ApJ 497: 464-482
Ryutova, M. et al 1997 ApJ 492:402 414
Bovelet and Wiehr 2003, A&A 412, p249
Weiss et al. "...Downward Pumping of Magnetic Flux” 2004, ApJ 600, p1073
McPherson, Lin, Kuhn, Solar Physics 1992, v139, p255 "Infrared Array Measurements of
   Sunspot Magnetic Fields"
Solanki, S. ‘Sunspots: an overview’, 2003, A&A, 11: 153-286
Foukal, Solar Astrophysics ‘Sunspots’

76 NSO REU/RET Annual Program Report 2004                            APPENDIX II – REU Student Reports
                         NSF Research Experiences for Undergraduates
                                    Summer 2004 Report

                                           Michelle T. McMillan
                                         Northern Arizona University
                               National Solar Observatory/Sacramento Peak
                      Advisors: Dr. K. Sankarasubramanian and Dr. H. Uitenbroek

                 On the Accuracy of Polarimetric Inversion Codes
                         for Spectroscopic Data Analysis
      Inversion codes are a common way to retrieve a physical solar model from spectroscopy. In this
      paper we test several inversion codes and initialization methods against theoretical two
      dimensional simulations of spectral formation to determine their accuracy. We found that with
      the correct initialization relatively good values, with correlation coefficients of greater that 50%,
      can be found over the height of formation of the profiles that are being inverted. Also, we found
      that codes that output an average value for the physical properties such as the URI inversion code
      are much faster than codes that return a gradient model. The former can be used for initialization
      of the latter to obtain the most accurate model possible.

1. Introduction
         In order to understand the processes that control solar activity, we must understand what magnetic fields
exist and how they interact in the atmosphere to cause the variety of process that are observed. The most accurate
way to find the values of the magnetic fields is through the Stokes polarization parameters. These parameters
quantify intensity as well as linear and circular polarization. Although other methods exist to extract the field
strength, the Stokes parameters can provide a vector map of the field including strength, angle of inclination, and
         It is important to have an accurate method to extract this vector map from the observed profiles for
several reasons. Besides satisfying our curiosity about the fields that exist in the sun, an accurate method to
determine magnetic field values can be used to compare with model calculations in order to understand the
processes that create those fields. In this paper we will test the accuracy of inversion codes to infer magnetic
fields as well as other physical parameters. We do this by comparing the model that was inverted from each
profile with the model that was used to calculate that profile in the first place. This assessment of accuracy
becomes possible only now that realistic simulations of the solar atmosphere exist.
         Inversion codes, from simple Milne-Eddington (M-E) codes to more complex models, are used to give an
idea of how certain parameters effect each of the Stokes profiles and in turn what parameters must exist to
produce the observed profiles. However, problems such as the risk of the χ2 merit function settling into a
secondary minimum instead of finding the best possible solution exist. This risk can be overcome by initializing
the code correctly.
         We tested four inversion codes, two M-E codes URI (modified Unno-Rachkousky Inversion code) and
MELANIE (Milne-Eddington Line Analysis using a Numerical Inversion Engine) as well as two other codes
which retrieve gradients, SIR (Stokes Inversion based on Response functions) and LILIA (LTE Inversion based on
the Lorien Iterative Algorithm). The URI and MELANIE codes return an average value for each parameter which
can be used as an average model or these values can be used to initialize SIR and LILIA.
         This paper discusses the theoretical calculations that we used and then presents the outcome of our

77 NSO REU/RET Annual Program Report 2004                                      APPENDIX II – REU Student Reports
2. Theoretical Calculations
          To calculate profiles for three Fe I lines at 630.15, 630.25, and 1564.85 nm we solved for the emergent
intensities of the four Stokes parameters I, Q, U, and V in two dimensions. We preformed this calculation twice,
once for the visible region containing the 630.15 and 630.25 lines and again for the infrared line. The central
wavelength, average formation height in a standard one-dimensional model, and effective Landé factor for each
line is shown in Table 1. Formation height is defined as the height above the photosphere at which optical depth
at the center of the line is unity. For more information on the calculations used in this study, see Uitenbroek 2003
and references therein.

                                                     TABLE 1
                                     λ0 [nm]          hFALC [km]        gLandé
                                     630.15             220.2           1.667
                                     630.25             141.7           2.500
                                     1564.85            144.0           3.000

        The model that we calculated is a cross-section through a single snapshot from a magneto-hydrodynamic
simulation of the solar convection. It has 64 depth points and covers 6 Mm with a resolution of 0.04 pixels per km.

                 Fig. 1: Fit of an abnormal Stokes V profile from MELANIE (dotted) and SIR (dashed).

3. Inversion Codes
         In this study we test two types of inversion codes, Milne-Eddington codes that return an atmospheric
model constant with depth and more complex codes that return line of sight gradients. All of the codes that we
tested used minimization of a merit function to find the best solution. Although M-E codes are much faster, the
simple structure that they impose cannot fit abnormal Stokes profiles that are regularly observed and that our
theoretical model produced. For instance, figure 1 shows a plot of an abnormal profile calculated by the radiative
transfer code compared with the profiles returned from the URI code which gives an M-E type model as well as
SIR which provides line of sight gradients. Although M-E codes are not practical for extracting accurate values
they can be used to initialize codes in order to find the best possible solution.

3.1 Milne-Eddington Codes
         A Milne-Eddington code makes the assumption that the source function varies linearly with τ, giving
analytic relations between the shapes of the Stokes profiles and the structure of the atom. When the values for
each parameter, which are constant over height, are substituted into the radiative transfer equations, the equations

78 NSO REU/RET Annual Program Report 2004                                        APPENDIX II – REU Student Reports
 Fig. 2: Scatter plot of RH calculated data (x-axis) and data from URI (y-axis) for magnetic field strength, velocity,
 inclination and azimuth.

can be solved analytically as shown by Jefferies, Lites, and Skumanich 1989. By solving the equations
analytically, computational time is drastically reduced. We test two codes based on these assumptions in this
          The URI (modified Unno-Rachkousky Inversion code) data analysis system returns magnetic field
strength, angle of inclination, azimuth, line strength, damping, velocity, source function and its gradient, and
filling factor. The returned values are very close to the average of each parameter over the height of formation as
calculated in the theoretical model. Figure 2 shows a comparison of theoretical values plotted on the x-axis and
URI values plotted on the y-axis. Perfect correlation is represented when every point falls on the 45o line. The
average inversion time is 1.2 seconds per profile on a Sun Workstation with a Sparc 3 750 MHz processor,
making this the fastest code that we tested.
          The other M-E code that we tested is MELANIE which is part of the Community Inversion Codes from
the High Altitude Observatory. MELANIE also provides reasonable averages; however, the code returned
unacceptable error bars that were plus or minus 50% of the value of each parameter at best. The reason for these
error bars is unknown, perhaps bugs in the code.

3.2 Inversion through minimization
         Unlike codes based on M-E models, SIR (Stokes Inversion based on Response Function) derives line of
sight gradients. SIR does this by starting with an atmosphere that is constant with depth and then perturbing each
parameter in turn to find the best fit to the observed profiles. The validity of each perturbation is evaluated using a
merit function, which is the sum of the squared differences between the observed and synthetic profiles weighted
by uncertainties (Bellot Rubio 1999 and Ruiz Cobo & del Toro Iniesta 1992). This causes the computational time
to increase to an average of 30 seconds per profile on the same sun workstation used in section 3.1.

79 NSO REU/RET Annual Program Report 2004                                        APPENDIX II – REU Student Reports
         LILIA, also part of the Community Inversion Codes, works on the same premise as SIR using slightly
different routines. Even thought SIR and LILIA are very similar, LILIA takes twice as much computational time.
The time for inversion averaged 2.8 minutes per profile. Even though the results are almost identical, this delay
when inverting a large number of profiles makes all the difference.

3.3 Initialization Methods
          We tried several initialization schemes in order to get the best possible results with SIR and LILIA.
Initialization is choosing the model that is used as a guess for the inversion code. If the model is too different
from the true atmosphere there is a risk that the minimization of the merit function will not find the correct
solution but will only find a secondary minimum. Because the merit function has so many variables (each
parameter that is being inverted times the number of nodes for each parameter) it is easy for the code to fall into a
secondary minimum.
          We first inverted our profiles with LILIA and SIR using a standard HSRA (Harvard Smithsonian
Reference Atmosphere) model for the quiet sun as the initializing model. The temperature was accurately
inverted over the height of formation of the line; however, the other parameters were much less accurate. Figure 3
shows a scatter plot of each parameter verses the calculated model. The magnetic field was very poorly
recovered, almost zero for all points.

   Fig. 3: Comparison of SIR and theoretical values with HSRA initialization. Black points represent SIR,
   green points represent LILIA

80 NSO REU/RET Annual Program Report 2004                                       APPENDIX II – REU Student Reports
 Fig. 4: Comparison of SIR and theoretical values with URI values as initialization. Black points represent SIR,
 green points represents LILIA.

          We then initialized using values from the URI inversion code for magnetic field, velocity, angle of
inclination, and azimuth. The resulting values were in much better agreement with the original values than with
the first initialization. Again the temperature was accurately recovered and the other parameters are closer to the
actual values from the simulation as seen in Figure 4. The azimuth is the least accurate of the parameters that we
compared. We believe this is due to an 180o ambiguity that can occur in inversion codes.
          In another attempt to better recover magnetic fields, we tested both inversion codes for the Fe I line in the
infrared at 1564 nm. This line has a Landé factor g=3 so the effects of the magnetic field should be more
pronounced than for the 630.15 and 630.25 nm lines which have smaller of the calculated values across the
formation height for magnetic field, velocity, inclination, and azimuth initialization. We used average values
because the URI code is not yet capable of inverting the IR line. Because the URI values in the visible lines and
the averages were close, we assumed that the values of URI for the IR would be equally close; however, URI
should be tested in the inversion of this line. Figure 5 shows very similar results to the visible case. SIR found
higher magnetic fields than we calculated in our simulation in this initialization, more than with other guess

81 NSO REU/RET Annual Program Report 2004                                        APPENDIX II – REU Student Reports
         Fig. 5: Inversion of IR line using average values for initialization. Black is SIR, green is LILIA.

         We then looked at the correlation coefficients for each height of the parameters. As was expected the
correlation peaks roughly over the height of formation as seen in Figure 6 for most parameters. However, in both
the LILIA and SIR inversions the temperature and electron density have low correlation over the height of
formation. To test if this came from poor initialization we ran SIR with a guess model containing all the actual
values from the calculation. Figure 7 shows that there is still a dip even when the best possible initialization is
         One idea on the cause of this is that SIR and LILIA assume a non-dynamic atmospheric model, guessing
that pressure is the only force supporting the atmosphere. Unfortunately, this is not the case in the real solar
atmosphere and in our model. The theoretical model takes into account dynamic interaction of particles with the
atmosphere. This difference in assumptions could be the cause of the dip in temperature and electron density.

3.4 Center of Gravity Method
         Another type of inversion that can be used to derive line-of-sight magnetic field and velocity values that
are constant with depth is the center-of-gravity method. The center of gravity (cog) wavelength of a line profile is
defined as the centroid of its residual intensity profile. Using the cog wavelength the Doppler shift can be
calculated and the velocity can be inferred. The magnetic field is found by finding the center of gravity for I±V
profiles. This method is studied in depth by Uitenbroek 2003.
         The center of gravity method could be used for initialization in the same way the URI code is used in this

82 NSO REU/RET Annual Program Report 2004                                       APPENDIX II – REU Student Reports
   Fig. 6: Correlation coefficients versus height for physical parameters initializing with URI. Solid lines are SIR,
   dashed lines represent LILIA. Shows lower correlation for temperature and electron density.

   Fig. 7: Correlation coefficients versus height for physical parameters initializing with complete model using
   SIR inversion. Temperature and electron density still show low correlation.

83 NSO REU/RET Annual Program Report 2004                                      APPENDIX II – REU Student Reports
4. Conclusions
         Inversion codes should be used with care when interpreting observed profiles. Even when good fits of the
Stokes profiles are obtained, the derived atmospheric values can still be inaccurate. In order to find more accurate
results good initialization is needed, such as provided by a M-E inversion code to find average values. However,
even with good initialization, values are still inaccurate at certain heights. As was discussed in section 3.3,
temperature and electron density have very poor correlation with the theoretical model even at the height of
         With new technology being developed that will take more and more data, the process of data reduction
using inversion codes will become heavily relied upon. In order for this to be possible faster and more dependable
codes must be developed. At the current rate of about two minutes per profile it would take on the order of
months for large data sets that are now being produced to be completely inverted. One method that is being
explored to quickly invert large data sets is the principle components analysis (PCA) technique (Socas-Navarro et
al. 2001). This technique uses a database of pre-computed profiles and their models to match the observed
profiles. Using FATIMA (fast analysis technique for the inversion of magnetic atmospheres) to implement the
PCA technique the computational time was cut down to minutes compared to the hours it would take using a M-E
         The uncertainty in the results obtained with codes like SIR and LILIA continues to be a problem. There
is always a danger when using minimization routines for the code to find a secondary minimum instead of the best
possible solution. Even when initializing with values from a M-E code the retrieved values can be unsatisfactory
depending on the level of accuracy needed.

This work was supported by the National Science Foundation (NSF) and the Association of Universities for
Research in Astronomy. This work is carried out through the National Solar Observatory Research Experiences
for Undergraduate (REU) site program, which is co-funded by the Department of Defense in partnership with the
National Science Foundation REU Program.


Bellot Rubio, L. R. (1999): A user guide to SIR, Instituto de Astrofísica de Canarias

Community Inversion Codes. (2000): <>.

Jefferies, J., Lites, B. W., & Skumanich, A. (1989): Astrophysical Journal, 343, p.920-935

Ruiz Cobo, B., & del Roro Iniesta, J. C. (1992): Astrophysical Journal, 398, p.375-385

Socas-Navarro, H., Lόpez Ariste, A., & Lites, B. W. (2001): Astrophysical Journal, 553, p.949-954

Uitenbroek, H. (2003): Astrophysical Journal, 592, p.1225-1233

84 NSO REU/RET Annual Program Report 2004                                      APPENDIX II – REU Student Reports
                                   To be submitted to Sol. Phys.– basic format is ApJ, but ref-
                                   erences are in Sol.

                   Solar Wind Forecasting with Coronal Holes
                                          Stuart Robbins1
      Department of Astronomy, Case Western Reserve University, Cleveland, OH 44106


                                   C. Henney2 and J. W. Harvey2
                          National Solar Observatory, Tucson, AZ 85719



            An empirical model for forecasting solar wind speed related geomagnetic
        events is presented here. The model is based on the location and size of so-
        lar coronal holes determined with Kitt Peak Vacuum Telescope He I 1083.0 nm
        spectroheliograms and photospheric magnetograms. This method differs from
        the Wang-Sheeley model that is based on photospheric magnetograms to esti-
        mate the open field line configuration. Solar wind and coronal hole data for the
        period between May 1992 and September 2003 is investigated. The new model is
        found to be accurate to within 4.5−5.7% (the range depends upon the number of
        days ahead forecast) of observed solar wind measurements for the best one-month
        periods within the time frame studied. Using coronal hole maps, the model can
        predict the solar wind velocity up to 8.5 days in advance with an average de-
        viation as low as 9.4 − 10.0% for a given one-month period. This is further in
        advance forecasting and up to a factor of 2 improvement over the Wang-Sheeley
        model. Its main features are a strong southern hemisphere bias, sunspot cycle
        dependence, and that a complete forecast can be made from a single solar image,
        as opposed to a full synoptic map required by the Wang-Sheeley method.

        Subject headings: Sun: coronal holes, Sun: solar wind

    REU 2004 summer intern with the National Solar Observatory. The REU program is funded by the
National Science Foundation.
      National Solar Observatory

                                         1.   Introduction

     The solar wind has an important impact upon technology today, for it influences the
longevity and the survival of all satellites outside and within Earth’s magnetosphere as well
as many ground-based electronics. In today’s world, communications, power grids, and even
trains are all at the mercy of space weather, governed by the ever-changing solar wind (Jansen
& Pirjola 2004).
     Because of the multi-billion dollar industries that rely upon these electronic devices and
the quasi-static nature of the atmosphere and magnetosphere, various methods have been
created to shield them from dramatic increases in solar wind output and resulting effects;
however, many of these rely upon a fore-knowledge of increased output, so an accurate
forecast of the solar wind velocity at Earth is essential.
     Presented here is an expansion of the 1970s work of Sheely and Harvey (e.g. Sheeley
& Harvey (1981)) to (1) show that coronal holes can be used to as accurately forecast
solar wind velocities as synoptic maps and magnetic field propagation, (2) demonstrate the
latitude and longitude weighting that leads to the best forecast for the entire 11-year period
studied, (3) illustrate the dependence upon the phase of the sunspot cycle, (4) show coronal
holes can also be used to forecast the interplanetary magnetic field polarity at Earth, and
(5) provide example 1-month period forecasts using this model.

                                       1.1.   Coronal Holes

     Krieger, Timothy, & Roelof (1972) used data from 1970 and traced solar wind streams
observed at Earth back to the Sun; they found a good agreement between high-speed solar
wind velocities and the positions of coronal holes (see fig. 7, Krieger, Timothy, & Roelof
(1972)). Extensive work by Sheeley and Harvey (e.g. Sheeley & Harvey (1981)) throughout
the 1970s further solidified the link between coronal holes, solar wind speeds, and geomag-
netic disturbances.
     Coronal holes are characterized by the following: Compared with the surrounding
corona, the holes are lower in density and temperature1 , and their magnetic fields are open,
weak, and divergent (Altschuler, Trotter, & Orrall 1972; Harvey & Recely 2002). They
are the main sources of intense solar wind that most seriously impact geomagnetic events
(Harvey & Recely 2002).

      Munro & Withbroe (1972) found the coronal temperature is lowered by ∼600,000 K.

     Coronal holes are mainly observed in x-ray, extreme ultra-violet, and radio wavelengths,
but they can also be seen in the two chromospheric spectral lines He I D3 and 1083 nm, and
also on the limb using Fe XIV 530.3 nm and K-coronagraphs (Harvey & Recely 2002). Work
here was done with images based upon 1083 nm images.
     The holes form when the remnants of active magnetic field regions, emerging in both
north and south hemispheres over several months, combine and form a relatively large region
of a unipolar field. The remnants of the opposite polarity fields surround the holes, forming
their divergent magnetic nature (Timothy, Krieger, & Vaiana 1975).
     There are three major classifications for coronal holes: Polar, non-polar, and transient
(Harvey & Recely 2002). In general, polar coronal holes exist for many years, forming from
non-polar holes that drift pole-ward and merge, generally soon after a solar maximum. They
reach their largest area around the time of solar minimum, where they can occupy as much
as 15% of the disk, and they wane and disappear around the time of the next maximum
(Harvey & Recely 2002). Polar coronal hole boundaries distort over their lifetimes, and this
is due to sunspot activity (Maravilla et al. 2001).
     Non-polar holes are more isolated, and they are limited to a latitude range of ±60◦ .
They are associated with an active or formerly active magnetic region and generally persist
for a few months - over one to several rotations (Harvey & Recely 2002). Transient holes
have briefer lifetimes, generally lasting for only one to two days; they are associated with
eruptive events such as flares and coronal mass ejections (Harvey & Recely 2002).

                               1.2.   Wang-Sheeley Model

     Wang & Sheeley (1990) discovered an empirical relation between solar wind velocity
and coronal flux-tube expansion. They found that both slow and fast solar winds have their
origins in coronal holes and that their velocities are related to the open field line divergence
rate close to the Sun. Their results suggested that synoptic maps can be run forward in time
under the influence of differential rotation, supergranular diffusion, and meridional flow; the
resulting simulated photospheric field can then be used to forecast the wind’s velocity at
     The Wang-Sheeley model was revised in 2000 by Arge & Pizzo (2000) with five main
modifications including the establishment of a continuous empirical function that relates
magnetic expansion factor to solar wind velocity at the source surface, propagation of the
solar wind from the source surface to Earth using assumptions about radial streams and
their interactions, and optimization of the technique for the construction of daily updated

synoptic maps from magnetogram data.
     Arge & Pizzo (2000) studied a three-year period centered about the May 1996 solar
minimum. Their three-year sample period had an overall correlation of ∼ 0.4 with observed
solar wind velocities and an average fractional deviation2 of 0.15. When excluding a 6-month
period with large data gaps, they correctly forecast the solar wind to within 10 − 15%.
Interplanetary magnetic field (IMF) polarity was correctly forecast ∼ 75% of the time.

                                        2.   Model Input Data

     Solar wind velocity data was obtained from the National Space Science Data Center’s
OMNIWeb website. Daily averages were used because that was the approximate sampling
rate of the coronal hole maps.
     Coronal hole data from May 28, 1992 through September 25, 2003 (JD 2448771 −
2452908) was used; this corresponds to the last half of Cycle 22 and the first half of Cycle
23. Hand-drawn maps based upon the KPVT 10 830 ˚ He I images and photosphere
magnetograms were used (Harvey & Recely 2002). The images were mapped into sine-
latitude to create square images of equal area.
     The images were keyed in such a way that the background solar disk was set to a value
of 0, and each coronal hole was labeled with an integer. Each hole was signed to correspond
to the polarity of the hole. In this manner, the third hole drawn on a map would be given a
value of 3. If it had a negative polarity it would be changed to −3. The total magnetic flux
of each hole was written into the header of the file.

                              3.    Solar Wind Correlation Analysis

     The solar wind velocity data was taken as-is from the OMNIWeb website. For the time
period studied, it had 91.8% coverage. Missing data was interpolated using a cubic spline
    Each coronal hole image was read into a data array. 14◦ -wide latitude regions from
90◦ N to 90◦ S were taken, masking out all the surrounding pixels; this corresponded to
approximately a 24-hour rotation on the Sun as seen from Earth. A 1-day-wide window
was selected because anything much wider resulted in mixed temporal signals. The absolute

  2                                               prediction−observed
      Average fractional deviation is defined as         observed        .

value of the signed mean of each swath was taken to yield a percent coverage of the area by
coronal holes. This gave a single number tied to the date of the coronal hole.
     The holes were not scaled by their magnetic flux. This was initially investigated, but it
was found that the eventual correlation was substantially increased when the flux was not
included and area was the only factor. Wang & Sheeley (1990) discovered the same effect.
     This data set was then interpolated into the time frame of the solar wind velocity data
in order to eliminate uncertainties generated by the temporal non-uniformity in which the
coronal hole data were sampled. Because there was 68.6% coverage by the coronal holes,
dates in the final time series that did not have a coronal hole map within 0.5 days were
masked for all of the analysis except when cross-correlations were calculated
    Cross-correlation analyses (see appendix A) between the two data sets were then per-
formed to determine if and how they were related.

                                  3.1.   Data Set Scaling

     Because the coronal hole data set was the percent of the solar disk covered by coronal
holes, it needed to be scaled in order to give a physical value for solar wind velocity. Scaling
does not affect the cross-correlation, so this allowed for the development of an appropriate
weighting scheme and then of a scaling mechanism.
     A scatter plot of the coronal hole versus the velocity data was constructed, and a linear
regression was then calculated. The coronal hole data set was multiplied by the slope of the
regression, and the vertical offset was then added into the set. The resulting time series was
mathematically the best fit, but the mathematical best scaling was non-physical in that it
did not result in the range of velocities observed.
     Therefore, scaling was chosen by a process of minimizing the average fractional deviation
and standard deviation while maintaining an accurate range for forecast values, assuming a
linear relationship between the two. The scaling that yielded the best agreement between
the two sets was nearly identical between the three forecasting windows, so the average was
taken. This scaling was then used in the sections that follow. With this scaling and the best
weighting (found in the following sections), the resulting average fractional deviation for the
entire 11-year data set is 1.2% with a 1-σ uncertainty of ±29.5%.

                 3.2.   Longitude Analysis and Forecasting Windows

    23 longitude swaths were examined, starting with the first centered at 77◦ E and going
through to the last centered on 77◦ W. Each of these were processed in the manner described
    One result is the quantification of how long solar features on the disk take to impact
Earth, based upon longitude:

                          d = (3.737 ± 0.022) − (0.07390 ± 0.0005) θ                         (1)

where d is the delay in days and θ is the center of the longitude swath in degrees as measured
from the central meridian (East is negative, West is positive); uncertainties are 1-σ. The
data are fit to this line in fig. 1.
     Figure 2 shows the maximum cross-correlation coefficient that corresponds to the best
time lag each longitude swath. Statistical significance is shown in table 2. This shows three
relatively good prediction windows that were later used to forecast. They are centered on
63◦ E, 14◦ E, and 7◦ W. These yield forecasts of 8.5 days, 5 days, and 3.5 days, respectively.

                           3.3.   Latitude Weighting Analysis

     To increase the cross-correlation, the Sun was further divided into latitude bins; 60◦
increments were used, generating Northern (90◦ N - 30◦ N), Equatorial (30◦ N - 30◦ S), and
Southern (30◦ N - 90◦ S) regions. A combination of Northern with Equatorial and Southern
with Equatorial were used in order to have a total of five bins. The cross-correlation of each
latitude bin were calculated for each of the three windows found in §3.2.
     Figure 3 shows there is a very strong bias towards the southern hemisphere. For the
entire 11-year time series, the correlation for the combination of the Southern with Equatorial
region is almost a factor of 2 greater than that for the Northern with Equatorial region.
     This southern bias was explored further by tiling 40◦ swaths with an overlap of 20◦ across
the disk for a total of 8 regions; three extra were calculated in order to verify the peak’s
position, and the results are illustrated in figure 4. It shows that the southern bias is centered
between 10◦ S to 50◦ S, and that the benefit of going farther South decays rapidly. Going
North, the fall-off in correlation is much more gradual, and it crosses into an anti-correlation
around 40◦ N.

     Based upon this, the final latitude window used spans 30◦ N to 70◦ S. This gave a
maximum cross-correlation coefficient of 0.382 for the 8.5-day forecast, 0.347 for the 5-day
forecast, and 0.377 for the 3.5-day forecast.

                     3.4.   Magnetic Activity Cycle Dependence

     Sheeley and Harvey’s work from the 1970s (e.g. Sheeley & Harvey (1981)) hinted at a
dependence upon the sunspot cycle for the correlation between the coronal holes and solar
wind. This was quantitatively explored for the time period the coronal hole data set spans -
the last half of Cycle 22 and the first half of Cycle 23. The full time series was divided into
six time series of approximately 690 days, shown in table 1 and illustrated in figure 5 with
sunspot counts overlaid.
     Figure 6 shows the results of this analysis. There is a strong dependence upon the phase
of the sunspot cycle, similar to the qualitative results of Sheeley and Harvey. The correlation
is best during and just after solar maximum, and it is worst during solar minimum and the
beginning ascending phase of the Cycle.

                               4.   IMF Polarity Forecast

    In addition to forecasting the solar wind velocity, coronal hole maps can be used to
predict the interplanetary magnetic field (IMF) polarity at Earth. The 3.5-day forecast
window at 7◦ W was used for this.
     The coronal hole maps were scaled so that each pixel with a positive polarity hole had a
value of +1 and each pixel with a negative polarity hole had a value of −1; the background
solar disk was set to 0. The mean of all the pixels in the 14◦ -wide region was taken, with a
latitude range of ±60◦ .
    Excluding days when the average was very close to 0 and the two data sets both did
not have actual measurements, 47.3% of the days in the 11-year time series were forecast; of
those, 73.6% of them had the IMF polarity correctly forecast. When only excluding days that
did not have both velocity and coronal hole data, the IMF polarity was correctly forecast
68.7% of the time, corresponding to 62.9% of the days having forecasts. This is a drop of
5% in accuracy, but an increase by 15% of the days forecast.
    This is slightly lower than the results Arge & Pizzo (2000) – they correctly forecasting
the polarity 75% of all the days in their time period.

                                 5.   Example Forecasts

    There are several ways to compare this empirical model with the modified Wang-Sheeley
(Arge & Pizzo 2000). The first is by looking at cross-correlations.
     Arge & Pizzo (2000) studied a three-year period. For this period, they found - using
their best forecast method - a correlation of 0.389. For the same time length of three years,
the model presented here has a maximum cross-correlation of 0.439 − 0.473 with the range
depending upon the length of the forecast made.
    The best one-month period presented by Arge & Pizzo (2000) had a cross-correlation
coefficient maximum of 0.813. The best one-month period within the 11-year series that was
studied for this paper had a correlation of 0.864 − 0.894 with the range again dependent
upon which of the three forecasts was used.
     A second way to compare this empirical model as an effective forecast tool is through
comparing the average fractional deviation (AFD), defined as prediction−observed , where a
value of 0 indicates a perfect match. Arge & Pizzo (2000) had an AFD that averaged
at 0.159 with a range of about 0.15 for the entire three-year period. The best AFD for
the best one-month period presented in their paper is 0.096. For the entire 11-year period
studied in this paper, the AFD had a mean of 0.00059 − 0.00740 with a standard deviation
of 0.293 − 0.306. The best AFD of a one-month period studied here is 0.00055 − 0.00126.
     Finally comes the practical test of forecasting the solar wind. Figures 7 and 8 present
two sample forecast periods spanning approximate one-month time frames. Fig. 7 shows
one of the best periods forecast while fig. 8 represents one of the worst periods for this
forecasting method.
     Besides the purpose of demonstrating the practicality of this method through sample
forecasts, an important statistical point was discovered. Arge & Pizzo (2000) argue that
the correlation coefficients are sometimes poor indicators of how well the model predicts the
solar wind velocity. However, when searching for sample periods to include in this by looking
at the lowest AFD, as they argue is the better indicator, it was found that the AFDs closest
to 0 were for time periods where the solar wind velocity was near a minimum, but still varied,
and the forecast was a flat line (such as illustrated around JD 2450350 in fig. 8). In order
to find periods where the forecast mirrored the observed values, the cross-correlations were
     This is illustrated within the two forecasts themselves. The correlations for the good
period are between 0.748−0.807 while those for the poor period are 0.069−0.207. Significance
levels are given in table 2. The AFD for the good period is in the range −0.167−−0.124 while

that for the poor period lies between −0.116 − −0.0935. When ±σ uncertainties are applied
(see §3.1; figures not shown due to size and readability constraints), the good period is
correctly predicted for every day forecast, while the poor period has over one dozen forecasts
that do not cover the observed values. It is thus concluded that, in general, correlation
coefficients are more telling of the goodness of forecast than AFD.

                                      6.   Discussion

     For the 11 years studied, this empirical model is a quantitative improvement over one of
the most widely used forecasting methods of solar wind velocity - the Wang-Sheeley. Both
are based in the empirical relation between coronal holes and the solar wind, but the Wang-
Sheeley treatment is to propagate the magnetic fields to create the forecast. The method
presented here is with simple weighting of the solar disk in order to create a forecast for the
solar wind velocity at Earth.
      There are several important and unique characteristics of this empirical model. First,
it relies upon a strong southern weighting that has yet to be explained. Weighting the south
while removing the northern regions of the Sun results in correlations over a factor of 2
higher. The final weighting shows higher cross-correlations than those found in the leading
modified version of the Wang-Sheeley model, presented in Arge & Pizzo (2000). Second,
there is a strong dependence upon the phase of the sunspot cycle. This has been found in
the past, and we hope to explore it further when more data becomes available in the future.
It can also be used to predict the magnetic field polarity at Earth from the sun.
     Another important feature of this model is its lack of dependence upon synoptic maps,
and that a forecast can be created from a single image of the solar disk. There are three
prime forecast windows centered at 63◦ E, 14◦ E, and 7◦ W that offer forecasts of 8.5, 5, and
3.5 days respectively, but any 14◦ -wide window can be used and the forecast in days can
be obtained for that via eq. 1. With multiple prime forecast windows, this method is less
sensitive to gaps in observations because of the overlap that is offered of forecast times with
gaps in coronal hole sampling.

     Solar wind data is available online at This work
is carried out through the National Solar Observatory Research Experiences for Undergrad-
uate (REU) site program, which is co-funded by the Department of Defense in partnership
with the National Science Foundation REU Program. This research was supported in part
by the Office of Naval Research Grant N00014-91-J-1040. The National Solar Observatory
is operated by AURA, Inc. under a cooperative agreement with the National Science Foun-
                                                  – 10 –

    Facilities: KPNO.

                                  A.        Cross-Correlation

    Equation A1 is the basic form of the cross-correlation used in this study.

                                     N−|L|−1
                                               (xk+|L| −¯)(yk −¯)
                                                         x      y
                                      k=0
                                  u                    !                    !    forL < 0
                                 u    P
                                      N−1                   P
                                 t              x
                                            (xk −¯)2
                                                                 (yk −¯)
                                                                      y 2
                                      k=0                  k=0
                      ρxy =                                                                              (A1)
                                                     x        y
                                                (xk −¯)(yk+L −¯)
                                 v
                                                       !                    !    forL ≥ 0
                                 u    P
                                      N−1                   P
                                 t              x
                                            (xk −¯)2
                                                                 (yk −¯)
                                                                      y 2
                                      k=0                  k=0

    In equation A2, PN (|r| ≥ |r0 |) gives the probability that with N observations, two
uncorrelated variables would give a coefficient r as large as r0 . Therefore, a correlation was
considered to be significant if the coefficient obtained fell in PN ≤ 5%. If it was less than
1%, then the correlation was highly significant.

                                     2Γ ((N − 1) /2)                    1
                                                                                         (N −4)/2
                 PN (|r| ≥ |r0 |) = √                                           1 − r2              dr   (A2)
                                      πΓ ((N − 2) /2)                 |r0 |

     For comparisons with other studies and to demonstrate features of this model, the
significance of several values for N were calculated, and they are listed in table 2.


Altschuler, M. D., Trotter, D. E., and F. Q. Orrall: 1972, Sol. Phys., 26, 354.

Arge, C. N. and V. J. Pizzo: 2000, J. Geophys. Res., 105:A5, 10,465.

Harvey, K. L. and F. Recely: 2002, Sol. Phys., 211, 31.

Jansen, F. and R. Pirjola: 2004, Eos, 85:25, 241.

Krieger, A. S., Timothy, A. F., and E. C. Roelof: 1972 Sol. Phys., 29, 505.
                                              – 11 –

Maravilla, D. et al.: 2001, Sol. Phys., 203, 27.

Munro, R. H. and G. L. Withbroe: 1973, ApJ, 176, 511.

Sheeley, N. R., Jr. and J. W. Harvey: 1980, Sol. Phys.70, 237.

Timothy, A. F., Krieger, A. S., and G. S. Vaiana: 1975, Sol. Phys., 42, 135.

Yang, Y.-M. and N. R. Sheeley, Jr.: 1990, ApJ, 355, 726.

   This preprint was prepared with the AAS L TEX macros v5.2.
                                                   – 12 –

                                                        Data Points
                     -8                                 Linear Regression
                                                          [Slope = +0.0739 Days/Degree]
       Day Offset




                          -70   -56   -42   -28   -14   0    14    28   42    56    70


Fig. 1.— Day Offset with Solar Disk Latitude: A linear relation was found in the affects
of the coronal holes within 14◦ -wide swaths of the solar disk upon the solar wind velocity.
Shown here are the data for this plotted along with the best linear fit. Negative longitude
corresponds to East and positive to West.
                                                                              – 13 –

       Maximum Cross-Correlation Coefficient






                                                     -70   -56   -42   -28   -14   0   14   28   42   56   70


Fig. 2.— Longitude Cross-Correlations: 14◦ -wide swaths of the solar disk were tiled
across the solar surface with 7◦ overlap between successive swaths. The maximum cross-
correlation (corresponding to a day offset described by eq. 1 and shown in fig. 1) of the
resulting data set with the solar wind velocity data are shown here. Negative longitude
corresponds to East and positive to West. This shows three relatively good forecast windows
centered on 63◦ E, 14◦ E, and 7◦ W.
                                                                     – 14 –

        Amplitude of Correlation Maxima

                                                     Centered at 63° E
                                          0.4        Centered at 14° E
                                                     Centered at 7° W


                                          0.1                                           Full Latitude
                                                                                         Centered at 63° E
                                          0.0                                            Centered at 14° E
                                                                                         Centered at 7° W
                                                 N         N+E                E           S+E                S

                                                                         Latitude Bin

Fig. 3.— Five Latitude Regions: The three longitude forecasts’ cross-correlations with
the solar disk divided into five latitude regions described in §3.3. This clearly indicates there
is a strong preference for the southern region of the Sun for an accurate forecast of the solar
wind. This is refined in fig. 4.
                                                                     – 15 –

       Amplitude of Correlation Maxima

                                                     Centered at 7° W - 40° Latitude Correlation
                                         0.4         Centered at 7° W - Full Latitude Correlation





                                         90°N-50°N       50°N-10°N                   10°S-50°S      50°S-90°S

                                                                      Latitude Bin

Fig. 4.— 40◦ Latitude Swath Correlations: Refined is the 7◦ W solar region for 11
different latitude regions, each 40◦ wide (described in §3.3). From this, the southern bias
that was discovered with fig. 3 is refined to show that it is centered around 30◦ S. The rapid
fall-off south of this led to the refinement of the final latitude window to 30◦ N through 70◦
S. This region was chosen in order to avoid too much tailoring to a specific time period.
                                            – 16 –

                                                Sunspot Number - Monthly Avg.
                         200                    Sunspot Number - Yearly Avg.
        Sunspot Number




                           1988   1992         1996            2000           2004


Fig. 5.— Sunspot Cycle: Both monthly and yearly sunspot number averages are shown
here along with six time bins into which the 11-year data set was divided; the six bins are in-
dicated by the rounded rectangles, are numbered 1−6 left to right, and lasted approximately
690 days each (see table 1).
                                                                   – 17 –

       Amplitude of Correlation Maxima

                                                                  Centered at 63° E - N+E
                                                                  Centered at 14° E - N+E
                                         0.4                      Centered at 7° W - N+E



                                               1       2          3            4            5   6

                                                                  Time Series Bin

Fig. 6.— Cycle Correlations: Divided into the six time bins (see fig. 5 and table 1),
the cross-correlations of the coronal hole data with the solar wind follows the sunspot cycle
well, with the correlation being very high during solar maximum and the beginning decline
from maximum, and the correlation reaching a minimum during solar minimum. Each of
the three forecast windows are shown here; the window centered on 7◦ W shows the widest
variation in correlation, going from nearly 0.5 to statistical insignificance (see Appendix A
and eq. A2), while the window centered at 14◦ E shows the least variation with the Cycle.
                                                             – 18 –

                                                               Solar Wind Velocity
        Solar Wind Velocity (km/s)

                                     800                       8.5-Day Forecast
                                     700                       5.0-Day Forecast
                                                               3.5-Day Forecast




                                           2451600       2451610         2451620     2451630

                                                     JD (39 Days) Around CR 1960

Fig. 7.— Good Forecast Period: JD 2451594 − 2451632 has a cross-correlation coefficient
of 0.748, 0.807, and 0.793 for the 8.5-, 5-, and 3.5-day forecasts respectively. For this number
of days, the correlation is considered significant if above 0.32 and highly significant if above
0.41. The AFD during this period is −0.167 ± 0.158, −0.124 ± 0.186, and −0.124 ± 0.161.
This period is from time bin 5 (see fig. 5).
                                                           – 19 –

                                                                     Solar Wind Velocity
        Solar Wind Velocity (km/s)

                                     800                             Interpolated SW Velocity
                                                                     8.5-Day Forecast
                                                                     5.0-Day Forecast
                                     600                             3.5-Day Forecast




                                       2450340   2450350       2450360        2450370

                                                   JD (39 Days) Around CR 1914

Fig. 8.— Poor Forecast Period: JD 2450339 − 3450377 has a cross-correlation coefficient
of 0.207, 0.121, and 0.0685 for the 8.5-, 5-, and 3.5-day forecasts respectively. For this number
of days, the correlation is considered significant if above 0.32 and highly significant if above
0.41. The AFD during this period is −0.0935±0.1771, −0.1161±1619, and −0.0976±0.2055.
This period is from time bin 3 (see fig. 5).
                                            – 20 –

Table 1: Time series   subdivisions of full JD 2448771 − 2452908 used in cross-correlations.
 Julian Date Range      Date                           Days Days with Data Completeness
 2448771 − 2449460      May 28, 1992 - Apr 17, 1994 690             456              66.1%
 2449461 − 2450150      Apr 18, 1994 - Mar 7, 1996      690         472              68.4%
 2450151 − 2450840      Mar 8, 1996 - Jan 26, 1998      690         498              72.2%
 2450841 − 2451530      Jan 27, 1998 - Dec 17, 1999     690         509              73.8%
 2451531 − 2452220      Dec 18, 1999 - Nov 6, 2001      690         446              64.6%
 2452221 − 2452908      Nov 7, 2001 - Sep 25, 2003      688         456              66.3%
 2448771 − 2452908      May 27, 1992 - Sep 25, 2003 4138           2838              68.6%
                                          – 21 –

Table 2: Significance for various time series lengths, as determined through eq. A2.
 Number of Days Significant if Above Highly Significant if Above
        30                  0.36                       0.46
        39                  0.32                       0.41
        365                 0.10                       0.14
        688                0.074                      0.098
        690                0.074                      0.098
       4138                 0.03                       0.04

To top