iconEuler Examples

Astronomical Functions

The astronomical functions in Euler have been programmed by Keith
Burnett. Extensions by R. Grothmann.
>load astro
Astronomical functions taken from 

Jean Meeus - 'Astronomical Algorithms'
Montenbruck and Pfleger - 'Astronomy on the Personal Computer'
Duffett-Smith - 'Practical astronomy with your calculator'

and other sources.

Written by Keith Burnett.
Click on the following link to see a list of functions.

See: astro.e.html

Use day to set the current date, here the January 2nd, 2008, 22:00 UTC
(23:00 local time). The function returns the days since 2000.
>now = day(2008, 1, 2, 22, 00)
2923.41666667
This is the Julian date.
>jday(2008, 1, 2, 22, 00)
2454468.41667
The Greenwich siderial time in degrees.
>gst(now)
71.9185589233
Let us compute the geocentric coordinates of the sun. We are in
Winter, so it is expected that the declination is negative. The
distance is in astronomical units.
>psun = sun(now)
[ 282.914502853  -22.9093960506  0.983280777656 ]
>gmoon = moon(now)
[ 216.062120546  -19.4451038467  405224.725114 ]
You can set a variable like this for your own location (longitude
in °, lattitude in °, altitude in m).
>locIngolstadt
[ 11.433  48.767  400 ]
Now, we add the temperature (in C°) and the barometric pressure
(in HP).
>here = [locIngolstadt,0,1020]
[ 11.433  48.767  400  0  1020 ]
We can then compute the altitude and azimuth of the sun. It is
now in NE and 60° below horizon. The function raltaz corrects
for refraction.
>raltaz(now, here, psun)
[ 322.09664678  -59.8739592678 ]
To compute a graph of the altitudes of the sun, we need a function,
which returns the altitude at any time.
>function altSun (now,loc=here) ...
r=raltaz(now,loc,sun(now));
return r[2];
endfunction
Now we can plot. Note, that we need to use map for day and altSun,
since both functions do not work for vector input, but plot2d
needs that for expressions.
>plot2d("altSunmap(daymap(2008,1,2,x,0))",a=0,b=24):

Astronomical Functions

Sunset and Sunrise

We can compute the sunrise with the bisection method now.
>hmsprint(secant("altSun(day(2008,1,2,x,0))",6,8))
07:09:04
And likewise the sunset.
>hmsprint(secant("altSun(day(2008,1,2,x,0))",14,16))
15:27:12
The highest point of the sun is also easily computed. Note: all
times are in UTC!
>hmsprint(fmax("altSun(day(2008,1,2,x,0))",10,12))
11:18:07
There is a function in astro, which does that automatically. It
computes the next rise after a given day.
>printday(rise("sun",now,here))
2008-01-03 07:08:59
Also the next sunset.
>printday(set("sun",now,here))
2008-01-03 15:28:13
Next we make a general function, which compute the next even "rise"
or "set" after a given date.
>function map computehour (month,year,planet;loc,compute) ...
now=day(year,floor(month),floor((month-floor(month))*30)+1);
r=compute(planet,now,loc);
t=getymdhms(r);
return t[4]+t[5]/60+t[6]/3600;
endfunction
Plot the sunrises over one year. This is a bit slow. Please be
patient!
>plot2d("computehour";2008,"sun",here,"rise",a=1,b=13,c=0,d=24,n=25,adaptive=0);
Add the sunsets over one year.
>plot2d("computehour";2008,"sun",here,"set",a=1,b=13,n=25,add=1,adaptive=0);
Add the times of highest altitude.
>plot2d("computehour";2008,"sun",here,"highest",a=1,b=13,n=25,add=1,adaptive=0);
>title("Sunrise and Sunset over one year"):

Astronomical Functions

>printday(rise("moon",now,here))
2008-01-03 02:34:08
>printday(set("moon",now,here))
2008-01-03 11:32:10

Examples Homepage