1

I'm trying to implement the Astronomical Almanac's approximation of the sun's position.

My code looks like this:

var now := OS.get_datetime()
var lat := deg2rad(player_latlon.x)
var lon := deg2rad(player_latlon.y)

### Astronomical Almanac, https://en.wikipedia.org/wiki/Position_of_the_Sun
# days since J2000.0
var n := Utils.get_julian_from_unix(OS.get_unix_time()) - 2451545.0
# mean longitude of sun
var L := fmod( deg2rad(280.46) + deg2rad(0.9856474)*n, 2*PI)
# mean anomaly of sun
var g := fmod( deg2rad(357.528) + deg2rad(0.9856003)*n, 2*PI)

### conversions, same page
# ecliptic longitude
var eclon := L + deg2rad(1.915)*sin(g) + deg2rad(0.020)*sin(2*g)
# ecliptic latitude (always < 0.00033deg)
var eclat := 0

### conversions, https://en.wikipedia.org/wiki/Celestial_coordinate_system
# obliquity of the ecliptic
var obl := deg2rad(23.4)
# equatorial right ascension
var eqra := atan( ( sin(eclon)*cos(obl) - tan(eclat)*sin(obl) )/cos(eclon) )
# equatorial declination
var eqde := asin( sin(eclat)*cos(obl) + cos(eclat)*sin(eclon)*sin(obl) )
# hour angle (assumes `now` in local sidereal time)
var eqha := deg2rad((now.hour as int) + (now.minute as int)) - eqra
# horizontal azimuth
var haz := atan( sin(eqha) / ( cos(eqha)*sin(lat) - tan(eqde)*cos(lat) ) )
# horizontal altitude
var hal := asin( sin(lat)*sin(eqde) + cos(lat)*cos(eqde)*cos(eqha) )

var sun_latlon := Vector2(rad2deg(hal),rad2deg(haz))

Reality says that close to where i'm located, the results should look like this:enter image description here

(azim,elev) == (292.7,6.09) in that picture are what i expect my code to return (at least approximately) as sun_latlon

but with the exact same data as the input, i get this:

sun_latlon: (3.291654, -63.823071)

which clearly can't be right, but i quadruple-checked every single line for days now, and it all looks correct to me.

nonchip
  • 1,084
  • 1
  • 18
  • 36
  • 1
    Without knowing the language you are writing in, I must say that my suspicions fall upon the line where you calculate `eqha`. Is `now` really in sidereal time? From the first line of your posted code, it looked to me like it would be in local civil time. And when you add `now.hour` and `now.minute` together, can you really treat them as being in units of degrees? One final thing - you use `atan()` to calculate `haz`. `haz` can vary from -pi to + pi, but `atan()` will only return a result in the range -pi/2 to +pi/2. Try `atan2()`. Maybe `eqra` is also affected by this. – DavidHoadley Sep 22 '20 at 02:09
  • @DavidHoadley sounds like you're right. in the meantime i've solved the issue in https://gitlab.com/-/snippets/1999214 but don't ask me how exactly, that's the result of lots of trial and error translating an algo i still don't quite understand to gdscript... the sidereal time issue was solved using *some* equation i found for "greenwich mean sidereal time for 0:00 at a given julian date" + the local hour offset, and yeah lets just say it took a WHILE to figure out all those angle conversions... – nonchip Oct 12 '20 at 16:21

0 Answers0