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.