﻿ Euler Math Toolbox - Examples

# Astronomical Functions

by Keith Burnett

The astronomical functions were taken from

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

and other sources.

Extensions by R. Grothmann.

```>load astro;
```

Click on the following link to see a 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; month=v; dmonth=v;
```
```[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;
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):
``` 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));
label(""+#,r,r);
end;
endfunction
```

For summer time, use 13 instead of 12.

```>plotSun12(year,12-utcoffset,here):
``` # 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+t/60+t/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"):
``` 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));
```>plotsky(now): 