5

I'm trying to develop an Excel macro to compute when a target will rise above a specified altitude given the date and location. I can get the correct answer with Stellarium and using a Python program I wrote that leverages the PyEphem library, but for some reason I cannot get the right answer when I try to implement the calculations myself using Excel. I'm hoping someone can point out the error of my ways.

These are the input variables:

Date: 10/18/2023

Altitude: 30°

Target: M74

RA:   01:36:43.1
Dec:  +15° 47' 08"

Location:

Lat:  32° 36' 50"
Long: -116° 20' 01"

According to Stellarium and the Python program using PyEphem, the correct answer is 20:16:55 PDT (03:16:55 UTC the following day). But I cannot arrive at the same answer and don’t know where I’m going wrong.

As I understand it, the steps are as follows:

  1. Compute the LHA

    Cos(LHA) = [SIN(Alt) - SIN(Lat) * SIN(Dec)] / (COS(Lat) * COS(Dec))

    Cos(LHA) =[0.5 - 0.538974977 * 0.272037656] / (0.842321776 * 0.962286607)

    Cos(LHA) = 0.435971065

    LHA = Acos(0.435971065) = 1.119679324 rad = 64.2° = 4.276860089 hours = 04:16:37

    This result is close to the LHA for M74 in the western sky according to Stellarium (04:16:54). Subtracting that value from 2*Pi to get the LHA when M74 is in the eastern sky, we get:

    5.163505973 rad = 295.8° = 19.72 hours = 19:43:23 (Stellarium has 19:43:05)

    I’m not sure why there is a difference of 17-18 seconds, but it’s close enough that I’m comfortable with the answer.

  2. Compute the LST

    LST(E) = LHA(E) + RA = 5.163505973 + 0.422013334 = 5.585519307 rad = 320.0° = 21.34 hours = 21:20:06

    Stellarium has 21:20:58.5 for the Apparent Sidereal Time so the difference is now 52 seconds, but again, it seems close enough to be on the right track.

  3. Compute the GMST

    GMST(E) = LST(E) – Long = 5.585519307 – (-2.03040451) = 7.61592393525783

    That value is greater than 2Pi; do I need to subtract 2Pi or leave it as is?

  4. Then, using the Meeus formulas, we have:

    JD0 = 2460235.5 (for 10/18/2023)

    T = (JD0 - 2451545) / 36524.22 = 0.237938003877975

    UT = (GMST(E-hours) - 6.697374558 - 2400.051336 * T - 0.000025862 * (T * T)) / 1.0027379093

    UT = -547.172006240208  Mod 24 = 4.82799375979209 = 04:49:41

    LT (PDT) = UT – 7 = 21:49:41

    That answer is incorrect by over 1-1/2 hours; Stellarium and PyEphem both show the correct time should be 20:16:54 PDT

It appears that the calculations are correct through Step 3, but it goes awry sometime during Step 4. Can anyone see where I’m getting off track?

BillR
  • 61
  • 3
  • 1
    For Step 1: difference might be precession/nutation - are you coordinates for M74 ICRS/J2000 ones or not ? In Step 3: Yes, you need to normalize this into the range 0...2PI – astrosnapper Oct 20 '23 at 20:24
  • 1
    The M74 coordinates are J2000. If I normalize the GMST values into the range of 0 to 2*PI in the results only change by ~4 minutes, so I'm doing something wrong elsewhere... – BillR Oct 20 '23 at 21:15
  • 1
    Since you have the Meeus book, use the Rise, Set, Transit algorithm, but set $ h_0 $ to the alt you want. – Greg Miller Oct 21 '23 at 04:34
  • 1
    @uhoh I did it. – Snack Exchange Oct 21 '23 at 10:27
  • 1
    @GregMiller - I don't have the book. I've just done a lot of Internet searches on the subject and thought I had the right equations based on what I've read. I was hoping someone could point out what I have wrong. – BillR Oct 21 '23 at 16:48
  • 1
    The equations are here: https://celestialprogramming.com/risesetalgorithm.html – Greg Miller Oct 21 '23 at 18:52
  • @GregMiller - Thanks for that reference. I decided to try something simple and just compute the meridian transit using those equations but am still not getting the right answer. One problem I encountered is that the web site appears to assume west longitude as a positive value, is that correct? Because if I enter a negative value for my west longitude, the result isn't anywhere near close to being correct. But if I enter it as a positive value, the calculator yield precisely the correct answer. But when I use the equations to compute the transit time, the result is ~30 minutes off. – BillR Oct 22 '23 at 16:44
  • Using the values in the OP: JD0 = 2460235 T = (JD0 - 2451545) / 36525 = 0.237932922655715 Theta0 = 3137426.22907707 % 360 = 26.2290770742111 Transit(deg) = 359.648 = 23.976 hours = ~23:58:34 But the correct transit time is 00:37:08, about 40 minutes later. I can't figure out what I'm doing wrong. – BillR Oct 22 '23 at 16:53
  • 1
    @GregMiller - I think I found the problem. There appears to be an error with one of the formulas on that website you referenced. It shows transit = (Dec + Long - GMST)/360, but in your C source code, you have transit - (RA + Long - GMST)/360. When I use RA instead of Dec in that equation, I get the right answer. – BillR Oct 26 '23 at 15:24

1 Answers1

1

Thanks to @GregMiller for pointing me to his website that has the formulas I needed.

enter image description here

BillR
  • 61
  • 3
  • 4
    Congratulations, problem solved! But to make this a proper Stack Exchange answer, please add the formula here or at least an explanation of it. "The answer is in this link" is called a link-only answer and these are not considered good answers because they're not answers. They require each reader to go offsite to actually find the answer, and if/when the link breaks they become totally useless. Thanks! – uhoh Oct 26 '23 at 21:37