mirror of
https://codeberg.org/Yael-II/Astrobs-Tools.git
synced 2026-03-15 11:26:27 +01:00
408 lines
14 KiB
Python
Executable File
408 lines
14 KiB
Python
Executable File
"""
|
|
* SCOPE: Sky Coordinates for Observations Python Estimator
|
|
* Version 2.1 - 2024-09-21
|
|
@ Yaël Moussouni
|
|
@ Observatory of Strasbourg
|
|
"""
|
|
import astropy.time as time
|
|
import astropy.coordinates as coord
|
|
import astropy.units as u
|
|
|
|
YES = ["YES", "Y", "1"]
|
|
OUT_DIR = "./Output/"
|
|
# Definitions
|
|
def get_sun(sun_time, sun_frame, angle, reverse=False):
|
|
i = 0
|
|
while coord.get_sun(sun_time).transform_to(sun_frame).alt.value > angle:
|
|
if reverse:
|
|
sun_time -= 1 * u.min
|
|
else:
|
|
sun_time += 1 * u.min
|
|
i += 1
|
|
if i > 60*24:
|
|
print("\033[93m"
|
|
+ "Warning: Sun did not reached the expected angle."
|
|
+ "\033[0m")
|
|
return time.Time("1970-01-01T12:00") # Error date
|
|
return sun_time
|
|
|
|
# Locations
|
|
available_obs = ["Observatoire astronomique de Strasbourg (ObAS)",
|
|
"Observatoire de Haute-Provence (OHP)"]
|
|
available_obs_short = ["ObAS",
|
|
"OHP"]
|
|
available_obs_lat = [coord.Latitude("""48d34m59s"""),
|
|
coord.Latitude("""43d55m53s""")]
|
|
available_obs_lon = [coord.Longitude("""07d46m07s"""),
|
|
coord.Longitude("""5d42m45s""")]
|
|
|
|
|
|
def main():
|
|
try:
|
|
print("\033[32m"+ "Please select a location:"+"\033[0m")
|
|
for i in range(len(available_obs)):
|
|
print("\033[32m"+ "\t{}. ".format(i+1)+"\033[0m"+available_obs[i]+".")
|
|
selection = int(input("\033[32m"+ "Choice: "+"\033[0m"))-1
|
|
if selection >= len(available_obs):
|
|
selection = -1
|
|
print("\033[93m"
|
|
+ ("Warning: Location unavailable, "
|
|
"selected the last one instead.")
|
|
+ "\033[0m")
|
|
except:
|
|
selection = -1
|
|
print("\033[91m"
|
|
+ ("Error! Input value not recognized, "
|
|
"selected the last location instead")
|
|
+ "\033[0m")
|
|
|
|
obs_name = available_obs[selection]
|
|
obs_short = available_obs_short[selection]
|
|
obs_lat = available_obs_lat[selection]
|
|
obs_lon = available_obs_lon[selection]
|
|
obs_loc = coord.EarthLocation(lon=obs_lon, lat=obs_lat)
|
|
|
|
# Date and time
|
|
try:
|
|
obs_date = input("\033[32m"
|
|
+ "Date of observation [YYYY-MM-DD]: "
|
|
+ "\033[0m")
|
|
test_date = time.Time(obs_date)
|
|
except:
|
|
obs_date = time.Time("1970-01-01", format="iso").now()
|
|
obs_date = obs_date.to_string()[0:10]
|
|
print("\033[91m"
|
|
+ "Error! Input value not recognized, selected today instead"
|
|
+ "\033[0m")
|
|
else:
|
|
del test_date
|
|
print("\033[34m" + "Computing Sun position..." + "\033[0m")
|
|
sun_frame = coord.AltAz(location=obs_loc)
|
|
sun_time = time.Time(obs_date + " 12:00:00", scale="utc")
|
|
|
|
sun_set = get_sun(sun_time, sun_frame, 0)
|
|
sun_civil = get_sun(sun_set, sun_frame, -6)
|
|
sun_nautical = get_sun(sun_civil, sun_frame, -12)
|
|
sun_astronomical = get_sun(sun_nautical, sun_frame, -18)
|
|
|
|
sun_time = time.Time(obs_date + " 12:00:00", scale="utc") + 1 * u.day
|
|
|
|
sun_rise = get_sun(sun_time, sun_frame, 0, reverse=True)
|
|
sun_rcivil = get_sun(sun_rise, sun_frame, -6, reverse=True)
|
|
sun_rnautical = get_sun(sun_rcivil, sun_frame, -12, reverse=True)
|
|
sun_rastronomical = get_sun(sun_rnautical, sun_frame, -18, reverse=True)
|
|
print("\033[34m" + "done!" + "\033[0m")
|
|
|
|
print("\033[36m"
|
|
+ "Night time [UTC]:"
|
|
+ "\033[0m")
|
|
print("\033[36m"
|
|
+ "\t === Evening === "
|
|
+ "\033[0m")
|
|
print("\033[36m"
|
|
+ "\t sunset: "
|
|
+ "\033[0m"
|
|
+sun_set.to_string())
|
|
print("\033[36m"
|
|
+ "\t civil twilight: "
|
|
+ "\033[0m"
|
|
+sun_civil.to_string())
|
|
print("\033[36m"
|
|
+ "\t nautical twilight: "
|
|
+ "\033[0m"
|
|
+sun_nautical.to_string())
|
|
print("\033[36m"
|
|
+ "\tastronomical twilight: "
|
|
+ "\033[0m"
|
|
+sun_astronomical.to_string())
|
|
print("\033[36m"
|
|
+ "\t === Morning === "
|
|
+ "\033[0m")
|
|
print("\033[36m"
|
|
+ "\tastronomical twilight: "
|
|
+ "\033[0m"
|
|
+sun_rastronomical.to_string())
|
|
print("\033[36m"
|
|
+ "\t nautical twilight: "
|
|
+ "\033[0m"
|
|
+sun_rnautical.to_string())
|
|
print("\033[36m"
|
|
+ "\t civil twilight: "
|
|
+ "\033[0m"
|
|
+sun_rcivil.to_string())
|
|
print("\033[36m"
|
|
+ "\t sunrise: "
|
|
+ "\033[0m"
|
|
+sun_rise.to_string())
|
|
|
|
|
|
obs_begin = input("\033[32m"
|
|
+ "Begin time of observation [hh:mm format; UTC]: "
|
|
+ "\033[0m")
|
|
obs_end = input("\033[32m"
|
|
+ " End time of observation [hh:mm format; UTC]: "
|
|
+ "\033[0m")
|
|
|
|
if obs_begin == "set" or obs_begin == "rise":
|
|
obs_begin_time = sun_set
|
|
elif obs_begin == "civil":
|
|
obs_begin_time = sun_civil
|
|
elif obs_begin == "nautical":
|
|
obs_begin_time = sun_nautical
|
|
elif obs_begin == "astronomical" or obs_begin == "auto" or obs_begin == "":
|
|
obs_begin_time = sun_astronomical
|
|
else:
|
|
try:
|
|
obs_begin_time = time.Time(obs_date + " " + obs_begin, scale="utc")
|
|
except:
|
|
obs_begin_time = sun_astronomical
|
|
print("\033[91m"
|
|
+ ("Error! Input value not recognized, "
|
|
"selected astronomical twilight instead.")
|
|
+ "\033[0m")
|
|
|
|
|
|
if obs_end == "set" or obs_end == "rise":
|
|
obs_end_time = sun_rise
|
|
elif obs_end == "civil":
|
|
obs_end_time = sun_rcivil
|
|
elif obs_end == "nautical":
|
|
obs_end_time = sun_rnautical
|
|
elif obs_end == "astronomical" or obs_end == "auto" or obs_end == "":
|
|
obs_end_time = sun_rastronomical
|
|
else:
|
|
try:
|
|
obs_end_time = time.Time(obs_date + " " + obs_end, scale="utc")
|
|
except:
|
|
obs_end_time = sun_rastronomical
|
|
print("\033[91m"
|
|
+ ("Error! Input value not recognized, "
|
|
"selected astronomical twilight instead.")
|
|
+ "\033[0m")
|
|
|
|
if obs_end_time.value < obs_begin_time.value:
|
|
obs_end_time += 1*u.day
|
|
|
|
obs_duration = obs_end_time - obs_begin_time
|
|
|
|
# Coordinates
|
|
|
|
print("\033[36m"
|
|
+ "Observation constraints"
|
|
+ "\033[0m")
|
|
try:
|
|
min_alt = coord.Angle(
|
|
input("\033[32m"
|
|
+ " Minimum altitude above horizon [deg]: "
|
|
+ "\033[0m")
|
|
+ "d")
|
|
except:
|
|
min_alt = coord.Angle(50 * u.deg)
|
|
print("\033[91m"
|
|
+ "Error! Input value not recognized, selected 50° instead."
|
|
+ "\033[0m")
|
|
try:
|
|
min_dec = coord.Angle(
|
|
input("\033[32m"
|
|
+ " Minimum declination around the north [deg]: "
|
|
+ "\033[0m")
|
|
+ "d")
|
|
except:
|
|
min_dec = coord.Angle(30 * u.deg)
|
|
print("\033[91m"
|
|
+ "Error! Input value not recognized, selected 30° instead."
|
|
+ "\033[0m")
|
|
try:
|
|
window_east = coord.Angle(
|
|
input("\033[32m"
|
|
+ "Observation window before crossing meridian [h]: "
|
|
+ "\033[0m")
|
|
+ " hours")
|
|
except:
|
|
window_east = coord.Angle(4 * u.h)
|
|
print("\033[91m"
|
|
+ "Error! Input value not recognized, selected 4h instead."
|
|
+ "\033[0m")
|
|
try:
|
|
window_west = coord.Angle(
|
|
input("\033[32m"
|
|
+ " Observation window after crossing meridian [h]: "
|
|
+ "\033[0m")
|
|
+ " hours")
|
|
except:
|
|
window_west = coord.Angle(2 * u.h)
|
|
print("\033[91m"
|
|
+ "Error! Input value not recognized, selected 2h instead."
|
|
+ "\033[0m")
|
|
|
|
obs_sky_rotation = coord.Angle(
|
|
"{} hours".format(obs_duration.to_value("h")))\
|
|
* (1*u.sday/u.day).decompose()
|
|
|
|
# * Date/time/location objects
|
|
obs = time.Time(obs_begin_time, location=obs_loc)
|
|
|
|
# * Local sidereal time
|
|
st = obs.sidereal_time("apparent") # ? apparent or absolute
|
|
|
|
# * Borders and corners coordinates
|
|
dec_upper = (90*u.deg-min_dec).wrap_at(360*u.deg)
|
|
dec_lower = (min_alt + obs_loc.lat - 90*u.deg).wrap_at(180*u.deg)
|
|
a_west = (st + obs_sky_rotation + window_west).wrap_at(24*u.h)
|
|
a_east = (st - window_east).wrap_at(24*u.h)
|
|
|
|
if a_east.degree <= a_west.degree:
|
|
comp_symb = "&&"
|
|
comp_text = "AND"
|
|
else:
|
|
comp_symb = "||"
|
|
comp_text = "OR"
|
|
# * Output
|
|
print("")
|
|
print("\033[36m"
|
|
+ "Location: "
|
|
+ "\033[0m"
|
|
+ obs_name)
|
|
print("\033[36m"
|
|
+ "\tlat: "
|
|
+ "\033[0m"
|
|
+ obs_lat.to_string())
|
|
print("\033[36m"
|
|
+ "\tlon: "
|
|
+ "\033[0m"
|
|
+ obs_lon.to_string())
|
|
print("\033[36m"
|
|
+ "Observation"
|
|
+ "\033[0m")
|
|
print("\033[36m"
|
|
+ "\tbegin date and time: "
|
|
+ "\033[0m"
|
|
+ obs_begin_time.to_string())
|
|
print("\033[36m"
|
|
+ "\t end date and time: "
|
|
+ "\033[0m"
|
|
+ obs_end_time.to_string())
|
|
print("\033[36m"
|
|
+ "\tbegin sidereal time: "
|
|
+ "\033[0m"
|
|
+ st.to_string(unit=u.hour))
|
|
print("\033[36m"
|
|
+ "\t end sidereal time: "
|
|
+ "\033[0m"
|
|
+ (st + obs_sky_rotation).wrap_at(24*u.h).to_string(unit=u.hour))
|
|
print("\033[36m"
|
|
+ "\t sky rotation: "
|
|
+ "\033[0m"
|
|
+ obs_sky_rotation.to_string(unit=u.hour))
|
|
print("\033[36m"
|
|
+ "Sky coordinates window"
|
|
+ "\033[0m")
|
|
print("\033[36m"
|
|
+ "\t east ra: "
|
|
+ "\033[0m"
|
|
+ a_east.to_string(unit=u.hour))
|
|
print("\033[36m"
|
|
+ "\t west ra: "
|
|
+ "\033[0m"
|
|
+ a_west.to_string(unit=u.hour))
|
|
print("\033[36m"
|
|
+ "\tupper dec: "
|
|
+ "\033[0m"
|
|
+ dec_upper.to_string(unit=u.degree))
|
|
print("\033[36m"
|
|
+ "\tlower dec: "
|
|
+ "\033[0m"
|
|
+ dec_lower.to_string(unit=u.degree))
|
|
print("\033[36m"
|
|
+ "Simbad query constraints: "
|
|
+ "\033[0m"
|
|
+ "WHERE (ra < {} {} ra > {}) AND (dec < {} AND dec > {}))".format(
|
|
a_west.degree,
|
|
comp_text,
|
|
a_east.degree,
|
|
dec_upper.degree,
|
|
dec_lower.degree))
|
|
print("\033[36m"
|
|
+ "Vizier query constraints:"
|
|
+ "\033[0m")
|
|
print("\033[36m"
|
|
+ "\tRA: "
|
|
+ "\033[0m"
|
|
+ "< "
|
|
+ a_west.to_string(unit=u.hour, sep=":")
|
|
+ " {} ".format(comp_symb)
|
|
+ "> "
|
|
+ a_east.to_string(unit=u.hour, sep=":"))
|
|
print("\033[36m"
|
|
+ "\tDE: "
|
|
+ "\033[0m"
|
|
+ "< "
|
|
+ dec_upper.to_string(unit=u.degree, sep=":")
|
|
+ " && "
|
|
+ "> "
|
|
+ dec_lower.to_string(unit=u.degree, sep=":"))
|
|
|
|
print("")
|
|
print("\033[32m"
|
|
+ "Write output to file ? [yes/no]"
|
|
+ "\033[0m")
|
|
write_output = input("\033[32m"
|
|
+ "Choice: "
|
|
+ "\033[0m").upper() in YES
|
|
if write_output:
|
|
try:
|
|
with open("{}{}_{}.cfg".format(OUT_DIR, obs_date, obs_short),
|
|
"w+") as file:
|
|
file.write("# File generated with SCOPE_v2\n")
|
|
file.write("LOCATION: ")
|
|
file.write(obs_short + "\n")
|
|
file.write("LAT: ")
|
|
file.write(obs_lat.to_string() + "\n")
|
|
file.write("LON: ")
|
|
file.write(obs_lon.to_string() + "\n")
|
|
file.write("SUN_SET: ")
|
|
file.write(sun_set.to_string() + "\n")
|
|
file.write("SUN_SET_CIVIL: ")
|
|
file.write(sun_civil.to_string() + "\n")
|
|
file.write("SUN_SET_NAUTICAL: ")
|
|
file.write(sun_nautical.to_string() + "\n")
|
|
file.write("SUN_SET_ASTRONOMICAL: ")
|
|
file.write(sun_astronomical.to_string() + "\n")
|
|
file.write("SUN_RISE_ASTRONOMICAL: ")
|
|
file.write(sun_rastronomical.to_string() + "\n")
|
|
file.write("SUN_RISE_NAUTICAL: ")
|
|
file.write(sun_rnautical.to_string() + "\n")
|
|
file.write("SUN_RISE_CIVIL: ")
|
|
file.write(sun_rcivil.to_string() + "\n")
|
|
file.write("SUN_RISE: ")
|
|
file.write(sun_rise.to_string() + "\n")
|
|
file.write("OBS_BEGIN: ")
|
|
file.write(obs_begin_time.to_string() + "\n")
|
|
file.write("OBS_END: ")
|
|
file.write(obs_end_time.to_string() + "\n")
|
|
file.write("ST_BEGIN: ")
|
|
file.write(st.to_string(unit=u.hour) + "\n")
|
|
file.write("ST_END: ")
|
|
file.write((st + obs_sky_rotation).wrap_at(
|
|
24*u.h).to_string(unit=u.hour) + "\n")
|
|
file.write("WINDOW_EAST: ")
|
|
file.write(a_east.to_string(unit=u.hour) + "\n")
|
|
file.write("WINDOW_WEST: ")
|
|
file.write(a_west.to_string(unit=u.hour) + "\n")
|
|
file.write("WINDOW_UPPER: ")
|
|
file.write(dec_upper.to_string(unit=u.degree) + "\n")
|
|
file.write("WINDOW_LOWER: ")
|
|
file.write(dec_lower.to_string(unit=u.degree) + "\n")
|
|
print("\033[34m"
|
|
+ "File saved in the {} directory ".format(OUT_DIR)
|
|
+ "with the name {}_{}.cfg".format(obs_date, obs_short)
|
|
+ "\033[0m")
|
|
print("\033[34m" + "done!" + "\033[0m")
|
|
except:
|
|
print("\033[91m"
|
|
+ "Error! File cannot be save in the output directory."
|
|
+ "\033[0m")
|
|
return 0
|
|
if __name__ == "__main__":
|
|
main()
|