3
$\begingroup$

I need to find the two times of zenith passage of the Galactic Equator for a given Earth Lat/Lon on any day between 3500 BC and 1500 AD. Since the Galactic Equator is a line and not a point, I try to use PyEphem to create a zenith point for my geographic location (see code below), then find the Galactic coordinates of that point. When the Galactic latitude of this zenith point is zero, that's the time of zenith passage of the Galactic Equator.

However, I think I need my zenith point to be in Equatorial coordinates (not Ecliptic) to accomodate the geographic location and the date/time, but I can't figure out how to make that work.

Here is my Python script based on PyEphem: 

import ephem
from math import *

D = (input("Enter date (yyyy/mm/dd): "))
d = ephem.Date(D)
print(d)

location = ephem.Observer()
location.lon = '-87.0'
location.lat = '20.0'
location.elevation = 12
location.pressure = 1013.0
import csv
my_list = []
with open(r'C:\Users\Mark\Desktop\PYTHON\gzt.csv', 'w', newline='') as f:
        for T in range(1,1440):
                dm = ephem.Date(d + ephem.minute*T)
                location.date = dm
                ra, dec = location.radec_of('0', '90')
                ec = ephem.Ecliptic(ra, dec, epoch = d)
                eq = ephem.Equatorial(ec)
                #print(eq.ra, eq.dec)
                z = ephem.FixedBody()
                z._ra = eq.ra
                z._dec =eq.dec 
                #z._epoch = d
                z.compute()
                #print(z.ra, z.dec)
                zg = ephem.Galactic(z.a_ra, z.a_dec)
                dmt = (('%s') % (dm))
                #print(dmt, 'Zlat = ', zg.lat)
                AZ = (('%s  %s') % (dmt, zg.lat))
                my_list = [AZ]
                writer = csv.writer(f, dialect='excel-tab')
                writer.writerow(my_list)
Improve this question
$\endgroup$
3
  • $\begingroup$ I cannot help you with the programming, but my suggestion would be to find the declination of the zenith (it equals the latitude). Then find the Right Ascension (RA) of the galactic equator where it cross the declination. Then find when that RA is at the zenith, that is, when the local sidereal time equals the RA. That will give you the 0, 1, 2, or 3 times in a day when the galactic equator is at the zenith. $\endgroup$
    – JohnHoltz
    Apr 10, 2020 at 12:38
  • $\begingroup$ Thank you! James K (I think) for fixing my post. $\endgroup$
    – Mark
    Apr 10, 2020 at 16:39
  • $\begingroup$ Thanks JohnHoltz. It's more involved than I imagined. I see now I must at the least pass the hour angle in Sidereal Time from each time step to the Equatorial coordinates. So I will try to adjust my code to do that part now. I'm not an Astronomer or a Programmer, but I need both of those fields for my research, so I teach myself what I can. But there is a "Swiss Cheese" effect to that. $\endgroup$
    – Mark
    Apr 10, 2020 at 17:09

1 Answer 1

Reset to default
2
$\begingroup$

If I eliminate ec and z, and do this:

    eq = ephem.Equatorial(ra, dec, epoch = d)
    zg = ephem.Galactic(eq)

the crossings I find in gzt.csv agree with Stellarium within a minute of time.

If you're open to another package by the same developer, skyfield.almanac can help with such calculations. In this case you would implement a function which takes a time argument and returns an integer indicating which side of the galactic equator the zenith is on. Then you would pass this function as an argument to find_discrete.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Mike, Very much appreciated! That simple change (and getting rid of 'epoch = d') makes this work just as I need. You saved my fanny from hours of looking into other things. I know of Skyfield, but have been resisting learning a new system. The 'find_discreet' seems like just what I need next though. $\endgroup$
    – Mark
    Apr 10, 2020 at 22:24

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .