How to calculate the position of the sun in long/lat?

Question

I have a world map, and I want to know, given a timestamp, where is the sun located right now, in terms of long/lat.

Here is a website that does that: https://www.timeanddate.com/worldclock/sunearth.html

But they do not share the calculations.

How can I calculate the sun's position in earth gps terms given a timestamp?

Edit: some context

I have a plane that is determined by 2 points on earth, and that plain is perpendicular to the earth in those points. I want to know where is the sun relative to that plane, so I can know if the sun is to the right of it, or to the left of it.

Instead of working in 3 dimensions, I decided to work in 2. I already have that plane in 2D (just a line on the map), now I need to know where the sun is

Edit2: Using the wiki page here, this is what I got so far

Note: The interesting part is getCurrentPosition

class Sun {
    private julianTime(timestamp: number) {
        // Source: http://stackoverflow.com/questions/466321/convert-unix-timestamp-to-julian
        return timestamp / 86400 + 2440587.5;
    }

    private convertEclipticEquatorial(lambda: number, beta: number, R: number, n) {
        const epsilon = 23.439 - 0.0000004 * n; // Obliquity of the ecliptic
        console.log("epsilon", epsilon);
        const alpha = Math.atan2(Degrees.cos(epsilon) * Degrees.sin(lambda), Degrees.cos(lambda));
        console.log("alpha", alpha);
        const delta = Math.asin(Degrees.sin(epsilon) * Degrees.sin(lambda)); // declination
        console.log("delta", delta);
    }

    getCurrentPosition(timestamp: number) {
        const JD = this.julianTime(timestamp);
        console.log("JD", JD);

        // Source: https://en.wikipedia.org/wiki/Position_of_the_Sun
        const n = JD - 2451545;
        console.log("n", n);
        const L = (280.460 + 0.9856474 * n) % 360;
        console.log("L", L);
        const g = (357.528 + 0.9856003 * n) % 360;
        console.log("g", g);
        const lambda = (L + 1.915 * Degrees.sin(g) + 0.020 * Degrees.sin(2 * g)) % 360;
        console.log("lambda", lambda);
        const R = 1.00014 - 0.01671 * Degrees.cos(g) - 0.00014 * Degrees.cos(2 * g); // Distance
        console.log("R", R);

        this.convertEclipticEquatorial(lambda, 0, R, n);
    }
}

class Degrees {
    static sin(angle: number) {
        return Math.sin(angle / 180 * Math.PI);
    }
    static cos(angle: number) {
        return Math.cos(angle / 180 * Math.PI);
    }
}

console.log(new Sun().getCurrentPosition(new Date().getTime() / 1000));

Running this now, gives:

JD 2457844.130512928 (can confirm, this is correct)

n 6299.1305129281245 (simple math)

L 9.181612328272422 (simple math)

g 85.95292328111373 (simple math)

lambda 11.09465308568091

R 0.9990992788207762

epsilon 23.436480347794827 (can confirm, that is roughly correct)

alpha 0.17801325468993906

delta 0.0766106706512887

Answer 1

By latitude and longitude of the sun, I think you mean latitude and longitude of a point on the earth's surface where the rays of the sun fall perpendicular to the surface.

In astronavigation, this point is called the Geographical Position (GP) of the sun and terminology to describe this point is slightly different:-

Latitude of the GP is called Declination, and it is numerically equal to what one might call "Latitude". It is measured in degrees North or South the equator.

Longitude of the GP is called Greenwich Hour Angle (GHA). Just like longitude, GHA is measured East or West of Greenwich Meridian.

Describing the sun's position in this way is called the Equatorial Coordinate System.

The Wikipedia page you have used gives the additional formulae below to convert from Ecliptic Coordinates (the steps you followed) to Equatorial Coordinates.

Follow these steps to determine Declination $\delta$ (Equivalent to Latitude).

To get the sun's GHA, you need to first calculate $RA$ (Right Ascension) as shown above and then convert this to GHA by these= rules...

$SHA° = 360 – 15 · RA[h] $

$GHA_{Aries}° = 15 · GST[h] $

$GHA = SHA + GHA_{Aries}$

For the second step, you will need $GST$ (Greenwich Sidereal Time) which is you can calculate if you know UT (UT can be obtained from the time shown by a clock). Use this link to convert directly from UT to GST

If you want a simple, concise way to be able to calculate the positions of all manner of celestial bodies, the sun, moon, planets, stars, then I highly recommended the book Astronomical Formulae for Calculators by Jean Meeus. It is the simplest, most comprehensive collection of astronomical formulae you are likely find anywhere on earth.

Answer 2

If you don't need this precisely, the easy way to figure this out would be the following:

You can assume that at noon UTC on one of the equinoxes, the Sun will be directly over (0N, 0E). The longitude of the sub-Solar point will then progress evenly over the course of a day: 4 hours later it's at 60W, 12 hours later at 180W, 18 hours later 90E, etc.

The latitude of the sub-Solar point will either increase or decrease (depending upon the starting equinox) as a sine curve with a period of year and an amplitude equal to the axial tilt of the Earth (23.43705 degrees).

That will be good to ~5 degrees over the course of the year.

To get a more precise answer, you'd need to calculate the Analemma over the course of the year, and position the sub-Solar point accordingly. Formally, this would mean calculating the Equation of Time and transforming the resulting time offset into a position offset for the sub-Solar point. E.g., what's happening in this image.

I'm not familiar with how to do that, so hopefully someone else can help - or the simple way above works.

Answer 3

I've implemented what Thomas suggested as an approximation. It seems to do its job for a small rendered globe toy in JavaScript.

function getSunEuler(date) {
        const now = date || new Date();

        // The boilerplate: fiddling with dates
        const soy = (new Date(now.getFullYear(), 0, 0)).getTime();
        const eoy = (new Date(now.getFullYear() + 1, 0, 0)).getTime();
        const nows = now.getTime();
        const poy = (nows - soy) / (eoy - soy);

        const secs = now.getUTCMilliseconds() / 1e3
                + now.getUTCSeconds()
                + 60 * (now.getUTCMinutes() + 60 * now.getUTCHours());
        const pod = secs / 86400; // leap secs? nah.

        // The actual magic
        const lat = (-pod + 0.5) * Math.PI * 2;
        const lon = Math.sin((poy - .22) * Math.PI * 2) * .41;

        return new THREE.Euler(0, lat, lon, 'YZX');
}

The returned Euler angle can be used to transform an X unit vector to find a point on a sphere where 0 N 0 E is at (1, 0, 0).

I got curious how accurate the solution is and fiddled with Skyfield a bit to get closer to the truth:

from skyfield.api import Topos, load
ts = load.timescale()
planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
t = ts.utc(2019, 6, 21, 12, 0, 0)

# A rather correct answer can probably be obtained from
earth.at(t).observe(sun).apparent().position.au
# but I ended up being too lazy to calculate the angles myself

# This serves as a convenient workaround for manual verification
np = earth + Topos('90 N', '180 W')
np.at(t).observe(sun).apparent().altaz()

The approximation seems to be indeed accurate within a few degrees, but is in an different coordinate system. I'd love it if someone could provide a clean solution using skyfield.