# Star visibility above horizon map likewise day & light map or moonlight world map

I am looking for an option to present the given star visibility world map likewise day & night map or moonlight world map like below:  Question: Is it possible to generate a map like this for some particular star or at least the point in the sky defined by celestial coordinates?

As long as we can live with a few tens of kilometers accuracy, we can treat the Earth's surface as a sphere.

The subpoint of the object (the point on Earth at which the object appears at the zenith, overhead, or the central point on Earth seen from the view of the object if they could see Earth clearly) is simply given by

$$\text{lat} = \text{Dec}$$

$$\text{lon} = \text{R.A.} - \text{GST}$$

where GST is Greenwich Sidereal time

To make your plot we draw a great circle around the subpoint, then make a simple Equirectangular projection or lat/lon projection.

In Python for example, calculate the latitude and longitude of the subpoint,

latsub = np.radians(dec)
lonsub = np.radians(RA - GST) # where these have already been converted to degrees


Pick a bunch of points to make a circle

theta = np.linspace(0, 2*np.pi, 10001)


Draw a unit circle around Earth's equator, assuming the subpoint is at the north pole (dec = +90 degrees).

x, y, z = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]


Tilt it by 90 - dec

tilt = np.pi/2 - np.radians(dec)
x, z = [x * f(tilt) for f in (np.cos, np.sin)]


calculate the latitude and longitude of the tilted circle by shifting the longitude (subtracting the longitude of the subpoint)

loncirc = np.arctan2(y, x) - lonsub
latcirc = np.arctan2(z, np.sqrt(x**2 + y**2))


Put it all together and convert units properly, and you get something like this. I did it for the Sun but you can do this for any object for which you have the R.A. and Dec. import numpy as np
import matplotlib.pyplot as plt

# for right now, from Google:
RA = (3 + 29/60 + 38/3600) * 360 / 24 # deg for the Sun
dec = +19. # deg for the Sun
GST = (16 + 15/60) * 360 / 24 # deg

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

lonsub = np.radians(RA - GST) # where these have already been converted to degrees

theta = np.linspace(0, twopi, 10001)
x, y, z = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]

x, z = [x * f(tilt) for f in (np.cos, np.sin)]

loncirc = np.arctan2(y, x) - lonsub
loncirc = np.mod(loncirc + pi, twopi) - pi
latcirc = np.arctan2(z, np.sqrt(x**2 + y**2))

dloncirc = np.abs(loncirc[1:] - loncirc[:-1])
jump = np.argmax(dloncirc)

latcirc, loncirc = [thing[:-1] for thing in (latcirc, loncirc)]
loncirc[jump] = np.nan
latcirc[jump] = np.nan

locd, lacd = [np.degrees(thing) for thing in (loncirc, latcirc)]
losd, lasd = [np.degrees(thing) for thing in (lonsub, latsub)]
losd = np.mod(losd + 180, 360) - 180 # put it within the maps -180 to +180 range

plt.plot(locd, lacd)
plt.xlim(-182, 182)
plt.ylim(-90, 90)
x, y = 121.5, 25 # Taipei
plt.plot([x], [y], 'or')
plt.text(x, y, ' taipei')
plt.plot([-losd], [lasd], 'oy')
plt.text(-losd, lasd, ' sun')

plt.title('Sunlit Earth 08:40 local time TPE, sunrise was 5:09')
plt.gca().set_aspect('equal')
plt.show()

• I don't understand this pattern completely. It's been demonstrated on the Sun. How about the Sirius star for example?
– MKR
Jul 2 at 11:11
• OK, I got it a bit. The biggest uncertainty is the declination. dec = +23. # deg for the Sun how will I know that declination is +23,37 ie? It looks like only integers come here. I need the declination more detailed. Could it be possible?
– MKR
Jul 2 at 11:17
• no rush. Thank you in advance.
– MKR
Jul 2 at 13:31
• are you able to help? Thank you in advance
– MKR
Nov 18 at 15:00