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

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.

• Is this something you need to do manually, or are tools acceptable? I'm pretty sure the skyfield python library has methods that could do this calculation from the time alone. Jul 18, 2022 at 20:16
• I'm actually hoping for a manual calculation that I can build into an Excel script function. I'm a MS Office VBA programmer and am not familiar with Python. I was unable to find any online tools that provide planetary zenith points by time and date. Timeanddate(dot)com offers sun and moon zeniths. Meanwhile I'll try to establish at least one zenith point and work backwards to my single location point to formulate a workable equation. Thanks Jul 19, 2022 at 13:00
• There's some issues with your data, the lat/lon specified is in Texas, not Detroit. And the az/elevation don't match for either of those places. Jul 19, 2022 at 15:50
• Google says it's in Atlanta, Georgia. goo.gl/maps/hE1QDFxMf8hkhLjv6 Jul 20, 2022 at 2:45
• @JohnMuggins the term for the point you're looking for is called the Geographic Position, that might help you find more information about it. It is most commonly used in Celestial Navigation. Jul 30, 2022 at 15:40

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.

• Thanks @JohnHoltz . When you say 'The latitude of the sub-Jupiter point is equal to the declination of Jupiter..' Are you talking about the declination of Jupiter from the observation coordinates or from the center of the Earth? Thank you kndly for this information. I'm going to take it home this weekend and study it. It's a little over my usual forte, but I'm sure to get it. I was toying with the idea of providing two observation points and calculating the point at which both their azimuth lines crossed as the sub-Jupiter Point. Jul 22, 2022 at 18:09
• @JohnMuggins Jupiter is far away, so far that the declination seen from the center of the Earth and from the surface of the Earth is the same. Jul 22, 2022 at 23:04
• Thank you for the advice. I'm neither a mathematician nor an astronomer and after studying your solutions it's apparent I can't do this. My solution will most likely be to use 2 observer locations and to find where their azimuth points to Jupiter converge. Say ObsPoint1(177.08437°E, 68.21444°N, 0 m, azimuth to Jupiter 251.448450) and ObsPoint2(-16.61562°E, -68.06560°N, azimuth to Jupiter 121.179359) Can you tell me how to find the coordinates where those 2 lines intersect? Jul 24, 2022 at 23:35
• The problem with your new approach of using two observers is that their observations intersect at about 500 million miles away (at Jupiter). Where the azimuths intersect is not a solution. Solving the one equation provided for H and adding it to the observer's longitude is an easy solution to program. Of course, this assumes you have the information that you provided: declination of Jupiter (one of 2 coordinates to locate an object on a star chart), observer's latitude and longitude, observed azimuth and altitude of Jupiter. If you don't have those, it takes more programming and math. Jul 25, 2022 at 1:24
• Thank you for your help @JohnHoltz. But I think there might be some confusion due to the way I worded my question. As an amateur I'm, certain that my wording was probably confusing. I hope my answer below will help the understanding. Jul 28, 2022 at 23:15

@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

• This is a good answer, but it answers a different question than your original question. Your original question provided the declination of Jupiter ("Declination: +17 29 55.9") and ONE observing location. That input is different than having observers at TWO (widely spaced) locations. In other words, you cannot use this solution with the information provided in the question. I suggest you create a new question based on two observing locations and cut/paste your answer to your new question. Jul 29, 2022 at 3:22
• Hi @JohnHoltz. Again, please forgive my words. I'm an amateur researcher of gravity and am not any kind of mathemetician. I checked your answer as the favored answer. I don't mean to upset the protocols. But what I really need is a way to calculate the zenith spot on Earth directly facing Jupiter at 90 degrees over time. What I found is that 2 intersecting azimuths is the easiest method for me due to my complete ignorance of the math symbols and methods involved with that. It's my process. Just simply how my brain works. Thank you for helping me. Jul 29, 2022 at 17:46

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

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

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

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


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.