AstroGeoVis v1.0: Astronomical Visualizations and Scientific Computing for Earth Science Education

Modern climate science, Earth system science, physical geography, oceanography, meteorology, and related disciplines have increasingly turned into highly quantitative, computational fields, dealing with processing, analysis and visualization of large numerical data sets. Students of these and many other disciplines thus need to acquire robust scientific computing and data analysis skills, which have universal applicability. In addition, the increasing economic importance and environmental significance of solar power and sustainable practices such as passive building design have recently increased the importance 5 of understanding of the apparent motions of the Sun on the celestial sphere, for a wider array of students and professionals. In this paper, I introduce and describe AstroGeoVis v1.0: open-source software that calculates solar coordinates and related parameters and produces astronomical visualizations relevant to the Earth and climate sciences. The software is written in MATLAB©; while its primary intended purpose is pedagogical, research use is envisioned as well. Both the visualizations and the code are intended to be used in the classroom in a variety of courses, at a variety of levels (targeting high school students 10 to undergraduates), including Earth and climate sciences, geography, physics, astronomy, mathematics, statistics and computer science. I provide examples of classroom use and assignment ideas, as well as examples of ways I have used these resources in my college-level teaching.

are listed in Sect. B. They can be used to study alternative approaches and further verify AstroGeoVis v1.0. The code has been extensively tested. The author welcomes suggestions for additions and improvements, and notifications of any possible 60 remaining errors, at the contact information provided.

Solar Declination Student Project
In this section, I present a semester-long student project in which students are asked to measure the altitude and azimuth of the Sun on a weekly basis. The measurements are performed using widely available materials, namely any straight stick-like 65 object (I use a 60 cm PVC pipe) that can serve as a gnomon to cast a shadow, a measurement tape, and a protractor. Students measure the shadow length and the azimuth of the shadow direction. At the end of the semester, students analyze their data and use their measurements to determine the declination of the Sun and track and discuss its progress throughout the semester.
Two variants of the project are presented. In the first variant, students make the measurements at any time of day and record both the altitude and the azimuth of the Sun. An additional requirement for this variant is that the direction of true North needs 70 to be predetermined and known to the students before the start of the measurements. In the second variant, students make the measurements at local solar noon and record only the shadow length. They are tasked to determine local solar noon for the dates and location(s) of their measurements. This variant greatly simplifies the relationship between declination and solar altitude and the interpretation of students' results.
The end of semester analysis for the first variant is accomplished in the MATLAB© script template 75 solar_declination_Exercise.m. Students import their data from a text file and either use the provided formulae to compute solar altitude and declination from their measurements, or they can be tasked with deriving or looking up the formulae.
The altitude of the Sun, h, is determined as: where g is the gnomon length , and s is the shadow length from the base of the gnomon to the tip of the shadow, expressed in 80 the same units as g. The gnomon has to be placed vertically on a horizontal level surface. Subsequently solar declination, δ, is determined using conversion from local horizontal to equatorial coordinates (Vincent (2003)): where h is solar altitude, φ is the latitude of the observer, and A is the measured azimuth of the Sun. Note that here and throughout AstroGeoVis v1.0 I define azimuth as measured clockwise from true North, unlike Meeus (1998) (see his Ch. 13), 85 and as in Vincent (2003).
The analysis for the second variant is achieved using an Excel© spreadsheet, a template for which is provided with As-troGeoVis v1.0. The second variant is more suitable for non-STEM students who may not have the required aptitude and background to deal with scientific code or more advanced trigonometry/geometry. Both variants are provided with comments   The daily path of the Sun on the celestial sphere at a given location and its annual variation have important implications for local insolation and climate, duration of daylight, and thus also for agriculture, solar renewable energy, passive building design and related applications. It is a function of the latitude of the observer and the declination of the Sun, and figures depicting its seasonal variation are common in introductory physical geography textbooks, e.g. Geosystems (Christopherson (2011)).
Here, I present the sun_path.m function, which plots one or more apparent daily Sun paths by first generating a time series of solar horizontal coordinates for a 24-hr period and then converting them to Cartesian coordinates for plotting. One way to call the function is to input a single latitude and one or more declination values (an example is given in Fig. 2A).
The observer is located at the center of the grey semi-transparent mesh sphere, which represents the celestial sphere with the local horizontal coordinate system grid. The local horizon and meridian are drawn, the zenith and nadir points are shown, as 110 is Earth's axis of rotation, indicating the directions of the North and the South Pole. Users can also choose to specify a given scalar solar declination and plot the path of the Sun for several latitudes instead -an example is shown in Fig. 2B. These 3D figures with PTZ capability are valuable to improve students' spatial understanding of the phenomena; inputs can be varied in real time.

115
The terminator is the curve on Earth's along which sunrise or sunset is occurring at a given moment. Here I present a collection of functions visualizing the terminator and creating animations of its progression on Earth's surface over a single day, or over an entire year. An example of a static plot of the terminator is shown in Fig.3. In addition to the line of the terminator, the locations of the sub-solar and anti-solar points are shown, and the meridian at which the Sun is transiting is drawn, as is the opposite meridian 180 • away. In addition, the altitude of the Sun is also mapped.
where α is the right ascension of the Sun. The variables in Eq. 3 must all be consistently expressed in either hours or degrees.
The altitude of the Sun, h, is then determined at every grid point using the conversion from equatorial to local horizontal 140 coordinates Meeus (1998): where φ is the latitude, δ is the declination of the Sun, and HA is the hour angle of the Sun.
Note that the location of the terminator determined by the functions presented here do not take into account atmospheric refraction and consider sunrise and sunset to be the moments when the center of the solar disk is at the mathematical horizon.

145
The terminator plot and some other plots presented subsequently purposefully avoid the use of the term declination. In my experience, this makes them more suitable for use in classes such as Introduction to Physical Geography for non-STEM majors, where the use of technical astronomy terms may be undesirable. Conversely, declination and the fact that it is a governing variable for many of the visualizations and phenomena discussed here should be discussed in astronomy and physics classes, and in more advanced Earth science classes and those focusing on the computational methods.

The Equation of Time and the Analemma
Humans are influenced by a strong circadian rhythm and have used the motions of the Sun on the celestial sphere for timekeeping since ancient times. However, the Sun exhibits irregularities that cause, for example, time-varying corrections to be needed for sundials. These irregularities are captured by the concept of the Equation of Time, E, which is defined as the difference between the hour angle of the true Sun, HA true and the hour angle of the mean Sun, HA mean , (a fictitious Sun that moves at 155 a uniform rate along the celestial equator) (Meeus (1998)): The temporal variability in E is caused by 1) the eccentricity of Earth's orbit, and 2) the obliquity of the ecliptic, resulting in the true Sun not moving uniformly along the celestial Equator, i.e. the right ascension of the true Sun doesn't change at a uniform rate, and there are times of the year when the true Sun is ahead of the mean Sun, and vice-versa.

160
The script analemma.m computes and plots E and some related quantities using the right ascension and declination of the Sun as well as local sidereal time as given by the Meeus (1998)  Hemisphere. When E is plotted against solar declination for an entire year, the analemma curve results (Fig. 5). It has the shape of a figure '8' and is often plotted on globes and in textbooks; it is also often demonstrated as a series of photographs of the Sun taken at the same civilian time at a given location throughout the year.
The time derivatives of the right ascension of the true Sun, α, and the Equation of Time, E, have instructional value. These 170 two quantities are related as follows: where k = ds dt = 1.00273790935 s/s is the time rate of change of sidereal time, i.e. it's the ratio of the duration of the mean solar day to the duration of the sidereal day. Eq. 6 is derived by differentiating Eqs. 3 and 5 with respect to time and substituting one in the other, and noting that the time derivative of HA mean is unity as long as the variables of Eq. 6 are all expressed in 175 the same time unit.
The time series of the two derivatives in Eq. 6 ( Fig. 4B, in minutes per day), illustrate that they are inversely proportional and offset by the values k − 1. The interpretation of these curves is as follows: the time derivative of E is zero when the solar right ascension is increasing at about 3 min 56 s per day (the amount of time by which the mean solar day is longer than the sidereal day), i.e. the true Sun is acting like the mean Sun. If the true Sun's dα dt were always equal to k − 1, i.e. if the true Sun moved 180 uniformly along the celestial Equator, then there would be no temporal variability in E. At times when dα dt > k − 1, the solar day is getting longer and the true Sun is lagging more and more behind the mean Sun, decreasing E (Eq. 5), and vice-versa.
The pedagogical value of the material and visualizations presented in this section is multi-faceted and is perhaps centered on illustrating the relationship between the orbital elements of Earth and Kepler's laws and their manifestation in the apparent motion of the Sun on the celestial sphere, resulting in practical applications, such as timekeeping. Details presented in this 185 section generally go beyond what is suitable in an introductory undergraduate course in Earth Sciences and are more suitable for more advanced classes and students of astronomy, physics, mathematics, and computer science.  As with all AstroGeoVis v1.0 code and visualizations, the gnomon plots and code are dynamic and adaptable and can be used in passive or active learning scenarios, e.g. instructors can just lecture using the figures, or they can ask students to come 200 up with the algorithm, or to investigate the effects of latitude on the traces. Additionally, students can be tasked with designing a sundial using the principles in gnomon_trace.m as a starting point and taking into account the need for correction given by the analemma (Fig. 5 and Fig. 6B). Further details on sundial construction mathematics are given in Meeus (1998), Chapter 58.  solar panel with any orientation, given its latitude and solar declination. Further, the functions compute total energy produced over the course of a year, and solve the optimization problem for maximizing energy production (Sect. 2.6.1). Irradiance is computed using first principles and the solar coordinates at high temporal resolution from the Meeus (1998) algorithms.  (2017)):

Irradiance on a Tilted Solar Panel
where t is time, emphasizing that S, r, and β are functions of time. The angle β represents the solar zenith angle in a coordinate system for which the x-y plane is the solar panel plane. The short-term and long-term variability of S o is not considered here,   Panel-Sun geometry and energy production over a given day is given by the function solar_panel_tilt_insolation.m, which computes Eq. 7 at chosen time steps and integrates it with respect to time over a day. Optionally, a 3D plot is produced showing the orientation of the panel, the normal to its plane, and the directions towards the Sun at local solar noon and in 1-hour steps from about sunrise to about sunset (Fig. 7). Also produced are time series plots of solar azimuth and altitude (Fig.   8A), and the angle β (Eq. 7), and instantaneous panel and horizontal irradiances (Fig. 8B). An annual time series of energy 225 production is computed and plotted by the function solar_panel_annual.m (Fig. 9).
The influence of clouds and atmospheric absorption, scattering and refraction is not taken into account, i.e. top-of-theatmosphere (TOA) irradiance is used, and the center of the solar disk is considered. A 100 W nominal solar panel rating of 100 W is used in the greatly simplified energy production calculations (no temperature effects on panel performance are considered, and panel is assumed to produce more than 100 W when irradiance exceeds the standard of 1000 W). Clear view of 230 the horizon in every direction is assumed. Students can be tasked to study the influence of and/or improve these assumptions.
A major improvement that could be given as assignments of various kinds and complexity to more advanced students would be to combine the code presented here with knowledge of local climatology (e.g. clouds, water vapor profiles) and an atmospheric radiative transfer model (e.g. SBDART -Ricchiazzi et al. (1998)) to improve the solar energy production estimates. Results can also be compared to the NREL solar resource maps ( NREL (2020); Sengupta et al. (2018)). Furthermore, the effect of 235 using simplified models of solar declination (rather than the full Meeus (1998) algorithms) (e.g. as given in Clack (2017) as a function of day of year) can be investigated. Finally, a natural development that can use the functions producing Fig. 8 is to task students with computing the times of local solar noon, sunrise and sunset, using various criteria such as the center of the solar disk at the mathematical horizon vs. at a given altitude below the horizon, to account for atmospheric refraction and the angular size of the solar disk.
240 Figure 9. Time series of annual energy production using the defaults as described in the captions to Fig. 7 and Fig. 8 and in the function producing this plot: solar_panel_annual.m. Solstices and equinoxes are plotted as vertical dash-dot red lines.

Solar Panel Orientation Optimization
An important practical problem that needs to be solved in the installation of solar panels is deciding what panel orientation will maximize energy production, e.g. over year for a fixed panel. Using the solutions described in Sect. 2.6, the solar_panel_orientation_optim.m function optimizes the orientation of a fixed-mount solar panel to maximize annual energy production. Specifically, the function finds the optimal value of two variables: 1) the angle between the direction 245 of the normal to the solar plane and the local vertical, i.e. the tilt of the solar panel, and 2) the azimuthal orientation of the solar panel in the horizontal plane, with respect to true North.
A hybrid optimization technique using genetic algorithms and and an interior-point algorithm is used, illustrating the implementation of a generic constrained optimization technique that is useful for a wide variety of problems and in itself is pedagogically valuable. The user needs to have the Optimization© and Global Optimization© toolboxes in MATLAB© for 250 this function to run. As an example, the optimization function returns a tilt of ∼31.35 • and azimuth of within less than 0.01 • of 180 • , using a default latitude of φ = 33 • 7 45 N (and takes ∼8 min to run in parallel (needs the Parallel Computing Tool-box© of MATLAB©) on a 6-core 4 GHz Intel© processor, using a time step of 5 min). This indicates that to maximize energy production over an entire year, the panel should be facing southwards (as expected), and be tilted at approximately the same angle away from the vertical as the latitude, but ∼2 • less, somewhat favoring better exposure during the summer months.

255
Multiple exercises can be assigned to students using this optimization framework -the most straightforward one being to investigate optimal orientation for various latitudes. Students can also investigate the effect of the various assumptions used by improving upon them or changing them, for example by incorporating atmospheric effects (See Sect. 2.6) that are likely to change the optimization results. Students can also be tasked to come up with an analytical gradient-based solution to the optimization (which would be computationally much more efficient), perhaps using simplified models of solar coordinates 260 (e.g. Clack (2017)). Further, students could be tasked to investigate the effects of changing the time step or taking or not taking into account Sun-Earth distance. More advanced students can investigate the effects of climatological atmospheric properties and the spectral response of a typical solar panel with respect to the spectral irradiance at the bottom of the atmosphere, as modeled by an atmospheric radiative transfer model (See Sect.2.6) and perhaps using meteorological Reanalysis products to obtain data on atmospheric properties and climate at a given location. Students can also explore the effect of having several 265 tilt options throughout the year (e.g. summer tilt vs. winter tilt -and when is it best to change the tilt), as well the effect of a tracking solar panel on energy production as compared to a fixed panel. All of these suggestions can be implemented using the existing code as starting framework. Conversely, students for whom delving into the code and computational methods is not appropriate can greatly benefit from the solar panel visualizations described in Sect. 2.6 and just running the code for various latitudes and panel orientations.

Solar Declination, Global Insolation and Daylength Plots
In this section, I present a few plots that are useful pedagogically in Earth science, physical geography, meteorology/climatology, astronomy, and related classes. A global visualization of mean daily top of the atmosphere (TOA) insolation as a function of date and latitude (Fig. 10)A is a frequently used plot that nicely illustrates seasonality on Earth and its variability with latitude, as well as the asymmetry of the hemispheres due to the eccentricity of Earth's orbit (the South Pole area receives higher 275 insolation at the December solstice than the North Pole area at the June solstice) and the temporal phase shift of the seasons between the hemispheres. The figure presented here is enhanced in ways specifically targeting instructional value, e.g.
students are taught to read a false color map of a variable that is a function of both a spatial and a temporal independent variables. Contours are also added to improve readability and can be used to teach isolines. Furthermore, the Equator and the solstices/equinoxes are drawn. This figure applies to contemporary Earth; instructors and students are referred to the Earth 280 Orbit Model v2.1 (Kostadinov and Gilb (2014)) if interested in producing this plot for times in the geologic past or future, e.g. when teaching/studying the Milankovitch cycles and their effect on insolation and climate, glacial/inter-glacial cycles and related paleoclimatology topics.
The data used in Fig. 10A can be used to compute zonal annually averaged insolation and its temporal standard deviation, illustrating seasonality strength (Fig. 10B). Importantly, these data can be numerically integrated to obtain the globally and 285 annually averaged TOA insolation, an essential variable controlling planetary energy balance and climate. This concept has where S o is TSI (as in Eq. 7), A is Earth's albedo, T is Earth's temperature (in Kelvin), and σ is the Steffan-Boltzmann constant. This is arguably the most important equation of climate science and by extension all geosciences (Kasting (2010)).
Note that this equation ignores the greenhouse effect and therefore significantly underestimates Earth's surface temperature, rather estimating Earth's "skin temperature" (Archer (2012)). The left-hand side of Eq. 8 is estimated numerically from the data in Fig. 10B by weighting the zonal averages by the cosine of the latitude, and is also computed directly. The small difference 295 between these two values (Fig. 10B) is useful pedagogically to illustrate the inaccuracies introduced by numerical integration and the influence of choices of a spatio-temporal sampling grid.
where φ is the latitude, and the inverse cosine result should be given in degrees. As previously discussed, Eq. 9 also considers the center of the solar disk being on the mathematical horizon as the sunrise/sunset times, and atmospheric refraction and twilight are not considered. It is derived by setting h = 0 • in Eq. 4. Students should be encouraged to make more sophisticated sunrise/sunset estimates (see Sect. 2.6) or obtain sunrise/sunset times from, e.g., the NOAA Solar Calculator (NOAA (2020)) 305 and compare the results. Fig. 11 also shows the solstices/equinoxes and solar declination, referred to in the plot as the subsolar point, making it suitable for use in introductory physical geography and Earth science classes for non-majors where introduction of additional astronomical terms such as declination may not be desirable. from simply using the provided figures in class to illustrate often challenging concepts, to using the code to give students more advanced physical geography, Earth science, remote sensing/GIS, scientific computing, astronomy, mathematics, and related tasks. The author has used many of the visualizations here in the classroom and in assignments. Students generally respond favorably to them and indicate in course evaluations that visual aids and especially video materials are very helpful for their learning. However, care has to be taken not to introduce plots, terminology and material that is too technical for the particular 320 classroom setting, e.g. a general education class with primarily non-STEM majors. A significant advantage of AstroGeoVis v1.0 is that the code is open source and commented, and is itself of significant instructional value, giving instructors and students opportunities to dynamically generate/modify figures and inspect some of them in 3D with PTZ capability. Importantly, the computational methods can be inspected and instructors and/or students can perform quality control and introduce further developments. To run the functions/scripts of AstroGeoVis v1.0, no installation is required. Simply download the folder containing all mfiles, and make sure your MATLAB© working directory is pointing to that folder, or put it on your path. Run the function/script of interest on the command prompt using the appropriate arguments as described in the code comments (some can be omitted 335 and defaults will be used).
The function date2jd_vec_v2020.m converts a date and time instant to Julian day (JD). The functions sidereal_time.m and solar_coord.m compute the sidereal time and solar right ascension, declination and the distance to the Sun for a given moment in JD. These functions use the Meeus (1998) algorithms. Note that figures produced by AstroGeoVis v1.0 will have a slightly different appearance (fonts, axes labels/titles/legends presence and exact wording, etc.) 340 than those in this manuscript for formatting reasons (to ensure figures fit in one column for compactness). No information or data of substance is different. comparison.
-The timeanddate.com website (timeanddate.com (2020)), contains a plethora of astronomical information and visualizations and related tools and is an excellent teaching resource.
Author contributions. T.S.K. designed the study, wrote the code, and wrote the manuscript.
Competing interests. The author declares no competing interests.

355
Disclaimer. The AstroGeoVis v1.0 software is supplied "as is". No warranty is given, express or implied, of fitness for any purpose. Under no circumstances shall the author or his institution be liable to anyone for direct, indirect, incidental, consequential, special, exemplary, or any other kind of damages (however caused and on any theory of liability, and including damages incurred by third parties), arising from or relating to this software, or user's use, inability to use, or misuse of the software, or errors of the software. This software is not guaranteed to be error-free and is not meant to be used in any mission-critical applications. Use at your own risk and verify with other sources when 360 appropriate.
Acknowledgements. The author thanks California State University San Marcos (CSUSM) for providing support for this research. He thanks the University of Richmond and the Andrew W. Mellon Foundation for providing support for development of earlier versions of some parts of this work. He also thanks the students of University of Richmond Fall 2012 GEOG 250 class for collecting data used here as an example.
Further, he thanks his CSUSM students in the 2018-2020 period for providing feedback on the use of some of these visualizations in the