<p>It's a bit late but hopefully it will help others in the future.</p>
<p>The calculations can be found on <a href="http://www.stargazing.net/kepler/moon2.html" rel="nofollow noreferrer">this site</a>, article from Keith Burnett (keith@xylem.demon.co.uk)</p>
<p>The first problem is that you use J2000 as time reference but it's not what the site takes as time reference!</p>
<p>&quot;Get the days to Dec 31st 0h 2000 - note, this is NOT same as J2000&quot;</p>
<p>So you have to do :</p>
<pre><code>let d = daysSinceJ2000(for: date) + 1.5
</code></pre>
<p>The second problem is you are missing one degrees to radians conversion when you calculate E, it should be</p>
<pre><code>let E = M.radians + e * sin(M.radians) * (1.0 + e * cos(M.radians))
</code></pre>
<p>The third problem is the calculation of the declination, it should be arctan, not arctan2 :</p>
<pre><code>var dec = atan(zE / sqrt(xE * xE + yE * yE) )
</code></pre>
<p>On a side note, when you do :</p>
<pre><code>var N = 125.1228 - 0.0529538083 * d // have to be normalized
    if N &lt; 0 {
      N += 360
    }
</code></pre>
<p>N is still negative, to put a negative value within (0,360) I did</p>
<pre><code>N = range(N)

func range(x: Double) -&gt; Double{
    if x &gt;= 0 {
        return x
    }
    else {
        return range(x: x + 360)
    }
}
</code></pre>
<p>For positive values, you can just use</p>
<pre><code>M = M.truncatingRemainder(dividingBy: 360)
</code></pre>
<p>Although I don't think it's really crucial for the calculations, I tried without and it worked fine</p>
<p>Cheers, Lucas</p>

CC BY-SA 4.0
2020-12-30T16:18:02.170Z
2020-12-30T16:18:02.170Z
lucascavatoni
<p>It took some effort, but in the end, I successfully calculated the geographical point of the Sun. I assumed that I could do the same for the Moon.</p>

<p>My plan is:</p>

<ol>
<li>Find RA and declination of the Moon</li>
<li>Find GHA of the Moon, using the same math as for the Sun</li>
</ol>

<p>However, after following <a href="http://njsas.org/projects/tidal_forces/altaz/pausch/ppcomp.html#6" rel="nofollow noreferrer">this article</a> almost to the letter (except for precession and perturbations), my calculated answer is much different from <a href="https://theskylive.com/moon-info" rel="nofollow noreferrer">this site</a>. I've examined the article and my code many times, but can't see an error.</p>

<p>I use Swift programming language, but here it's merely a way to represent my calculations and results.</p>

<p><a href="https://i.sstatic.net/eNVEJ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/eNVEJ.png" alt="my code and results"></a></p>

<p>If anyone would like to try it for yourself, here's my code in text format:</p>

<pre><code>import Foundation

extension Double {
  var radians: Double {
    return self * .pi / 180
  }

  var degrees: Double {
    return self  * 180 / .pi
  }

  var rangedUpTo360: Double {
    let multiples = (self/360).rounded(.down)
    return self - 360*multiples
  }

  var rangedUpTo86400: Double {
    let multiples = (self/86400).rounded(.down)
    return self - 86400*multiples
  }
}

extension Date {
  var greenwich: Date {
    let greenwichDate = self.addingTimeInterval(TimeInterval(-TimeZone.current.secondsFromGMT()))
    return greenwichDate
  }
}

class Astronomy {

  static func julianDays(since: Date) -&gt; Double {
    return since.timeIntervalSince1970 / 86400 + 2_440_587.5
  }

  static func daysSinceJ2000(for date: Date) -&gt; Double {
    return julianDays(since: date) - 2_451_545
  }

  static func sunAxialTilt(for daysSinceJ200: Double) -&gt; Double {
    return 23.4393 - 3.563e-7 * daysSinceJ200
  }

  static func greenwichHourAngle(rightAscension: Double, date: Date) -&gt; Double {
    let greenwichTimeZone = TimeZone(secondsFromGMT: 0)!
    var calendar = Calendar.current
    calendar.timeZone = greenwichTimeZone
    let components = calendar.dateComponents([.hour, .minute, .second], from: date)
    let secondsOfDay = Double(components.hour!)*3600 + Double(components.minute!)*60 + Double(components.second!)

    let centuriesSinceJ200 = daysSinceJ2000(for: date.greenwich) / 36525
    let greenwichSiderealTime = 24110.54841 + 8640184.812866*centuriesSinceJ200 + 0.093104*pow(centuriesSinceJ200, 2) - 0.0000062*pow(centuriesSinceJ200, 3)
    let earthSiderealRotationRate = 1.00273790935 + 5.9e-11*centuriesSinceJ200
    let earthRotation = earthSiderealRotationRate*secondsOfDay
    let gmstSeconds = greenwichSiderealTime + earthRotation
    let gmstSecondsNormalized = gmstSeconds.rangedUpTo86400
    let gmst = gmstSecondsNormalized / 3600 * 15
    return gmst.rangedUpTo360 - rightAscension
  }

  static func moonPosition(on date: Date) -&gt; (declination: Double, RA: Double, distance: Double) {
    date
    // day number
    let d = daysSinceJ2000(for: date)
    // longitude of the ascending node
    var N = 125.1228 - 0.0529538083 * d // have to be normalized
    if N &lt; 0 {
      N += 360
    }
    // inclination to the ecliptic (plane of the Earth's orbit)
    let i = 5.1454
    // argument of perihelion
    var w = 318.0634 + 0.1643573223 * d // have to be normalized
    w = w.rangedUpTo360
    // semi-major axis, or mean distance from Sun
    let a = 60.2666 // (Earth radii)
    // eccentricity (0=circle, 0-1=ellipse, 1=parabola)
    let e = 0.054900
    // mean anomaly (0 at perihelion; increases uniformly with time)
    var M = 115.3654 + 13.0649929509 * d // have to be normalized
    M = M.rangedUpTo360
    // eccentric anomaly
    let E = M + e * sin(M.radians) * (1.0 + e * cos(M.radians))
    let xv = a * (cos(E.radians) - e)
    let yv = a * (sqrt(1.0 - e*e) * sin(E.radians))
    // true anomaly (angle between position and perihelion)
    let v = atan2(yv, xv).degrees
    let r = sqrt(pow(xv, 2) + pow(yv, 2))
    // geocentric ecliptical position
    let xG = r * (cos(N.radians) * cos((v + w).radians) - sin(N.radians) * sin((v + w).radians) * cos(i.radians))
    let yG = r * (sin(N.radians) * cos((v + w).radians) + cos(N.radians) * sin((v + w).radians) * cos(i.radians))
    let zG = r * (sin((v + w).radians) * sin(i.radians))
    // equatorial coordinates
    let ecl = sunAxialTilt(for: d)
    let xE = xG
    let yE = yG * cos(ecl.radians) - zG * sin(ecl.radians)
    let zE = yG * sin(ecl.radians) + zG * cos(ecl.radians)
    //right scension
    var rA = atan2(yE, xE)
    rA = rA.degrees
    // declinatino
    var dec = atan2(zE, sqrt(xE * xE + yE * yE) )
    dec = dec.degrees
    // geocentric distance
    let rG = sqrt(xE * xE + yE * yE + zE * zE)
    return (declination: dec, RA: rA, distance: rG)
  }
}

Astronomy.moonPosition(on: Date().greenwich)
</code></pre>

CC BY-SA 4.0
2019-02-09T18:03:42.693Z
How big are the errors that you are seeing ? I compared theskylive.com link you gave and the apparent Moon position it gave for Greenwich at 2019-02-09 16:52 UTC was off by (2.8,-37.2) arcminutes in RA,Dec respectively compared to JPL HORIZONS and SLALIB. As noted in the njsas.org article, you need the perturbations included to get the position better than 2 degrees - the Moon is so close and moves relatively quickly that small errors in time etc get magnified. [Lunar motion theory](https://en.wikipedia.org/wiki/Lunar_theory) is surprisingly complicated
astrosnapper
CC BY-SA 4.0
2019-02-09T20:57:09.103Z
@astrosnapper My code returns these values for Feb 9, 2019 at 8:16 PM:
- RA: -0.5199
- Dec: -5.1769
While the website returns:
- RA: 15.56
- Dec: 1.073
If it were to the lack of perturbation calculations, shouldn't the difference be in 2 degrees range?
Somewhere I read that if atan2 returns a negative number in the calculations, I need to add .pi to it. Is it so? And does it mean I have to add 3.14 or 180? I tried though, but it results in a much bigger error.
bitemybyte
CC BY-SA 4.0
2019-02-10T19:28:07.820Z
Here's what I get using njsas's formulae for your datetime (assuming it's UTC , for the geocenter and without the perturbations): `N=115.460502, i=5.1454, w=25.416300,a=60.2666,e=0.0549, M=240.048858` (`N,i,w,M` are in degrees, `a` in Earth radii, `e` unitless). `xg,yg,zg` are `(60.938805003268435, 9.606098833264198, -5.326251364267174)` and obliquity of the ecliptic, `ecl=23.436812725124444` degrees. This gives `RA=0.17750561458751776`, `Dec= -0.017218433367563225` (both radians)
or 00h40m40.88s -00d 59' 11.56". HORIZONS gives: 01 03 06.19 +01 09 42.9 so ~2deg out in Dec
astrosnapper
CC BY-SA 4.0
2019-02-10T19:36:47.453Z
One thing to be careful is whether your quantities are in the right units. In the njsas formulae many things are done in degrees initially but there seems to be some missing conversions to radians when used in trig functions. For the eccentric anomaly, `E`, both the `sin(M)` and `cos(M)` need conversions. Similarly the use of `E` in the `xv,yv` formulae (Sec. 6) need converting as does the use of `N` and `i` in the formula for 3D position (Sec. 7) and `ecl` in the conversion from ecliptic to equatorial (Sec. 12).
astrosnapper
CC BY-SA 4.0
2019-02-11T18:27:47.097Z
Hm, for the same date (Feb 9 20:16) I get these results:
d=6979.2194;
N=115.5465;
i=same;
w=25.1492;
a=same;
e=same;
M=218.8182;
E=218.7853;
xg,yg zg=(62.5865, -2.7821, -4.9766);
ecl=23.4368;
RA=-0.0091580905076487;
declination=-0.09038631707987221;
Which is different on more then 2 degrees from yours. Would you mind sharing what was your d, E, xV, yV and v? I believe I've handled the conversions properly, as all the sin/cos inputs in my code get converted to radians. ".degrees" and ".radians" in my code convert the value to from rad to degrees and from degrees to rad respectfully
bitemybyte
CC BY-SA 4.0
2019-02-11T19:02:06.987Z
Your value for `d` isn't right. I don't know how Swift works but are you sure it's doing integer division as required? In Python, for `y=2019 ; m=2; D=9; UT=20.0+(16/60.0)`, `d = 367*y - 7 * ( y + (m+9)//12 ) // 4 + 275*m//9 + D - 730530` gives '6980' and adding the fractional UT part produces `d=6980.844444444445`. This agrees with what you get if you use [USNO JD calc](https://aa.usno.navy.mil/data/docs/JulianDate.php) and calc. for 2019-02-09T20:16 and J2000 (1999-12-31T00:00:00) and subtract the two
astrosnapper
CC BY-SA 4.0
2019-02-16T10:27:31.507Z
@astrosnapper I corrected my formula for `d` so that for 9.02.19 20:16 it produces exactly the same value as your calculations, thank you. I think you might have a typo in calculating UT part, as for it to be 0.844444444445 I had to: `UT = (hours*3600 + minutes*60.0 + seconds)/86400`. My `xg, yg, zg` are still off `(58.60271962971655, 19.22395067561122, -5.508631049423485)`. Considering my N, i, w, M are exactly the same as yours, the only possible place for an error would be true anomaly `v`, `-122.6862201425946` in my case. (`237.3137798574054` if I add 360)
bitemybyte
CC BY-SA 4.0
2019-02-16T10:31:39.253Z
I get `v` by `atan2(yv, xv)`, where `yv = a * (sqrt(1.0 - e*e) * sin(E)) = -52.11505497081653`, `xv = a * (cos(E) - e) = -33.43957369651925` and `E = M + e * sin(M) * (1.0 + e * cos(M)) = 240.0025936288414`
bitemybyte
<the-moon><declination><right-ascension>
Moon equatorial coordinates
bitemybyte
# Moon equatorial coordinates

(Asked by: username_0 on 2019-02-09)

It took some effort, but in the end, I successfully calculated the geographical point of the Sun. I assumed that I could do the same for the Moon.
My plan is:

Find RA and declination of the Moon
Find GHA of the Moon, using the same math as for the Sun

However, after following this article almost to the letter (except for precession and perturbations), my calculated answer is much different from this site. I've examined the article and my code many times, but can't see an error.
I use Swift programming language, but here it's merely a way to represent my calculations and results.

If anyone would like to try it for yourself, here's my code in text format:
import Foundation

extension Double {
  var radians: Double {
    return self * .pi / 180
  }

  var degrees: Double {
    return self  * 180 / .pi
  }

  var rangedUpTo360: Double {
    let multiples = (self/360).rounded(.down)
    return self - 360*multiples
  }

  var rangedUpTo86400: Double {
    let multiples = (self/86400).rounded(.down)
    return self - 86400*multiples
  }
}

extension Date {
  var greenwich: Date {
    let greenwichDate = self.addingTimeInterval(TimeInterval(-TimeZone.current.secondsFromGMT()))
    return greenwichDate
  }
}

class Astronomy {

  static func julianDays(since: Date) -> Double {
    return since.timeIntervalSince1970 / 86400 + 2_440_587.5
  }

  static func daysSinceJ2000(for date: Date) -> Double {
    return julianDays(since: date) - 2_451_545
  }

  static func sunAxialTilt(for daysSinceJ200: Double) -> Double {
    return 23.4393 - 3.563e-7 * daysSinceJ200
  }

  static func greenwichHourAngle(rightAscension: Double, date: Date) -> Double {
    let greenwichTimeZone = TimeZone(secondsFromGMT: 0)!
    var calendar = Calendar.current
    calendar.timeZone = greenwichTimeZone
    let components = calendar.dateComponents([.hour, .minute, .second], from: date)
    let secondsOfDay = Double(components.hour!)*3600 + Double(components.minute!)*60 + Double(components.second!)

    let centuriesSinceJ200 = daysSinceJ2000(for: date.greenwich) / 36525
    let greenwichSiderealTime = 24110.54841 + 8640184.812866*centuriesSinceJ200 + 0.093104*pow(centuriesSinceJ200, 2) - 0.0000062*pow(centuriesSinceJ200, 3)
    let earthSiderealRotationRate = 1.00273790935 + 5.9e-11*centuriesSinceJ200
    let earthRotation = earthSiderealRotationRate*secondsOfDay
    let gmstSeconds = greenwichSiderealTime + earthRotation
    let gmstSecondsNormalized = gmstSeconds.rangedUpTo86400
    let gmst = gmstSecondsNormalized / 3600 * 15
    return gmst.rangedUpTo360 - rightAscension
  }

  static func moonPosition(on date: Date) -> (declination: Double, RA: Double, distance: Double) {
    date
    // day number
    let d = daysSinceJ2000(for: date)
    // longitude of the ascending node
    var N = 125.1228 - 0.0529538083 * d // have to be normalized
    if N < 0 {
      N += 360
    }
    // inclination to the ecliptic (plane of the Earth's orbit)
    let i = 5.1454
    // argument of perihelion
    var w = 318.0634 + 0.1643573223 * d // have to be normalized
    w = w.rangedUpTo360
    // semi-major axis, or mean distance from Sun
    let a = 60.2666 // (Earth radii)
    // eccentricity (0=circle, 0-1=ellipse, 1=parabola)
    let e = 0.054900
    // mean anomaly (0 at perihelion; increases uniformly with time)
    var M = 115.3654 + 13.0649929509 * d // have to be normalized
    M = M.rangedUpTo360
    // eccentric anomaly
    let E = M + e * sin(M.radians) * (1.0 + e * cos(M.radians))
    let xv = a * (cos(E.radians) - e)
    let yv = a * (sqrt(1.0 - e*e) * sin(E.radians))
    // true anomaly (angle between position and perihelion)
    let v = atan2(yv, xv).degrees
    let r = sqrt(pow(xv, 2) + pow(yv, 2))
    // geocentric ecliptical position
    let xG = r * (cos(N.radians) * cos((v + w).radians) - sin(N.radians) * sin((v + w).radians) * cos(i.radians))
    let yG = r * (sin(N.radians) * cos((v + w).radians) + cos(N.radians) * sin((v + w).radians) * cos(i.radians))
    let zG = r * (sin((v + w).radians) * sin(i.radians))
    // equatorial coordinates
    let ecl = sunAxialTilt(for: d)
    let xE = xG
    let yE = yG * cos(ecl.radians) - zG * sin(ecl.radians)
    let zE = yG * sin(ecl.radians) + zG * cos(ecl.radians)
    //right scension
    var rA = atan2(yE, xE)
    rA = rA.degrees
    // declinatino
    var dec = atan2(zE, sqrt(xE * xE + yE * yE) )
    dec = dec.degrees
    // geocentric distance
    let rG = sqrt(xE * xE + yE * yE + zE * zE)
    return (declination: dec, RA: rA, distance: rG)
  }
}

Astronomy.moonPosition(on: Date().greenwich)

**Tags:** the-moon, declination, right-ascension
> **Comments:**
> **username_1**: How big are the errors that you are seeing ? I compared theskylive.com link you gave and the apparent Moon position it gave for Greenwich at 2019-02-09 16:52 UTC was off by (2.8,-37.2) arcminutes in RA,Dec respectively compared to JPL HORIZONS and SLALIB. As noted in the njsas.org article, you need the perturbations included to get the position better than 2 degrees - the Moon is so close and moves relatively quickly that small errors in time etc get magnified. [Lunar motion theory](https://en.wikipedia.org/wiki/Lunar_theory) is surprisingly complicated  
> **username_0**: @astrosnapper My code returns these values for Feb 9, 2019 at 8:16 PM: - RA: -0.5199 - Dec: -5.1769 While the website returns: - RA: 15.56 - Dec: 1.073 If it were to the lack of perturbation calculations, shouldn't the difference be in 2 degrees range? Somewhere I read that if atan2 returns a negative number in the calculations, I need to add .pi to it. Is it so? And does it mean I have to add 3.14 or 180? I tried though, but it results in a much bigger error.  
> **username_1**: Here's what I get using njsas's formulae for your datetime (assuming it's UTC , for the geocenter and without the perturbations): `N=115.460502, i=5.1454, w=25.416300,a=60.2666,e=0.0549, M=240.048858` (`N,i,w,M` are in degrees, `a` in Earth radii, `e` unitless). `xg,yg,zg` are `(60.938805003268435, 9.606098833264198, -5.326251364267174)` and obliquity of the ecliptic, `ecl=23.436812725124444` degrees. This gives `RA=0.17750561458751776`, `Dec= -0.017218433367563225` (both radians) or 00h40m40.88s -00d 59' 11.56". HORIZONS gives: 01 03 06.19 +01 09 42.9 so ~2deg out in Dec  
> **username_1**: One thing to be careful is whether your quantities are in the right units. In the njsas formulae many things are done in degrees initially but there seems to be some missing conversions to radians when used in trig functions. For the eccentric anomaly, `E`, both the `sin(M)` and `cos(M)` need conversions. Similarly the use of `E` in the `xv,yv` formulae (Sec. 6) need converting as does the use of `N` and `i` in the formula for 3D position (Sec. 7) and `ecl` in the conversion from ecliptic to equatorial (Sec. 12).  
> **username_0**: Hm, for the same date (Feb 9 20:16) I get these results: d=6979.2194; N=115.5465; i=same; w=25.1492; a=same; e=same; M=218.8182; E=218.7853; xg,yg zg=(62.5865, -2.7821, -4.9766); ecl=23.4368; RA=-0.0091580905076487; declination=-0.09038631707987221; Which is different on more then 2 degrees from yours. Would you mind sharing what was your d, E, xV, yV and v? I believe I've handled the conversions properly, as all the sin/cos inputs in my code get converted to radians. ".degrees" and ".radians" in my code convert the value to from rad to degrees and from degrees to rad respectfully  
> **username_1**: Your value for `d` isn't right. I don't know how Swift works but are you sure it's doing integer division as required? In Python, for `y=2019 ; m=2; D=9; UT=20.0+(16/60.0)`, `d = 367*y - 7 * ( y + (m+9)//12 ) // 4 + 275*m//9 + D - 730530` gives '6980' and adding the fractional UT part produces `d=6980.844444444445`. This agrees with what you get if you use [USNO JD calc](https://aa.usno.navy.mil/data/docs/JulianDate.php) and calc. for 2019-02-09T20:16 and J2000 (1999-12-31T00:00:00) and subtract the two  
> **username_0**: @astrosnapper I corrected my formula for `d` so that for 9.02.19 20:16 it produces exactly the same value as your calculations, thank you. I think you might have a typo in calculating UT part, as for it to be 0.844444444445 I had to: `UT = (hours*3600 + minutes*60.0 + seconds)/86400`. My `xg, yg, zg` are still off `(58.60271962971655, 19.22395067561122, -5.508631049423485)`. Considering my N, i, w, M are exactly the same as yours, the only possible place for an error would be true anomaly `v`, `-122.6862201425946` in my case. (`237.3137798574054` if I add 360)  
> **username_0**: I get `v` by `atan2(yv, xv)`, where `yv = a * (sqrt(1.0 - e*e) * sin(E)) = -52.11505497081653`, `xv = a * (cos(E) - e) = -33.43957369651925` and `E = M + e * sin(M) * (1.0 + e * cos(M)) = 240.0025936288414`  

---

## 1 answers

### Answer by username_2 on 2020-12-30. Score: 3

It's a bit late but hopefully it will help others in the future.
The calculations can be found on this site, article from Keith Burnett (keith@xylem.demon.co.uk)
The first problem is that you use J2000 as time reference but it's not what the site takes as time reference!
"Get the days to Dec 31st 0h 2000 - note, this is NOT same as J2000"
So you have to do :
let d = daysSinceJ2000(for: date) + 1.5

The second problem is you are missing one degrees to radians conversion when you calculate E, it should be
let E = M.radians + e * sin(M.radians) * (1.0 + e * cos(M.radians))

The third problem is the calculation of the declination, it should be arctan, not arctan2 :
var dec = atan(zE / sqrt(xE * xE + yE * yE) )

On a side note, when you do :
var N = 125.1228 - 0.0529538083 * d // have to be normalized
    if N < 0 {
      N += 360
    }

N is still negative, to put a negative value within (0,360) I did
N = range(N)

func range(x: Double) -> Double{
    if x >= 0 {
        return x
    }
    else {
        return range(x: x + 360)
    }
}

For positive values, you can just use
M = M.truncatingRemainder(dividingBy: 360)

Although I don't think it's really crucial for the calculations, I tried without and it worked fine
Cheers, Lucas


---

