Most of the programs I make at work are not suitable or interesting to post here, but I have also made some reusable libraries, and some curious side projects.
I'll start with one of the latter, made in
Maple. It calculates the equations of the position of the Sun, relative to any point on the Earth, as a function of the time of the year.
Code:
# speed of Earth's orbit in rad/h
oo_t:= evalf(2*Pi/(365.25*24)):
# speed of Earth's rotation in rad/h
oo_r:= evalf(2*Pi/24):
# tilt of Earth's rotation axis
tilt:= evalf(23.44*Pi/180):
# latitude: 41º 38' N
lat:= evalf(( 41 + 38 /60)*Pi/180):
# longitude: 4º 43' W
lon:= evalf(( -4 - 43 /60)*Pi/180):
x:= evalf(Vector(3, [cos(oo_t*t), -sin(oo_t*t), 0])):
# changes of basis
B:= simplify(evalf(Matrix(3, 3, [
cos(tilt), 0, sin(tilt),
0, 1, 0,
-sin(tilt), 0, cos(tilt)
]).Matrix(3, 3, [
cos(oo_r*t), -sin(oo_r*t), 0,
sin(oo_r*t), cos(oo_r*t), 0,
0, 0, 1
]).Matrix(3, 3, [
cos(lon), -sin(lon), 0,
sin(lon), cos(lon), 0,
0, 0, 1
]).Matrix(3, 3, [
sin(lat), 0, cos(lat),
0, 1, 0,
-cos(lat), 0, sin(lat)
]))^%T):
# finally, parametric equation of position
r:= simplify(B.x):
The algebra behind it is sketched here:
(click)
The purpose is to know the position of the sun at any hour of any day of the year, but I think the most visual result is had by plotting the Sun's trajectory along a whole year.
Code:
plots[spacecurve]( r, t=0..365.25*24, numpoints=87660, axes=normal, labels=["S-N", "E-W", "z-n"], tickmarks=[0, 0, 0], shading=zhue, thickness=2, orientation=[130, 60]);
For example in Valladolid, Spain:
(click)
One axis goes from South to North, another from East to West, and the other is the verticality. The colours colder than blueish green are below the horizon--it's night.
Another example: the Equator, you can see why there's no Summer or Winter over there:
(click)
And Murmansk, Russia, above the Arctic Circle. Notice how the Sun doesn't set during some Summer days, and it doesn't dawn during some Winter ones.
(click)