[Maxima-discuss] proof of concept Go to Maxima translator applied to astronomical formulas

[Robert D. <rob...@gm...>]

I was looking at some stuff about determining the position of the sun
the sky (in order to tell time without a clock, just for fun) and it
seems like "Astronomical Algorithms" by Jean Meeus is a commonly-used
reference. There is an implementation in Go:
https://github.com/soniakeys/meeus

Go, as you may know, is an Algol-family language, apparently intended
to be somewhat more predictable than C (in the sense of enforcing more
constraints on the programmer). Go is a not-very-complex language for
not-very-complex programming -- the lead developer has said as much --
but anyway it turns out it's not that hard to translate to Maxima.

I have a proof of concept of a Go to Maxima translator, see the folder
robert-dodier/maxima_from_go at:
https://github.com/maxima-project-on-github/maxima-packages

The translator is not too strongly organized, having been accreted by
adding rules one after another in the best Maxima development style.
It works well enough to translate a few of the subunits of the meeus
project.

The starting point for the translator is a syntax tree produced by a
tool named asty (https://github.com/asty-org/asty).Then successive
rules are applied to translate specific constructs to a Maxima
equivalent.

The result is not too far from idiomatic Maxima code, although
translating Go's 'return' and 'continue' requires some kind of
nonlocal goto -- I have implemented that with throw/catch in Maxima.
Maybe there are other ways to do it. Also, one could try to clean up
functions which don't have any throw inside them, via some additional
rules.

I wrote some stuff by hand to fill out stuff around the edges --
providing implementations of built-in functions and what-not. Probably
some of that stuff could be replaced by automatically translated code
at some point.

I can say more about any of this if there's interest.

Robert

PS. Here's a first example. I'll call a function
HorizontalByDeclination which takes arguments φ (latitude) and a
(height of sundial gnomon). I am just north of 45 degrees north
latitude, so φ = pi/4 approximately.

Input:

linel: 72;
fpprintprec: 6;
foo: HorizontalByDeclination (math@Pi/4, 1);
displayfoo (l) := (display (l@Declination), display (l@Description),
display (l@Points));
[second (foo), third (foo)];
for x in first (foo) do displayfoo (x);

Output:

                     l@Declination = Angle(- 23.44)

                    l@Description = winter solstice

l@Points = [PointWithHour(H = 8, X = - 18.436, Y = 14.0529),
PointWithHour(H = 9, X = - 3.65579, Y = 4.17006),
PointWithHour(H = 10, X = - 1.63509, Y = 3.00513),
PointWithHour(H = 11, X = - 0.687555, Y = 2.62886),
PointWithHour(H = 12, X = 0, Y = 2.53087),
PointWithHour(H = 13, X = 0.687555, Y = 2.62886),
PointWithHour(H = 14, X = 1.63509, Y = 3.00513),
PointWithHour(H = 15, X = 3.65579, Y = 4.17006),
PointWithHour(H = 16, X = 18.436, Y = 14.0529)]

                     l@Declination = Angle(- 20.15)

               l@Description = winter solstice + 1 month

l@Points = [PointWithHour(H = 8, X = - 9.20431, Y = 6.51529),
PointWithHour(H = 9, X = - 2.93972, Y = 3.15739),
PointWithHour(H = 10, X = - 1.4168, Y = 2.47043),
PointWithHour(H = 11, X = - 0.611073, Y = 2.22519),
PointWithHour(H = 12, X = 0, Y = 2.15925),
PointWithHour(H = 13, X = 0.611073, Y = 2.22519),
PointWithHour(H = 14, X = 1.4168, Y = 2.47043),
PointWithHour(H = 15, X = 2.93972, Y = 3.15739),
PointWithHour(H = 16, X = 9.20431, Y = 6.51529)]

                     l@Declination = Angle(- 11.47)

                   l@Description = equinox - 1 month

l@Points = [PointWithHour(H = 7, X = - 24.4317, Y = 8.25809),
PointWithHour(H = 8, X = - 4.12243, Y = 2.36595),
PointWithHour(H = 9, X = - 1.98334, Y = 1.80487),
PointWithHour(H = 10, X = - 1.06634, Y = 1.61198),
PointWithHour(H = 11, X = - 0.479707, Y = 1.53185),
PointWithHour(H = 12, X = 0, Y = 1.50912),
PointWithHour(H = 13, X = 0.479707, Y = 1.53185),
PointWithHour(H = 14, X = 1.06634, Y = 1.61198),
PointWithHour(H = 15, X = 1.98334, Y = 1.80487),
PointWithHour(H = 16, X = 4.12243, Y = 2.36595),
PointWithHour(H = 17, X = 24.4317, Y = 8.25809)]

                        l@Declination = Angle(0)

                        l@Description = equinox

l@Points = [PointWithHour(H = 6, X = - 2.30959e+16, Y = 1.0),
PointWithHour(H = 7, X = - 5.27792, Y = 1.0),
PointWithHour(H = 8, X = - 2.44949, Y = 1.0),
PointWithHour(H = 9, X = - 1.41421, Y = 1.0),
PointWithHour(H = 10, X = - 0.816497, Y = 1.0),
PointWithHour(H = 11, X = - 0.378937, Y = 1.0),
PointWithHour(H = 12, X = 0, Y = 1.0),
PointWithHour(H = 13, X = 0.378937, Y = 1.0),
PointWithHour(H = 14, X = 0.816497, Y = 1.0),
PointWithHour(H = 15, X = 1.41421, Y = 1.0),
PointWithHour(H = 16, X = 2.44949, Y = 1.0),
PointWithHour(H = 17, X = 5.27792, Y = 1.0),
PointWithHour(H = 18, X = 2.30959e+16, Y = 1.0)]

                      l@Declination = Angle(11.47)

                   l@Description = equinox + 1 month

l@Points = [PointWithHour(H = 6, X = - 6.96976, Y = - 1.0),
PointWithHour(H = 7, X = - 2.95852, Y = 0.121093),
PointWithHour(H = 8, X = - 1.7424, Y = 0.422663),
PointWithHour(H = 9, X = - 1.09888, Y = 0.554057),
PointWithHour(H = 10, X = - 0.661507, Y = 0.620356),
PointWithHour(H = 11, X = - 0.313155, Y = 0.652804),
PointWithHour(H = 12, X = 0, Y = 0.662639),
PointWithHour(H = 13, X = 0.313155, Y = 0.652804),
PointWithHour(H = 14, X = 0.661507, Y = 0.620356),
PointWithHour(H = 15, X = 1.09888, Y = 0.554057),
PointWithHour(H = 16, X = 1.7424, Y = 0.422663),
PointWithHour(H = 17, X = 2.95852, Y = 0.121093),
PointWithHour(H = 18, X = 6.96976, Y = - 1.0)]

                      l@Declination = Angle(20.15)

               l@Description = summer solstice - 1 month

l@Points = [PointWithHour(H = 5, X = - 12.6345, Y = - 5.78768),
PointWithHour(H = 6, X = - 3.8541, Y = - 1.0),
PointWithHour(H = 7, X = - 2.183, Y = - 0.172781),
PointWithHour(H = 8, X = - 1.41273, Y = 0.153485),
PointWithHour(H = 9, X = - 0.93106, Y = 0.316718),
PointWithHour(H = 10, X = - 0.573502, Y = 0.404787),
PointWithHour(H = 11, X = - 0.274616, Y = 0.449399),
PointWithHour(H = 12, X = 0, Y = 0.463124),
PointWithHour(H = 13, X = 0.274616, Y = 0.449399),
PointWithHour(H = 14, X = 0.573502, Y = 0.404787),
PointWithHour(H = 15, X = 0.93106, Y = 0.316718),
PointWithHour(H = 16, X = 1.41273, Y = 0.153485),
PointWithHour(H = 17, X = 2.183, Y = - 0.172781),
PointWithHour(H = 18, X = 3.8541, Y = - 1.0),
PointWithHour(H = 19, X = 12.6345, Y = - 5.78768)]

                      l@Declination = Angle(23.44)

                    l@Description = summer solstice

l@Points = [PointWithHour(H = 5, X = - 7.81708, Y = - 3.96219),
PointWithHour(H = 6, X = - 3.26181, Y = - 1.0),
PointWithHour(H = 7, X = - 1.97292, Y = - 0.252386),
PointWithHour(H = 8, X = - 1.3119, Y = 0.0711595),
PointWithHour(H = 9, X = - 0.876674, Y = 0.239805),
PointWithHour(H = 10, X = - 0.544099, Y = 0.332764),
PointWithHour(H = 11, X = - 0.261541, Y = 0.380393),
PointWithHour(H = 12, X = 0, Y = 0.395121),
PointWithHour(H = 13, X = 0.261541, Y = 0.380393),
PointWithHour(H = 14, X = 0.544099, Y = 0.332764),
PointWithHour(H = 15, X = 0.876674, Y = 0.239805),
PointWithHour(H = 16, X = 1.3119, Y = 0.0711595),
PointWithHour(H = 17, X = 1.97292, Y = - 0.252386),
PointWithHour(H = 18, X = 3.26181, Y = - 1.0),
PointWithHour(H = 19, X = 7.81708, Y = - 3.96219)]