Study of the determination of the topocentric lunar librations by a best fit estimate of the plate constants from digital images with a small telescope

Study of the determination of the topocentric lunar librations by a best fit estimate of the plate constants from digital images with a small telescope Fabrizio Pinto EKOSPACE and Department of Aerospace Engineering, Izmir University of Economics Teleferik Mahallesi, Sakarya Cd. No: 156, 35330 Balçova/İzmir, Türkiye. fabrizio.pin[email protected].tr Abstract The determination of the lunar optical librations in longitude and latitude from direct observation of the Moon represents an exciting, though possibly challenging, student project with a profound pedagogical potential. In principle, it is possible to compute the optical librations from the small number of observations needed to obtain the plate constants. In practice, however, the plate constants are more effectively obtained by a best fit of the positions of a larger number of lunar features on a digital image by means of the well-known relations to the selenographic coordinates available from a lunar atlas. Firstly, here we report on an analysis of the relationship between the final accuracy in the determination of the librations by this methodology and the several factors that affect the result. For this purpose, the fundamental photogrammetric equations are studied by means of simulations to investigate the effect of the number of features used, their position on the lunar face, and the resolving power of the optical system. Secondly, actual measurements are carried out on digital images acquired by the author by means of a Canon 500D at the prime focus of a Meade ETX-90. These measurements are also further augmented by the estimates of a class of university sophomore students working independently in an Introduction to Remote Sensing class. Additional refinements are considered for future work, such as parallactic effects, the diurnal parallax and the evolution of librations during the lunar orbit. 1. Introduction The optical librations of the Moon are, technically speaking, sufficiently large to be visible by the unaided eye. Since multiple observations separated in time are needed, historical research has focused on the possibility that keen sky observers and gifted artists, including Leonardo himself, may be among those who first recorded sufficiently detailed images of the lunar disk, although likely without recognizing this phenomenon (Tucci 2022). After the introduction of the astronomical telescope, reports about librations both in latitude and in longitude rapidly appeared (Włodarczyk 2011), commencing with an early handwritten note by Thomas Harriott (c. 1560-1621) in December 1611, followed by the observations by Florentius van Langren (1600–75), Galileo Galilei (1564–1642), and, of course, by Johannes Hevelius (1611-1687). The pedagogical potential of the observation and measurement of lunar librations in astronomy classes is now widely recognized and full discussions of possible projects are available (Buchheim 2015). This topic has been a long time interest for the present author (Hughes 1993), who has devoted efforts to include it in the curriculum since the dawn of the commercial digital camera market (Pinto 1995). It is immediate to develop an awareness of the phenomenon of librations by viewing images of the Moon taken at different times, as dramatically shown in a very popular video developed at NASA Goddard (2011). However, the goal of carrying out reliable measurements of librations from individual images of the Moon taken with small telescopes offers some challenges of a technical and mathematical nature well worth the time of astronomy and space engineering undergraduate students. Geometrical considerations show that the positions on a digital image of 4 lunar features of known selenographic coordinates are sufficient, in principle, to derive the information needed to compute the librations. However, as we shall see, attempts to naively follow this minimalist strategy lead to unsatisfactory results, due to the limited accuracy available when using ordinary, small telescope systems. That this should be of no surprise 2 Figure 1. The geometry of the measurement of the position of a point P upon a camera sensor. In addition to an obvious scale factor, the sensor coordinates (X, Y) must be connected by transformation laws to the standard selenographic Cartesian coordinates, (, , ). This involves all angles that cause departures from the idealized geometry of this Figure, that is, librations in latitude and longitude, and rotation of the field of view (see text). Notice the definition of longitude P, positive eastward, towards Mare Crisium, shown by the generic angle , and of the latitude P, positive towards the hemisphere containing Mare Imbrium (not indicated above to avoid clutter) (see Kopal (1966), p. 169). can be best appreciated by considering that such a famous observer as Saunder (1900, 1901, and 1905) reported continued challenges measuring the selenographic coordinates of lunar features, as difficult as detecting the parallaxes of stars. Saunder was using the most remarkable of the French Coudé equatorials, the Grand Coudé reflector at Paris, an impressive instrument with a 60 cm objective and a focal length of 18 m. The fascinating article by Lequeux (2011) includes historic photographs of those instruments, including of the sophisticated controls at the focal point. The images of the Moon obtained by such Coudé reflectors were so exquisite as to be produced at 80 cm of diameter “on special plates supplied for free by the Lumière brothers,” and they were “still found to be useful when choosing the landing sites for the Apollo project” (Lequeux 2011). In this paper, we analyze such challenges and focus on a specific issue within the broader goal of carrying out student lunar libration measurement projects. This is the error analysis that must accompany any final estimate. In order to extract more reliable information from available digital images, a statistical approach based on the measurement of the coordinates of a much larger number of surface features has been developed. By a best fit, the Plate Constants can be estimated, leading to values of the librations, accompanied by an understanding of the sources of error. One further element that has enabled bringing such activities within reach of students in undergraduate classes has been the development of computer algebra systems (CAS), such as that available in the Mathematica language. Indeed, personal discovery of Mathematica by this author occurred specifically during the search for a simple approach to sharing data analysis of potentially hundreds of coordinates from a single image for libration calculations. The inherent power of the Mathematica language (Wolfram 2013), along with the availability of the notebook platform, also introduced by Steve Wolfram, have represented personal ‘game changers’ in the pedagogical treatment of this problem. In recent years, the Free Wolfram Engine has become available for use on the highly flexible Jupyter Notebooks, a development of historic importance that has further immensely increased the pedagogical potential of this approach (Romer 2018, Somer 2018, Rule 2019).. In the spirit of sharing the present approach via such powerful, modern tools, a Jupyter Notebook was prepared and is available as Supplementary Materials, with LaTeX-formatted equations, worked examples, and plots (Pinto 2025). This Supplementary Materials file, stored in the Open Science Framework, is written in Jupyter Notebook format (ipynb), fully executable and editable upon installation of the Free Wolfram Engine; it can be also downloaded in html language (not executable) for simple viewing. The workflow leading to the image provided and to the measurements, including other software used, is briefly described below. However, our emphasis is on the analysis of the sensitivity of libration calculations to various important factors, explored by using synthetic data; later, one sample image of the Moon is used along with data produced by almost 500 student measurements, which are analyzed by means of the algorithms presented in the Supplementary Materials. 2. Theoretical background Following the standard approach (Saunder 1900, 1901, and 1905; Kopal 1966, Kopal and Carder 1974), the fundamental equations defining the selenographic Cartesian coordinates of a point P on the lunar surface) in terms of the selenographic longitude and latitude, (P, P) are: P = cos P sin P , (1)  P = sin P , P = cos P cos P, P 2 +  P2 + P2 = 1, assuming that the Moon is spherical. With reference to Fig. 1, it is clear that he origin of the (P, P) reference frame would be represented by a point at (X=0, Y=0) on the sensor in the absence of librations and field rotation (we define the field rotation angle as the angle made by the Y-axis with the polar axis of the Moon, that is, the departure from an ideal North-up orientation). If the disk of the Moon is not centered on the origin of the sensor reference frame, further terms must be added. In general, apart from a scale factor equal to the ratio of the physical radius of the Moon to its radius on the sensor, the above equations simply imply that, if the Moon is centered on the sensor, the coordinates of the image of point P on the sensor, point, Q, are connected to its selenographic Cartesian coordinates as: XQ = P , (2) YQ =  P . In the presence of librations (lib, lib), however, these equations become: XQ = P cos lib - P sin lib , (3) YQ = - P sin lib sin lib +  P cos lib - P cos lib sin lib. Finally, if the field is rotated by an angle rot during the exposure, the above equations must be further modified by the standard rotation transformation, R(rot), around the -axis. X’Q = XQ cos rot - YQ sin rot , (4) Y’Q = XQ sin rot +YQ cos rot . By substituting Eqs. (3) into Eqs. (4), we obtain the full result (neglecting terms that describe noncentering). For succinctness, here we do not reproduce the full expression (Supplementary Materials, Sec. 2.1). The conclusion to be drawn from this treatment is that Eqs. (2) must be completely generalized to consider the coordinates of the sensor as linear functions of all selenographic Cartesian coordinates. Including constants that will also describe both the scale factor and non-centering, we have: X’Q = A P + BP + C P + D, (5) Y’Q = E P + FP + G P + H. where (A, …. H) are referred to as the Plate Constants. As apparent from the above discussion and proven numerically in the Supplementary Materials, a minimum of 4 measurements is needed to determine each of the two sets of constants, which allow for the computation of the optical lunar librations. The Plate Constants can then be used, for instance, to compute the radius of the Moon on the sensor, in pixels, by using the two equivalent equations: RMoon, X = (A2 + B2 + C2)1/2, (6) RMoon, Y = (E2 + F2 + G2)1/2, RMoon, X = RMoon, Y . Also, the center of the lunar disk in sensor reference frame in the absence of librations and field rotation is given by: XCenter = D , (7) YCenter = H . Let us now proceed to extract the desired information from the Plate Constants. From Eqs. (3), we see that the coordinate  P appears only in the expression for YQ. Therefore, by dividing the Plate Constants of the coordinate  P in Eqs. (5), B and F, by each other and with a sign change due to the negative sign at the righthand-side of the former of Eqs. (4), we can eliminate rot: rot = - arctan 𝐵 𝐹 (8) Since we now have the rotation angle, we can apply the inverse transformation to both Eqs. (5), so as to return to the unrotated expressions at Eq. (3): R-1(rot)(𝑋’ 𝑌’) = R(-rot) (𝑋’ 𝑌’) . (9) The longitude libration can then be obtained from the observations in various, equivalent ways by recognizing that, for instance: lib, OBS = - arcsin C’ , (10) where C’ is the -coordinate coefficient after the rotational back-transformation. Since this libration only depends on the P–coefficient, unnecessary complications may be avoided by letting the R(-rot) operator act only upon a column vector, (C, G). Once lib, OBS is known, for the libration in latitude we have: lib, OBS = - arcsin (G’/cos lib, OBS), (11) 4 Figure 2. The complete optical system (see text). where 𝐶′ and 𝐺′ are the ‘unrotated’ coefficients: (𝐶’ 𝐺’) = R(-rot) (𝐶 𝐺). (12) 3. The hardware In order to create applicable synthetic data, the following information regarding the hardware available for actual image acquisition is provided. The telescope used is a Meade ETC-90  a MaksutovCassegrain reflector with a nominal aperture listed as D = 90 mm (3.5”) and focal length FL = 1,250 mm, corresponding to an f/13.8 aperture (for technical specs, following the demise of Meade Instruments, see, for instance, Meade Instruments, Meade ETX Telescopes, and OpticsStar, 2025, in the References). Imaging is carried out by a Canon 500D camera body (Pinto 2024, and References therein) placed at the prime focus by a TelescopeAdapters ring mount (https://www.telescopeadapters.com/ ) and powered by an external power supply. The system is assembled on an equatorial mount but no tracking is used in this case (Fig. 2). The camera sensor is a CMOS Type APS-C (14 bit) with a maximum resolution of (4,752  3,168) pixels and a pixel pitch of spix = 4.68 m. Focusing is achieved by visually monitoring the image on the back LCD screen magnified electronically at up to 10 times, while providing adjustments by means of the Meade electric focuser (Fig.3). Figure 3. LCD screen at 10 magnification (see text). By using the above data and using an average wavelength equal to light = 500 nm, the theoretical resolving power of the system (Rayleigh criterion) is: res  1.22 light/D = 1.40”. (13) 4. Synthetic data analysis In this Section, a few examples are discussed in order to estimate the effect of uncertainties on the final results by generating synthetic data for assumed field rotation and libration angles and by monitoring the accuracy of recovery of those assumptions. As mentioned in the Introduction, synthetic data analysis was carried out by means of the Wolfram Mathematica language. The Free Wolfram Engine (WFE), v. 11.3, was installed on Jupyter Notebook, server v. 6.5.4, running on a relatively old but fully operational hp Pavilion laptop under Windows 8.1. Commencing in 2014, the Mathematica language has featured the powerful GeoGraphics function, which can potentially aid in carrying out synthetic analyses using actual lunar terrain (see Supplementary Materials for examples). However, here we shall not associate synthetic data to existing formations but shall instead generate the selenographic coordinates by using a random number generator. For this purpose, the equations of Sec. 2 were coded in Mathematica and used to produce simulated measurements, (X’i, Y’i) corresponding to hypothetical surface features of selenographic coordinates (i, i), with i = 1 … Nsynth. The arbitrary inputs of the simulation are the total number Nsynth of such entries, the ranges in selenographic longitude and latitude Figure 4. Randomly extracted positions on the lunar disk for Nsynth = 4. The axes represent the sensor reference frame in normalized coordinates. (min  i  max, min  i  max) in which measurements are assumed to have been taken, the field rotation angle, rot, and the librations, (lib , lib ). The selenographic coordinates (i, i), are randomly drawn from uniform distributions given by the desired coordinate ranges. The corresponding values (X’i, Y’i), of the measured coordinates, are computed as we have previously discussed. In order to explore the effect of measurement uncertainty, noise is injected into such measurement set by adding random errors distributed according to a normal distribution (Taylor 1997) of expected value  = 0 and standard deviation expressed in terms of a dimensionless multiplicative parameter of the theoretical resolution of the optical system found at Eq. (13), that is,  =  res. After computing the pixel coordinates (X’i, Y’i), an arbitrary constant is added to each, representing the off-center position of the lunar disk. In order to connect this simulation to the optical system used, the radius of the Moon on the sensor was obtained from the distance of the Moon at the epoch of the observation to be analyzed in the following Sections (2020-04-30 22:54:04). The ephemerides used are those of the Virtual Moon Atlas, which gives a distance rEM = 3.72939506 105 m. Since the physical radius of the Moon, if assumed spherical, is RM = 1.7373 103 km, the radius in pixels is computed to be: Rpix = (RM /rEM ) (FL /spix) = 1469.77 pix. (14) It is important to consider that, in practice, departures of the actual focal length and pixel pitch from their nominal values do exist and a study must be carried out to determine the effective focal length of the instrument at the time of observation, which should be considered a time variable even during one session (as indeed is repeatedly mentioned by Saunder). This fact can be captured by another arbitrary factor to multiply the result of the previous equation and to be obtained from the measurements. In what follows, we briefly survey the results. Full algorithms and plots are in the Supplementary Materials. 4.1 Exact measurements In this first example, we consider the ideal case of 4 exact measurements ( = 0) and we verify that, in this case, we recover the assumed arbitrary parameters. Therefore, Nsynth = 4, lib = - 7.4 deg, lib = -3.7333 deg, rot = 21.7048 deg. The longitudes and latitudes are assumed to be -40  i  40, -40  i  40. The off-center constants are Xc = 1250.6 and Yc = 2538.2. By using a random seed (SeedRandom) equal to 1234, for repeatability, we find, for instance (these values are rounded below for practicality, all angles are in degrees): X1= 2037.4, Y1= 2981.5, 1 = 30.13, 1 = 1.76 X2= 1489.6, Y2= 2688.9, 2 = 3.500, 2 = 1.65 X3= 1027.0, Y3= 3345.3, 3 =- 3.278, 3 = 30.78 X4= 2141.1, Y4= 3156.6, 4 = 39.98, 4 = 7.035. The location of the randomly extracted points is shown at Fıg. 4. Upon execution of the algorithm to estimate the field rotation angle and the librations, all values are recovered with absolute numerical errors 10-13–10-12. Also the condition that the two radii estimated from Eqs. (6) be equal to each other is satisfied within an absolute value of  10-11. The remarkable stability of the process can be appreciated by repeating the computation for 4 points clustered in s a very small area of the lunar disk. For instance, even by choosing coordinates bound as 70 deg  i  75 deg, 80  i  85)  to be considered extreme from the observational point of view  the Mathematica algorithm recovers the assigned parameters with the same accuracy. 4.2 Random errors In this Section, we explore the effect of injecting random errors into the measurements. For this purpose, we shall repeat the computations of previous 6 Figure 5. Showing the log-log plot of the absolute value of the error  in the libration in longitude as a function of the parameter . Measurements spanning 70 deg  i  75 deg, 80  i  85 (black) lead to larger errors than those in a much wider area, -40  i  40, -40  i  40 (blue). cases but assume that   0. The cases considered are shown in Fig. 5 in a few illustrative examples. As can be seen, the choice of measurement points within a smaller area leads to numerical errors orders of magnitude larger than in the case of a narrow area. Interestingly, however, if the number of points employed is much larger, the error again decreases in any case, as intuitively expected. For instance, by employing Nsynth = 100 points, the libration in longitude error in the case of a wide survey area (blue curve) for  = 10-3 decreases from  2.7  10-3 deg to 5.6 10-4 deg. This approximately corresponds to a factor of 25 times, equal to the ratios of points employed, 100/4. These numerical estimates do not rise to the level of a statistics theorem but offer guidance as to whether the measurement we seek to carry out is possible with the smaller instrumentation available. From the practical point of view, choosing such small values of  is clearly unphysical but, in our case, it helps to present the dramatic behavior of the error for extremely small survey areas. As one last example, we consider the case of a large sample (Nsynth = 497 points) the author secured over time by presenting the same image to students enrolled in an Introduction to Remote Sensing class. By choosing parameters appropriate to the first quarter Moon (0 deg  i  90 deg, -90  i  +90), we find, for a physically permissible, though quite optimistic, parameter  = 1.0, errors |lib|  0.012 deg, |lib|  0.011 deg, respectively. Even for much higher values of  = 102, the errors over such a wide survey area as half of the lunar face remain as low as  0.05 deg in both longitude and latitude. As we shall see, this is incompatible with the results from student measurements, which yielded errors  0.6 deg in both longitude and latitude. This is not physically consistent with the optical system available, as it Figure 6. Analogously to Fig. 4, showing randomly extracted positions on the lunar disk for Nsynth = 497 with (0 deg  i  35 deg, -60  i  +60). The axes represent the sensor reference frame in normalized coordinates. would require unrealistically high values of the  parameter. On the other hand, for  = 20, and considering a restricted area defined by (0 deg  i  35 deg, -60  i  +60), leads to synthetic samples that produce the observed errors. In this Section, for brevity, we have focused on the libration in longitude. However, our conclusions also apply to the libration in latitude, to field rotation, and to the similarity of the radii along the two sensor axes. Such data can be obtained by accessing the Supplementary Materials. 5. Image Acquisition and processing For simplicity of execution, unlike deeper sky objects (Pinto 2024), no calibration frames were taken. All images were saved to the proprietary Canon RAW CR2 format. This project was carried out both with single images and with much higher quality images derived by stacking a large number of frames. No attempt has yet been made to determine whether processing of lunar images leads to improved estimates of the librations. Although that appears reasonable, photogrammetry also requires experience on the side of the observer. In this case, students without previous astronomical experience open a FITS file, or create one by image stacking programs, and carry out the measurements described by using ASTAP and the Visual Moon Atlas as a guide. Naturally, particular attention is drawn to the terminator area whereas areas where the Sun is high on the horizon on the Moon are neglected. Figure 7. The image of the Moon analyzed in this work. The image considered herein is part of a set acquired on 2020-04-30:19:54:04 (UTC) from the author’s observatory in the Izmir region. The exposure time was 1/20 s at 200 ISO, while the Moon was approximately at an altitude of 43o. 6. Photogrammetric image analysis The image shown at Fig. 7 was provided to students in FITS format and analyzed by them by means of the ASTAP program and the Virtual Moon Atlas. This allowed them to build a file equivalent to that created synthetically by Mathematica and discussed earlier. Every student worked on approximately 20-100 features but it was noticed that some students shared choices of the same features, being predictably attracted to clearly identifiable formations near the terminator, as shown at Fig. 8. All files were then merged without particular curation and analyzed by the same algorithm used to produce the simulations presented in Sec. 4. The result found is: rot = 21.7048 deg, lib = - 7.99 deg, lib = -4.40 deg. The errors on the librations, estimated from the model, are  0.6 deg. Importantly, the statistical analysis carried out by Mathematica also shows that the confidence level intervals of the plate constants obtained from the model are of the same order of magnitude as those from the observed data (1-10). Figure 8. The lunar formations chosen by the student analysts for libration determination. Notice the much higher density of entries along the terminator (this is an unrotated image, North up, East to the right). 7. Conclusions Possibly the most important quantitative conclusion of this paper is that the crucial factor to determine the accuracy of the libration estimates is the distribution of the features to be measured. The synthetic model shows that a large number of measurements can indeed lead to very accurate results. However, such potential disappears if parts of the lunar disk are not properly sampled. In the future, this author intends to further explore these conclusions in order to achieve detection of the diurnal topocentric librations and of the physical librations of the Moon. An additional improvement of the model, already being tested, is the inclusion of the finite distance of the Moon in the data analysis. The pedagogical potential of this activity is being exploited by transitioning interested students to mapping of small bodies of the solar system, such as asteroids and moonlets, for the purpose of modeling their gravitational field. 3. Acknowledgements I thank my astronomy students at Boise State University who, over three decades ago, attempted the first steps of this project, particularly Jennifer Sype, Debbie Street, and Carolyn Hughes (Hughes, 1993). I 8 am grateful to all students in my Introduction to Remote Sensing class (AE-202) at the Izmir University of Economics in the years 2017-2025, mostly sophomores in our space engineering program, and to some in my Astrophysical Systems class (AE415), mostly seniors, for their work on digital lunar photogrammetry. My heartfelt gratitude also goes to all programmers who developed and made available the remarkable software used in this work under various licenses. 4. References Buchheim, R. K., Astronomical Discoveries You Can Make, Too! (Springer, Cham, Switzerland, 2015). Ch.2, Project 14. Hughes C. L. “Photogrammetric determination and interpretation of geocentric and topocentric lunar librations,” (1993). Announcer: Joint APS/AAPT April Meeting Highlights, Boise, ID, 23, 45 (F. Pinto, Supervisor). Kopal, Z., An Introduction to the Study of the Moon (Springer-Science+Business Media, B. V., Switzerland AG, 1966). Ch. 3. Kopal Z., and Carder, R. W. Mapping of the Moon (Springer-Science+Business Media, B. V., Dordrecht, 1974); Ch. 3. C. Lagrande and P. Chevalley, Virtual Moon Atlas, Version 8 (2023). Accessed: 3 May 2025. URL: http://ap-i.net/avl/en/documentation Lequeux, J. “The Coudé Equatorials,” (2011). JAHH, 14, 191-202. Meade ETX telescope. Wikipedia. Accessed 4 May 2025. URL: https://en.wikipedia.org/wiki/Meade_ETX_telescope. Meade Instruments. In Wikipedia. Accessed 4 May 2025. URL: https://en.wikipedia.org/wiki/Meade_Instruments. Meeus, J., Mathematical Astronomy Morsels (William Bell, Inc., Richmond, VA, 1997), Ch. 6. NASA Goddard, “Moon Phase and Libration,” uploaded by NASA, 15 June 2011. Accessed: 2 May 2025. URL: https://www.youtube.com/watch?v=3f_21N3wcX8 OpticsStar, “Meade ETX90 Observer (ETX-90),” (2025). Accessed: 4 May 2025. URL: https://www.meadeuk.com/Meade-ETX90Observer.html Pinto, F. “CCDs in the mechanics lab - A competitive alternative? (Part I),” (1995). The Physics Teacher, 33, 436-441. Pinto, F.. “Astrometry with a DSLR camera. Observation planning, image acquisition, and data analysis,” 2024. In Proceedings for the 43rd Annual Conference of the Society for Astronomical Sciences, Ontario, CA, pp. 29-38. URL: https://socastrosci.org/ https://socastrosci.org/wpcontent/uploads/2024/06/2024Proceedings_Ver1.3c.pdf Pinto, F., “Study of the determination of the topocentric lunar librations by a best fit estimate of the plate constants from digital images with a small telescope. Supplementary materials,” (2025). Accessed: 27 April 2025. URL: https://osf.io/jfmyb/?view_only=053faaadffe84190b4 78824503767624. Romer, P. “Jupyter, Mathematica, and the Future of the Research Paper,” 13 April 2018. Accessed: 20 May 2023. URL: https://paulromer.net/jupytermathematica-and-the-future-of-the-research-paper/ Rule, A. “Ten simple rules for writing and sharing computational analyses in Jupyter Notebooks,” (2019). PLoS Comput. Biol., 15, e1007007. Saunder, A. “The determination of selenographic positions and the measurement of lunar photographs,” (1900). MNRAS, 60, 174-201. Saunder, A. “The determination of selenographic positions and the measurement of lunar photographs. [Second Paper.],” (1901). MNRAS, 62, 41-61. Saunder, A. “The determination of selenographic positions and the measurement of lunar photographs. [Fourth Paper.],” (1905). MNRAS, 65, 458-473. Somers, J. “The Scientific Paper is Obsolete,” The Atlantic, 5 April 2018. Taylor, J. R. An Introduction to Error Analysis (University Science Books, Sausalito, California, 1997). Tucci, P. “The Moon’s ashen light and libration in Leonardo and Galileo,” (2022). Quaderni di Storia della Fisica, 26, 21-60. S. Wolfram, “There Was a Time before Mathematica …,” 6 June 2013. Accessed: 2 May 2025. URL: https://blog.stephenwolfram.com/2013/06/there-wasa-time-before-mathematica/ Włodarczyk, J. “Libration of the Moon, Hevelius' theory, and its early reception in England,” (2011). J. Hist. Astron., 42, 495-519.