1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ What tools do you have access to? How are you calculating this? If you're trying to do this from scratch, you're in for a rough time. $\endgroup$ – zephyr Mar 29 '17 at 15:07
  • $\begingroup$ @zephyr Well, I have internet access, not much more. I love cracking tough problems. Perhaps can you point me in the right direction? $\endgroup$ – Amit Mar 29 '17 at 15:16
  • $\begingroup$ This is a really tough problem though. There are people at NASA whos entire job is to do just stuff like this. There's no single source I could point you to that would help you do this yourself. Basically though, you'll need an ephemeris for a particular time period, the ability to progress that ephemeris forward in time to your specific date, and then lots of math to convert the ephemeris info into what you need. Your best bet might be to look at some code by someone who's already done something like this (e.g., clearskyinstitute.com/xephem) and copy what they've done. $\endgroup$ – zephyr Mar 29 '17 at 15:21
  • 2
    $\begingroup$ Glad I could help. I didn't know there was a wiki page with the relevant equations so that's good to know. Feel free to answer your own question and accept it if that information is sufficient for you. $\endgroup$ – zephyr Mar 29 '17 at 15:35
  • 1
    $\begingroup$ @zephyr I'll answer it once I have some code taking a timestamp and returning the sun's position $\endgroup$ – Amit Mar 29 '17 at 15:36
1
$\begingroup$

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.

enter image description here

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks for the answer! That indeed is a way to calculate this, and I will use it in case I fail the current route (as this is more or less a rough approximation). I have edited the question, with some code, and the result, however, I'm failing to understand what do I need - ecliptic/equatorial/horizontal coordinates. Given the numbers in the result, what do I need to do to get long/lat? $\endgroup$ – Amit Mar 31 '17 at 15:17
  • $\begingroup$ Amit, once you get the equatorial coordinates, you should be there. Equatorial coords are for objects in the sky, but you can "project" them down onto the surface of the Earth to find the sub-Solar point. That is: alpha = long., delta = lat. $\endgroup$ – Thomas Mar 31 '17 at 15:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.