Euler Math Toolbox - Examples

iconEuler Examples

Astronomical Functions

by Keith Burnett

The astronomical functions were taken from

and other sources.

Extensions by R. Grothmann.

>load astro;

Click on the following link to see a list of functions.

List of Functions

Use day to set a specific date. The function returns the days since 2000.

>now = day(2008, 1, 2, 22, 00)
2923.41666667

For later use, we save the current year, month and day in month. The time is in UTC (universal time coordinated).

>v=getnow(1), year=v[1]; month=v[2]; dmonth=v[3];
[2014,  8,  1,  9,  2,  58,  165]

This is the current Gregorian day number in UTC.

>now = daynow(1)
5325.87706228

This is the offset in hours relative to utc.

>utcoffset = round(24*(daynow(0)-daynow(1)))
2

This is the Julian date.

>long jday(2008, 1, 2, 22, 00)
2454468.41667
>long jdaynow()
2456870.87706

The Greenwich sidereal time in degrees.

>gst(now)
85.6397470722

Let us compute the geocentric coordinates of the sun. In Winter, the declination is negative. The distance is in astronomical units.

>psun = sun(now)
[131.488,  17.989,  1.01497]
>gmoon = moon(now)
[187.323,  -4.44876,  399321]

You can set a variable like this for your own location (longitude in °, latitude 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. The function raltaz() corrects for refraction.

>raltaz(now, here, psun)
[125.717,  48.5564]

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 the altitude of the sun today depending on the hour.

>plot2d("altSun(day(year,month,dmonth,x,0))",a=0,b=24,n=40, ...
   title=year+"-"+month+"-"+dmonth):

Astronomical Functions

We want to plot the sun at 12 o'clock local time over the year.

>function plotSun12 (year,hour=12,loc=here) ...
 plot2d(none,140,220,0,80);
 loop 1 to 12
   d=day(year,#,1,hour,0);
   r=raltaz(d,loc,sun(d));
   plot2d(r[1],r[2],>points,style="[]",>add);
   label(""+#,r[1],r[2]);
 end;
 endfunction

For summer time, use 13 instead of 12.

>plotSun12(year,12-utcoffset,here):

Astronomical Functions

Sunset and Sunrise

We can compute the sunrise with the bisection method now.

>hmsprint(secant("altSun(day(year,month,dmonth,x,0))",0,12,y=-0.5))
03:44:57

The time marks the point, where the sun is half a degree under the horizon.

And likewise the sunset.

>hmsprint(secant("altSun(day(year,month,dmonth,x,0))",12,24,y=-0.5))
18:55:24

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, which does that automatically. It computes the next rise after a given day.

>printday(rise("sun",now,here))
2014-08-02 03:50:30

Also the next sunset.

>printday(set("sun",now,here))
2014-08-01 18:51:11

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 in UTC"):

Astronomical Functions

Moon rise and set in UTC.

>printday(rise("moon",now,here)), ...
 printday(set("moon",now,here))
2014-08-01 09:21:24
2014-08-01 20:59:51

Plot the Sky

The following plots the sky at the current time

>function plotsky (now,loc=here,r=1.5) ...
 plot2d("cos(x)","sin(x)",xmin=0,xmax=2pi,=r,<grid);
 sp=["sun","moon","venus","mars","jupiter","saturn"];
 loop 1 to cols(sp)
   r=raltaz(now,loc,sp[#](now));
   x=cos(rad(270-r[1]))*((90-r[2])/90);
   y=-sin(rad(270-r[1]))*((90-r[2])/90);
   plot2d(x,y,>points,>add,style="o");
   label(sp[#],x,y,pos="cr");
 end;
 label("S",0,-1.5,pos="cc"); label("N",0,1.5,pos="cc");
 label("E",-1.5,0,pos="cc"); label("W",1.5,0,pos="cc");
 endfunction
>plotsky(now):

Astronomical Functions

Examples