-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_calfacs_wave.py
117 lines (96 loc) · 3.79 KB
/
plot_calfacs_wave.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
import argparse
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import QTable, join
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--png", help="save figure as a png file", action="store_true")
parser.add_argument("--pdf", help="save figure as a pdf file", action="store_true")
args = parser.parse_args()
# read in the calfactors for each filter
waves = np.array([5.6, 7.7, 10.0, 11.3, 12.8, 15.0, 18.0, 21.0, 25.5])
filters = [
"F560W",
"F770W",
"F1000W",
"F1130W",
"F1280W",
"F1500W",
"F1800W",
"F2100W",
"F2550W",
]
atab = None
for cfilter in filters:
ctab = QTable.read(f"CalFacs/miri_calfactors_{cfilter}.fits")
# merge multiple measurements of the same star
gvals = np.full(len(ctab), False)
for cname in np.unique(ctab["name"]):
sindxs, = np.where(ctab["name"] == cname)
if len(sindxs) > 1:
ctab[f"calfac_{cfilter}"][sindxs[0]] = np.average(ctab[f"calfac_{cfilter}"][sindxs])
ctab[f"calfac_{cfilter}_mean_dev"][sindxs[0]] = np.average(ctab[f"calfac_{cfilter}_mean_dev"][sindxs])
ctab[f"calfac_{cfilter}_med_dev"][sindxs[0]] = np.average(ctab[f"calfac_{cfilter}_med_dev"][sindxs])
gvals[sindxs[0]] = True
ctab = ctab[gvals]
if atab is not None:
atab = join(atab, ctab, join_type="outer")
else:
atab = ctab
nstars = len(atab)
atab.write("miri_calfactors_all.fits", overwrite=True)
# these two do not have F1500W obs, so nothing in the model column
mindx, = np.where(atab["name"] == "GD 71")
atab["modflux_F1500W"][mindx[0]] = 0.0
mindx, = np.where(atab["name"] == "2MASS J17571324+6703409")
atab["modflux_F1500W"][mindx[0]] = 0.0
# sort by F1150W model flux as all are measured in this band
aindxs = np.argsort(atab["modflux_F1500W"])
atab = atab[aindxs]
# print(atab)
# atab.write("test.dat", format="ipac", overwrite=True)
# exit()
# make plot
fontsize = 14
font = {"size": fontsize}
plt.rc("font", **font)
plt.rc("lines", linewidth=2)
plt.rc("axes", linewidth=2)
plt.rc("xtick.major", width=2)
plt.rc("ytick.major", width=2)
ignore_names = ["HD 167060", "16 Cyg B", "HD 37962", "del UMi", "HD 180609"]
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(14, 8))
colormap = plt.cm.gist_ncar
plt.gca().set_prop_cycle(plt.cycler('color', plt.cm.jet(np.linspace(0, 1, nstars))))
lstyle = ["solid", "dotted", "dashed", "dashdot"]
for k, cname in enumerate(atab["name"]):
pfacs = []
pwaves = []
for j, cfilter in enumerate(filters):
pwaves.append(waves[j])
pfacs.append(atab[k][f"calfac_{cfilter}_med_dev"])
pwaves = np.array(pwaves)
pfacs = np.array(pfacs)
if np.sum(np.isfinite(pfacs)) > 1:
gvals = np.isfinite(pfacs)
ax.plot(pwaves[gvals], pfacs[gvals], linestyle=lstyle[k % 4], label=cname)
else:
ax.plot(pwaves, pfacs, "ko", label=cname)
#if cname in ignore_names:
# determine the calfac/mean calfac ratio for bands <= F1500W
cvals = (pwaves < 16.) & np.isfinite(pfacs)
print(cname, np.average(pfacs[cvals]), np.std(pfacs[cvals]))
#print(pwaves[cvals])
#print(pfacs[cvals])
# ax.set_xlim(5.0, 35.0)
ax.set_xlabel(r"$\lambda$ [$\mu$m]")
ax.set_ylabel("calfac / (ave calfac)")
ax.legend(fontsize=0.8 * fontsize, ncol=2)
plt.tight_layout()
fname = "Figs/miri_calfactors_allwaves"
if args.png:
fig.savefig(f"{fname}.png")
elif args.pdf:
fig.savefig(f"{fname}.pdf")
else:
plt.show()