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)]
latsub = np.radians(dec)
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)]
tilt = halfpi - np.radians(dec)
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()