#!/usr/bin/env python
"""

28 June 2023 (das), read the FITS files that Paul made from Lucas' ascii files.
25 Julye 2023 (psr), Cleaned up code a bit

"""
import sys, os
import astropy.io.fits as fits
from matplotlib import pyplot as plt

if len(sys.argv) < 2:
    print("Usage: Example_3PC_PlotProfileData.py <JHHMM[+/-]DDMM>")
    sys.exit(1)

thisFile = "3PC_ProfileData/%s_3PC_data.fits" % sys.argv[1]
thisPlot = "%s_3PC_data.png" % sys.argv[1]

hasFit = False
hasRadio = False
with fits.open(thisFile) as f:
    dataphmin = f["GAMMA_LC"].data["Ph_Min"]
    dataphmax = f["GAMMA_LC"].data["Ph_Max"]
    dataCounts = f["GAMMA_LC"].data["GT100_WtCnt"]
    dataCountsUnc = f["GAMMA_LC"].data["Unc_GT100_WtCnt"]

    if "BEST_FIT_LC" in f:
        fitphmin = f["BEST_FIT_LC"].data["Ph_Min"]
        fitphmax = f["BEST_FIT_LC"].data["Ph_Max"]
        fit = f["BEST_FIT_LC"].data["Intensity"]
        hasFit = True

    if "RADIO_PROFILE" in f:
        radiophmin = f["RADIO_PROFILE"].data["Ph_Min"]
        radiophmax = f["RADIO_PROFILE"].data["Ph_Max"]
        radio = f["RADIO_PROFILE"].data["Norm_Intensity"]
        hasRadio = True

myplots = plt.figure(figsize=(10, 10))
ax1 = myplots.add_subplot(1, 1, 1)

middle = (dataphmin + dataphmax) / 2.0
ax1.plot(middle, dataCounts, "k-", label="LAT Data")

if hasRadio:
    middle = (radiophmin + radiophmax) / 2.0
    ax1.plot(middle, radio, "r-", label="Radio")

if hasFit:
    middle = (fitphmin + fitphmax) / 2.0
    ax1.plot(middle, fit, "b-", label="Profile Fit")

ax1.set_xlabel("Phase")
ax1.set_ylabel("Weighted counts")
ax1.set_title(f"3PC Profile for {sys.argv[1]}")
ax1.grid(True)
ax1.legend()
ax1.set_xlim((0.0, 2.0))
myplots.savefig(thisPlot)
plt.show()
plt.close(myplots)