from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation
from astropy.coordinates import FK4
from astropy.constants import c
import astropy.units as u
import numpy as np
import sys

if len(sys.argv) < 2:
    print("usage: python " + sys.argv[0] + " YYYY-MM-DD")
    sys.exit()

date_str = sys.argv[1]

# ----------------------------
# Date (UTC)
# ----------------------------
t0 = Time(date_str+"T00:00:00", scale="utc")
t0_tdb = Time(date_str+"T00:00:00", scale="tdb")

# time duration : 36 hours
dt = np.arange(0, 36*3600) * u.s
# time array
t_array = t0 + dt

# Day of year
doy = t0.datetime.timetuple().tm_yday

year_str = f"{t0.datetime.year:04d}" 
doy_str  = f"{doy:03d}"

# ----------------------------
# Observer location (IPRT)
# Geographycal latitude, longitude, altitude
# ----------------------------
location = EarthLocation(
    lat=37.68 * u.deg,
    lon=140.74 * u.deg,
    height=600.0 * u.m
)
RX_name = "IttP" 

# ----------------------------
# Target source position (RADEC) (3C144）
# ----------------------------
target = SkyCoord(
    ra=(5.0 + 31.0/60.0 + 31.4/3600.0) / 24.0 * 360.0 * u.deg,
    dec=(21.0 + 58.0/60.0 + 55.0/3600.0) * u.deg,
    frame=FK4(equinox=Time("B1950"))
)

# ----------------------------
# Light time correction (barycentric)
# ----------------------------
ltt_bary = t_array.light_travel_time(target, location=location)

# ----------------------------
# BJD_TDB
# ----------------------------
t_tdb = t_array.tdb
bjd_tdb = t_tdb + ltt_bary

# ----------------------------
# line-of-sight velocity (barycentric correction)
# ----------------------------
vcorr = target.radial_velocity_correction(
    kind="barycentric",
    obstime=t_array,
    location=location
)

# ----------------------------
# Doppler factor
# ----------------------------
f2f = ((1 + vcorr/c) / (1 - vcorr/c))**0.5

# ----------------------------
# Output data
# ----------------------------
dir_save = "./"
file_tdb = dir_save + year_str + doy_str + "_" + RX_name + "_TDB"

# ascii
np.savetxt(file_tdb+".txt", (bjd_tdb-t0_tdb).sec, fmt="%.9f")
np.savetxt(file_tdb+"_fact2Fism.txt", f2f.value, fmt="%.9f")

# IEEE 754 double（Binary64）
arr = np.array((bjd_tdb-t0_tdb).sec, dtype=np.float64)
arr.tofile(file_tdb+".bin")
arr = np.array(f2f.value, dtype=np.float64)
arr.tofile(file_tdb+"_fact2Fism.bin")