Calculate Earth Zenith Point of Another Planet with Azimuth and Elevation of another Observer Point

Question

I need to find the Longitude of the closest point on Earth (zenith) to Jupiter with the Azimuth and Elevation data from another observation location at a given point in time.

Below is the source of my data points for one specific second in time:

Target:  Jupiter
Observer Point: Detroit, Mi
Center geodetic : 275.606400,33.7528000,1.581E-12 
{E-lon(deg),Lat(deg),Alt(km)}
Point in Time
2014-Sep-14 00:00
Target
Right Ascension:  09 00 45.15
Declination:      +17 29 55.9
Azimuth:      311.279081
Elevation:        -22.059168

Below is the calculation I use to get the zenith latitude. I'm not 100% sure it's correct, but logically it seems correct.

Calulate Target Zenith Latitude
90° - Elevation = Target Zenith Latitude
90° - (-22°) = 78°
- The Zenith latitude of Jupiter at 
  this point in time is 78° 

Though I'm not sure how to find the correct zenith longitude. Since the elevation is negative, -22.059168, I assume Jupiter is not visible from Detroit as it is below the horizon at that time.

Logically it seems that one would first need to find the NWN horizon longitude at the Azimuth of 311. And then add 22° to that latitude to find the zenith longitude of Jupiter to Earth given the data from the Detroit observer point..

I don't know if this logic is correct. And if it is, I don't know how to calculate that absolute horizon longitude at the 311 azimuth at a particular moment in time.

I'm hoping someone can help me calculate the correct longitude of the Earth's Zenith point to Jupiter.

Answer 1

Let's call the point on the Earth's surface where Jupiter is at the zenith as the "sub-Jupiter point".

The latitude of the sub-Jupiter point is equal to the declination of Jupiter. Therefore, latitude 17° 29' 55.9" N is the answer for the latitude. (This is true if the flattening of the Earth can be ignored which would introduce an error of a fraction of a degree.)

There are two methods of calculating the longitude of the sub-Jupiter point:

• Based on the known right ascension of Jupiter and the Greenwich Mean Sidereal Time (GMST, which can be calculated from the date and time), calculate the longitude where the meridian has the same right ascension as Jupiter.

The GMST can be calculated from a number of posts such as How to find Greenwich Mean Sideral Time?. If the GMST were 2 hour, then the meridian at 15 E longitude would be 3 hour right ascension, 15 W longitude would be at 1 hour right ascension, and so on. (1 hour right ascension = 15 degrees of longitude.) The difference between GMST and Jupiter's right ascension (converted from hours to degrees as needed) gives the longitude of the sub-Jupiter point.

• Based on the observed altitude and/or azimuth of Jupiter, calculate how many degrees Jupiter is east or west of the meridian. This is the hour angle of Jupiter. Add the hour angle to the current longitude to get the longitude of the sub-Jupiter point.

The relationship between celestial coordinates (hour angle and declination) and horizon coordinates (altitude and azimuth) can be calculated using posts such as Translating a zenith position to the nadir. This formula from that page gives the hour angle H directly (positive to the west of the observer): $$\mathrm{tan}\ H = \frac {\mathrm{sin}\ A}{\mathrm{cos}\ A\ \mathrm{sin}\ \phi + \mathrm{tan}\ h\ \mathrm{cos}\ \phi }$$

where $\phi$ is the latitude (+ North, − South) in degrees, A is the azimuth in degrees, and h is the height of the object in degrees. Please note that azimuths in this formula are measured from the South heading East (90°) then North (180°) and West (270° = -90°).

If H were +30 degrees, then the sub-Jupiter longitude would be 30 degrees further west of the observer's longitude.

Answer 2

@JohnHoltz provided an excellent answer above for professional scientists. To an amateur like myself though it's like learning a new language. So I did some research on dual observation coordinates each with different azimuths pointing toward Jupiter at a given point in time. With the help of this excellent website I was able to find the intersecting point of two azimuths from two separate Earth coordinates.

I found that by using the triangulation on the above link that the resulting coordinates were a very-near exact point of Jupiter's zenith on Earth. (I use the word zenith to mean the spot on Earth where Jupiter is currently at 90 degrees directly above at all possible azimuths and in the very center of the night sky.)

Below is the data I used for date and 2 locations

'   Data Date:  01-01-1970 00:00:00
'   Detroit Lat = 42°19'53.1"N
'   Detroit Lon = 83°02'44.7"W
'   Detroit Lat = 42.331429
'   Detroit Lon = -83.045753
'   Detroit Azimuth to Jupiter = 335.578522
'   Antarctic Lat = 67°12'50.2"S  - Randomly chosen from the other side of earth
'   Antarctic Lon = 88°02'30.3"E
'   Antarctic Lat = -67.21394
'   Antarctic Lon = 88.04174
'   Antarctic Azimuth to Jupiter = 26.025038

When I plug the two coordinates into the webpage the resulting coords are:

Intersection point: 11° 08′ 50″ S, 110° 19′ 49″ E

Now you can plug in the intersecting coordinates into Nasa's Jet Propulsion Horizon App to find azimuth and elevation of that location: (First I converted DDMMMSS to decimal coords)

Ephemeris Type: Observer Table
Target Body: Jupiter
Observer Location -> Specify Coordinates
Lon:   110.330277777778  (Converted from 110° 19′ 49″ E)
Lat:  -11.1472222222222  (Converted from 11° 08′ 50″ S)
Date Start Time:  1970-01-01
Date Stop Time:   1970-01-02
Step Size3 "1" and Type "hours"
Finally click "Generate Ephemeris"
Results:  Azimuth 182.190765, Elevation 89.999201

As you can see the intersection coordinates are nearly 180° Azimuth and 90° elevation, which would be Jupiter's zenith point in my words. At that location Jupiter would be (almost) perfectly centered in the sky. I used the words "night sky" before but that was wrong. That is where Jupiter will be at that particular date and time, regardless if it is night or day at that location.

The website creator provides both mathematical triangulation formulas and also his own JavaScript code for use in web pages. But as an amateur astronomer I'm having some difficulty translating the formula math to Microsoft Office VBA programming scripting language. I will keep chugging along at that. If any VBA programmers out here are more knowledgeable on the math of it all then I would certainly accept any help on that.

In the mean time I did manage to create a set of VBA coordinate conversion functions if anyone is interested in them.

Sub test_DMS_Coordinates_To_Decimal()
Dim myCoordString As String
myCoordString = DMS_Coordinates_To_Decimal(Sheet1.Range("H8").Value)


Debug.Print Sheet1.Range("H8").Value & " = " & myCoordString


End Sub
Function DMS_Coordinates_To_Decimal(dmsCoords As String) As String
'/////////////////////////////////////////////////////////////////////////////////////////////////////
'
'   This macro is built to accept any of the following formats of DMS coordinates to convert to decimal
'
'   38° 53' 55" N
'   38°53'55"N
'   38 53 55 N
'
'   USAGE:  anyStringVariable = DMS_Coordinates_To_Decimal(sheet1.Range("G2").value)
'      OR:  anyStringVariable = DMS_Coordinates_To_Decimal("38 53 55 N")
'      OR:  anyStringVariable = DMS_Coordinates_To_Decimal("38° 53' 55" & chr(34) & "N")
'
'
        Dim degreesString As String
        Dim minutesString As String
        Dim secondsString As String
        Dim finalProduct1 As String
        Dim finalProduct2 As String
        Dim finalProduct3 As String
        Dim degreesBooleanStart As Boolean
        Dim degreesBooleanStop As Boolean
    Dim minutesBooleanStart As Boolean
    Dim minutesBooleanStop As Boolean

    Dim secondsBooleanStart As Boolean
    Dim secondsBooleanStop As Boolean

    degreesBooleanStop = False
    minutesBooleanStop = True
    secondsBooleanStop = True


    For i = 1 To Len(dmsCoords)


getDegrees:
        If Not degreesBooleanStop And IsNumeric(Mid(dmsCoords, i, 1)) Then
            degreesBooleanStart = True
            degreesString = degreesString & CStr(Mid(dmsCoords, i, 1))
        Else
            If degreesBooleanStart And Not degreesBooleanStop Then
                degreesBooleanStop = True
                minutesBooleanStop = False
                GoTo getMinutes
            End If
        End If


getMinutes:
        If Not minutesBooleanStop And IsNumeric(Mid(dmsCoords, i, 1)) Then
            minutesBooleanStart = True
            minutesString = minutesString & CStr(Mid(dmsCoords, i, 1))
        Else
            If minutesBooleanStart And Not minutesBooleanStop Then
                minutesStart = i + 1
                minutesBooleanStop = True
                secondsBooleanStop = False
                GoTo getSeconds
            End If
        End If


getSeconds:
        If Not secondsBooleanStop And IsNumeric(Mid(dmsCoords, i, 1)) Then
            secondsBooleanStart = True
            secondsString = secondsString & CStr(Mid(dmsCoords, i, 1))
        Else
            If secondsBooleanStart And Not secondsBooleanStop Then
                secondsBooleanStop = True
                GoTo do_The_Math
            End If
        End If
    Next i


do_The_Math:
    finalProduct1 = degreesString
    finalProduct2 = finalProduct & CStr((CDbl(minutesString) / 60) + CDbl(secondsString) / 3600)

    finalProduct3 = CStr(CDbl(finalProduct1) + CDbl(finalProduct2))

    If InStr(1, UCase(dmsCoords), "S") > 0 Or InStr(1, UCase(dmsCoords), "W") > 0 Then
        finalProduct3 = CStr(CDbl(finalProduct3) * -1)
    End If

    DMS_Coordinates_To_Decimal = finalProduct3


End Function
Sub TEST_convert_Decimal_To_Degrees_Minutes_Seconds()
Debug.Print convert_Decimal_To_Degrees_Minutes_Seconds(-67.21394, -88.04174)


End Sub
Function convert_Decimal_To_Degrees_Minutes_Seconds(ddLat As Double, ddLong As Double) As String
'        Dim ddLat As Double
'        Dim ddLon As Double
    Dim dmsLatDeg As Long
    Dim dmsLatMin As Long
    Dim dmsLatSec As Double
    Dim dmsLatHem As String ' "N" or "S"

    Dim dmsLongDeg  As Long
    Dim dmsLongMin  As Long
    Dim dmsLongSec  As Double
    Dim dmsLongHem  As String ' "E" or "W"

    Dim myLatSplitArr
    Dim myLonSplitArr
    Dim myLatMinuteRemainderSplitArr
    Dim myLonMinuteRemainderSplitArr
    Dim myLatSecondRemainderSplitArr
    Dim myLonSecondRemainderSplitArr

    Dim ddLatMinuteRemainder As Double
    Dim ddLonMinuteRemainder As Double

    Dim ddLatSecondRemainder As Double
    Dim ddLonSecondRemainder As Double


    ' Is it negative number
    If ddLat < 0 Then          ' if decimal lat has a negative sign "-""  Set South option "1"
        dmsLatHem = "S"
    Else
        dmsLatHem = "N"
    End If

    If ddLong < 0 Then
        dmsLongHem = "W"
    Else
        dmsLongHem = "E"
    End If


    '////////   Degrees
    myLatSplitArr = Split(ddLat, ".")    '  Split Lat decimal by periiod (.)
    dmsLatDeg = myLatSplitArr(0)            '  ddLatDeg = part before period

    myLonSplitArr = Split(ddLong, ".")   '  Split Lon decimal by periiod (.)
    dmsLongDeg = myLonSplitArr(0)           '  dmsLongDeg = part before period


    '///////    Minutes
    ddLatMinuteRemainder = CDbl("0." & CStr(myLatSplitArr(1))) * 60
    myLatMinuteRemainderSplitArr = Split(ddLatMinuteRemainder, ".")
    dmsLatMin = myLatMinuteRemainderSplitArr(0)

    ddLonMinuteRemainder = CDbl("0." & CStr(myLonSplitArr(1))) * 60
    myLonMinuteRemainderSplitArr = Split(ddLonMinuteRemainder, ".")
    dmsLongMin = myLonMinuteRemainderSplitArr(0)


    '////////   Seconds
    ddLatSecondRemainder = CDbl("0." & CStr(myLatMinuteRemainderSplitArr(1))) * 60
    dmsLatSec = Round(ddLatSecondRemainder, 1)
    ddLonSecondRemainder = CDbl("0." & CStr(myLonMinuteRemainderSplitArr(1))) * 60
    dmsLongSec = Round(ddLonSecondRemainder, 1)

    myString = Replace(dmsLatDeg & " " & dmsLatMin & " " & dmsLatSec & " " & dmsLatHem & vbNewLine, "-", "") & _
                Replace(dmsLongDeg & " " & dmsLongMin & " " & dmsLongSec & " " & dmsLongHem, "-", "")

    connvert_Decimal_To_Degrees_Minutes_Seconds = myString



End Function

Answer 3

As promised I am including the VBA code below for finding the coordinates of the intersecting point of two separate azimuths.

I adapted this VBA code from Chris Veness' calculation functions found at http://www.movable-type.co.uk/scripts/latlong.html. This is a must have sight for any serious astronomer or gravitational scientist.

I only tested a few random points with this code but so far have not found any errors. My only concern is that sometiems the longitude comes back as a double precision number (stripped of neg or pos by abs) larger than 180 degrees. Of course we have to manipulate that by subtracting it from 360. And the problem arises when the lon is negative and or positive. I'm certain my code addresses it properly when the longitude is negative and higher than 180. But I don't know if it works properly when the longitude is > 180 and is positive. I'd appreciate it if yu let me know if you find an answer to that. I ran this code to calculate the geographic subpoint of jupiter through one year of data broken down hourly in about 4 seconds. That's 8,761 calculations made in under 5 seconds. The bearings are from two separate coordinates on Earth with azimuth bearing pointing at the same extraneous body of mass. In this case Jupiter. The resulting coordinates usually prove to be 89.99 or higher elevation on the given date and time.

Function getIntersectCoords(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double, brng1 As Double, brng2 As Double) As String
'   This code finds the geographic intersecting coordinate points
'   where 2 azimuth bearings cross paths
'   This VBA code was adapted from Chris Veness' calculation functions
'   at http://www.movable-type.co.uk/scripts/latlong.html
'   It Was Adapted to VBA by John Richter August 2022.
'   https://johnallenrichter.wordpress.com/
        Dim PI As Double
        PI = WorksheetFunction.PI

        Dim q1 As Double
        Dim y1 As Double
        Dim q2 As Double
        Dim y2  As Double

        q1 = degrees_to_radians(lat1)
        y1 = degrees_to_radians(lon1)
        q2 = degrees_to_radians(lat2)
        y2 = degrees_to_radians(lon2)

        Dim O13 As Double
        Dim O23 As Double
        Dim change_q As Double
        Dim change_y As Double

        O13 = degrees_to_radians(brng1)
        O23 = degrees_to_radians(brng2)
        change_q = q2 - q1
        change_y = y2 - y1

        Dim d12 As Double


'           angular dist. p1–p2,  Value with current positions and bearings = 2.69916649962162
            d12 = 2 * WorksheetFunction.Asin(SqRt(Sin(change_q / 2) * Sin(change_q / 2) + Cos(q1) * Cos(q2) * Sin(change_y / 2) * Sin(change_y / 2)))
        Dim cosOa As Double
        Dim cosOb As Double
        Dim Oa As Double
        Dim Ob As Double

        cosOa = (Sin(q2) - Sin(q1) * Cos(d12)) / (Sin(d12) * Cos(q1))
        cosOb = (Sin(q1) - Sin(q2) * Cos(d12)) / (Sin(d12) * Cos(q2))

        Oa = WorksheetFunction.Acos(WorksheetFunction.Min(WorksheetFunction.Max(cosOa, -1), 1)) '      protect against rounding errors
        Ob = WorksheetFunction.Acos(WorksheetFunction.Min(WorksheetFunction.Max(cosOb, -1), 1)) '      protect against rounding errors


        Dim O12 As Double
        If Sin(y2 - y1) > 0 Then
            O12 = Oa
        Else
            O12 = 2 * PI - Oa
        End If

        Dim O21 As Double
        If Sin(y2 - y1) > 0 Then
             O21 = 2 * PI - Ob
        Else
            O12 = Ob
        End If



        Dim a1 As Double
        Dim a2 As Double

        a1 = O13 - O12
        a2 = O21 - O23


        Dim cosa3 As Double

        cosa3 = -Cos(a1) * Cos(a2) + Sin(a1) * Sin(a2) * Cos(d12)

        Dim d13 As Double
        d13 = WorksheetFunction.Atan2(Cos(a2) + Cos(a1) * cosa3, Sin(d12) * Sin(a1) * Sin(a2))

        Dim q3 As Double
        q3 = WorksheetFunction.Asin(WorksheetFunction.Min(WorksheetFunction.Max(Sin(q1) * Cos(d13) + Cos(q1) * Sin(d13) * Cos(O13), -1), 1))



        Dim change_y13 As Double
        change_y13 = WorksheetFunction.Atan2(Cos(d13) - Sin(q1) * Sin(q3), Sin(O13) * Sin(d13) * Cos(q1))


        Dim y3 As Double
        y3 = y1 + change_y13

        Dim lat As Double
        Dim lon As Double

        lat = radians_to_degrees(q3)
        lon = radians_to_degrees(y3)

        If Abs(lon) > 180 Then
            If lon < 0 Then
                lon = 360 - Abs(lon)
            Else
                lon = -1 * (360 - Abs(lon))
            End If
        End If


        getIntersectCoords = CStr(lat) & "," & CStr(lon)





End Function
Sub testgetIntersectCoords()
Dim myString As String
Dim mySplitArray
    ' bearings are on sheet1 columns B and C
    ' resulting coordinates are put into columns D and E of the same row

    ' The pts Lat1, lon1 are constant because I used the same pts for all azimuth bearings

    Const lat1 As Double = 42.331429        '   Position 1 Lat (Detroit)
    Const lon1 As Double = -83.045753       '   Position 1 Lon (Detroit)
    Const lat2 As Double = -67.21394        '   Position 2 Lat (Antarctic)
    Const lon2 As Double = 88.04174         '   Position 2 Lon (Antarctic)
    Dim brng1 As Double
    Dim brng2 As Double

    For i = 2 To Sheet1.Range("A" & Rows.Count).End(xlUp).Row
        brng1 = CDbl(Sheet1.Range("B" & i).Value)
        brng2 = CDbl(Sheet1.Range("C" & i).Value)

        myString = getIntersectCoords(lat1, lon1, lat2, lon2, brng1, brng2)

        mySplitArray = Split(myString, ",")
        latitude = mySplitArray(0)
        longitude = mySplitArray(1)
        Sheet1.Range("D" & i).Value = latitude
        Sheet1.Range("E" & i).Value = longitude
    Next i


End Sub

Answer 4

I found a slight anomaly with the VBA conversion code above. This will detail how I plan to work around the anomaly.

Methodology: My goal was to find and track the geographic subpoint coordinates over time of several of our Solar system's planets, the Sun, and the Moon. These subpoints are sometimes referred to as "zeniths." At least by me. It's my process. They are simply a coordinate point on Earth which at a given point in time - exact second in time, that the target planet or body will be at 90° elevation. This means the observer at those coords at that specific time will find the target planet directly in the center of his observable sky. And he, trees or buildings around him will be the closest things in distance to that target planet than any other point on Earth. In these example's I'm using Jupiter as the target planet. I have used two definite locations on Earth to derive the azimuth directions to Jupiter. One coord point is Detroit and the second is in the eastern hemisphere in the Antarctic. I think the exact coordinates are listed above somewhere. The important thing to remember for now is this: The bearing between the Detroit coords and the Antarctic coords is 171.94722222°. And conversely the bearings between the Antarctic coordinates and Detroit is -195.516388888889.

I downloaded an entire year of azimuths pointing toward Jupiter in hourly increments from both coordinates in Detroit and in the Antarctic. (See Nasa or jet propulsion lab Horizons API for getting azimuths) This resulted in over 8,500 rows of azimuth data on an Excel spreadsheet. I had hypothesized that the intersecting point of the two separate azimuth bearings could result in fairly accurate geographic subpoint coords for any target planet or body. The hypothesis was overwhelmingly correct. Using formulas devised by Chris Veness at this website I was able to adapt his formulas to MS VBA scripting code. (See above) I used the above code to calculate the zenith subpoints for all 8,500 sets of azimuths. Then I used the horizon app on each of the 8,500 subpoint coordinates to find the elevation point for each one. In all but 10 coords the resulting elevations were over 89.9° for the specific time given.

Why did 10 coords fail to produce accurate elevations?

Earlier I said it is important to remember the bearings between the two original coordinate points Detroit and Antarctic. As it turns out, whenever the planet's calculated coordinate pair, lat and lon, are similar to the inter coordinate bearings between Detroit and Antarctic, the sub point position coordinates will be wrong. Basically, whenever the actual subpoint is very near one of the fixed original coordinates, the calculations have a hiccup. As I said it happened in 10 out of over 8,500 times. Seemingly small percentage but very significant for my purposes.

An eyeball estimate is that if a resulting latitude falls within one degree of the Detroit to Antarctic bearing, 171.??????, and the resulting longitude is within 1 degree of the Antarctic to Detroit bearing, -195.???, the resulting sub-point position coordinates will be wrong. In the ten errors I found elevations between 50° and 88.6°.

The fix?

The fix is very simple. While the above VBA code is calculating sub-position coordinates, when it encounters any lat within 1 degree of 171.??? (the first original location) and/or long of -195.???? (bearings to the second coord in Antarctic) then I am going to add a function to call an xmlHTTP request to horizon for azimuths from two completely different locations other than Detroit and the Antarctic. It would have to be polar opposites of Detroit and Antarctic. Perhaps a point in Brazil for the western hemisphere and a point in Poland for the eastern hemisphere. Then calculating those coordinates with the proper observation coordinates should produce sub-point coordinates that result in 89.9° elevation. Theoretically that should result in 100% reliable calculations of planetary/body geographic sub-points.