When I was looking into a few things like angles of the sun and the moon (and uranus) from the earth, I really got annoyed going to suncalc tools to get a particular angle on this and that date and place. So I thought: Why not build my own sunapp?
And so I started..

And here it is: Solar System Animation App

First of all, it is actually not very complicated. The basic rules for the path of the planets around the sun are based on six parameters for each body:

The primary orbital elements are:

N (Ω)= longitude of the ascending node
i = inclination to the ecliptic (plane of the Earth's orbit)
w (ω) = argument of perihelion or argument of periapsis (which is the same for the sun),
a = semi-major axis, or mean distance from Sun (the sum of the distance at periapsis and apoapsis divided by 2.
e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
M = mean anomaly (0 at perihelion; increases uniformly with time)

Related orbital elements are:

w1 = N + w = longitude of perihelion
L = M + w1 = mean longitude
q = a*(1-e) = perihelion distance
Q = a*(1+e) = aphelion distance
P = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)
T = Epoch_of_M - (M(deg)/360_deg) / P = time of perihelion
v = true anomaly (angle between position and perihelion)
E = eccentric anomaly

That is a lot of difficult terminology, which are explained on https://en.wikipedia.org/wiki/Orbital_elements and I will try to explain them in laymen terms:

The reference plane is the ecliptic, this is the plane of earth's orbit around the sun. In the picture this is the dark grey plane.

The inclination is the angle to this reference plane, and thus the angle the orbit has to the ecliptic.

The sun is off course in the ecliptic, so the inclination is 0. The moon is slightly (5.14^o) tilted in the ecliptic, so the inclination (i) is 5.14^o.

The orbital elements are specified for a specific date and time, i.e. January 1^st 2000 at midnight and then calculated for a specific date.

i = 0.0
N = 0.0
a = 1.00000011 (AU)
e = 0.016709
w = 282.9404 
M = 356.0470

i is zero, as the sun and the earth are (per definition) in the same reference plane:The ecliptic.

N: The longitude of the ascending node (the angle of the direction of the sun and the intersection of the plane of reference and the orbit of the earth is undefined, as the intersection is on the whole plane. We therefore can assign the value of 0 to N.

a is 1,00000011 (and thus nearly 1 Astronomical Unit

e is an indication of the shape of the orbit. for e=0 it is a perfect circle, between 0 and 1 the shape is an ellipse, and for e>1 the orbit becomes a parabola. At 1-1-2000 the shape of the orbit gets a 0.016709, so

it is close to a circle. This shape changes over (a very long) time.

So we are left with only 2 interesting values w, and M.

w is the argument of perihelion. The perihelion is the position of the earth in it's orbit when it is closest to the sun. In 2019 this happened on January 3 at 05:20 UTC, when the sun was at 147.099.761 km distance). The argument of the perihelion is the angle between the sun, at 1-1-2000, to the perihelion.

M is the mean anomaly. This is the projection of the angles on a true circle from the ellipsis. This is necessary, because the speed of the earth around the sun changes during the year., It is fastest when the sun is furthest away, and slowest at perihelion. Projected on a circle, the speed never changes, and therefore it is easier to calculate with. The mean anomaly is 356.0470^o.

So these are the orbital elements explained, for 1-1-2000. Now we have to calculate them for a specific date.
For that we are going to use these formulas, where d is the number of days since/before 1-1-2000.Now d is something special. We cannot just take the days since 1-1-2000, because a solar year is longer than our year of 365 days, hence we have leap years, etc. That is why we use Julian Day Numbers.

The formula I used to convert to Julian days =

$$d = 367*Y - \frac{7*(Y + ((M+9)/12))}{4} + \frac{275*M}{9} + D - 730530$$

 i = 0.0
N = 0.0
w = 282.9404 + 4.70935E-5 * d
a = 1.000000 (AU)
e = 0.016709 - 1.151E-9 * d
M = 356.0470 + 0.9856002585 * d.

There are 3 parameters that are depended on d. The argument of the perihelion, the shape of the orbit, and the mean anomaly. The factors are all describing the shape of the orbit, and the speed of the earth along this orbit.
M is easy, as this show a single degree per Julian Day.
The eccentricity changes over time, this is a very small number, so the changes are small, but they are there. They are caused by the gravitational influence of Jupiter and Saturn.
The argument of the perihelion is changing too over time, at a factor of 0,00005^o per day. (-0.1668523^o in 10 years). This is called the Apsidal Precession. The orbit of the earth itself tilts

 

Calculating the path of the sun, the final calculations.

So now we can calculate the true distance and angle on a specific day:

We can calculate the eccentric anomaly (or true anomaly) from Kepler's equation

M = E - e sin(E)

E = eccentric anomaly = M + e sin(E) =  M + (180/pi) * e * sin(M) * (1 + e * cos(M))

x = cos(E) - e
y = sin(E) * sqrt(1 - e*e)

r = sqrt(x² + y²)
v = arctan2( y, x )

 

The moon is influenced by a lot of things: The graviational forces of the sun, the planets and of course the earth itself. But even by the tides on earth, or the shape of the earth and the position of the moon at a particular moment.

Add these terms to the Moon's longitude (degrees):

    -1.274 * sin(Mm - 2*D)          (the Evection)
    +0.658 * sin(2*D)               (the Variation)
    -0.186 * sin(Ms)                (the Yearly Equation)
    -0.059 * sin(2*Mm - 2*D)
    -0.057 * sin(Mm - 2*D + Ms)
    +0.053 * sin(Mm + 2*D)
    +0.046 * sin(2*D - Ms)
    +0.041 * sin(Mm - Ms)
    -0.035 * sin(D)                 (the Parallactic Equation)
    -0.031 * sin(Mm + Ms)
    -0.015 * sin(2*F - 2*D)
    +0.011 * sin(Mm - 4*D)

Add these terms to the Moon's latitude (degrees):

    -0.173 * sin(F - 2*D)
    -0.055 * sin(Mm - F - 2*D)
    -0.046 * sin(Mm + F - 2*D)
    +0.033 * sin(F + 2*D)
    +0.017 * sin(2*Mm + F)

Add these terms to the Moon's distance (Earth radii):

    -0.58 * cos(Mm - 2*D)
    -0.46 * cos(2*D)

All perturbation terms that are smaller than 0.01 degrees in longitude or latitude and smaller than 0.1 Earth radii in distance have been omitted here. A few of the largest perturbation terms even have their own names! The Evection (the largest perturbation) was discovered already by Ptolemy a few thousand years ago (the Evection was one of Ptolemy's epicycles). The Variation and the Yearly Equation were both discovered by Tycho Brahe in the 16'th century.

The computations can be simplified by omitting the smaller perturbation terms. The error introduced by this seldom exceeds the sum of the amplitudes of the 4-5 largest omitted terms. If one only computes the three largest perturbation terms in longitude and the largest term in latitude, the error in longitude will rarley exceed 0.25 degrees, and in latitude 0.15 degrees.

 

What we just did was calculating the expected position of the earth in the orbit around the sun. When we want to calculate the angle we see the sun at, from the surface of the earth, we have to deal with one extra item: The tilt of the earth itself. Moving through the ecliptic plane we are tilted at a 23.4^o angle: The obliquity, and the main reason we have seasons. With the formula oblecl = 23.4393 - 3.563E-7 * d  we can calculate this for a specific day.

This is pretty cool, is it not? And doing some more test the data is very exact. So we show here that with some relative simple formulae we get an accurate calculation of what we see when we look outside.