<jupyter_start><jupyter_text>False Star Tests<jupyter_code>x_pos = 3150 + 120*np.random.rand(50) #2524 + 120*np.random.rand(50)
y_pos = 3380 + 90*np.random.rand(50) #2530 + 120*np.random.rand(50)

f435w_mags = 17 + 12*np.random.rand(1000) # Range of 17 to 29 in mag
Iax_dir = "/Users/tktakaro/Documents/Type-Iax-HST"
file = open(Iax_dir + "/new_files/sn2014dt/fake_stars_f438w.list", "w+")
for j in range(x_pos.size):
    for i in range(f435w_mags.size):
        file.write('0 1 {x} {y} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4} {f4}\n'.format(x=x_pos[j], y=y_pos[j], f4=f435w_mags[i]))
file.close()

x_pos = 3150 + 120*np.random.rand(50)
y_pos = 3380 + 90*np.random.rand(50)

f435w_mags = 17 + 12*np.random.rand(1000) # Range of 17 to 29 in mag
Iax_dir = "/Users/tktakaro/Documents/Type-Iax-HST"
file = open(Iax_dir + "/new_files/sn2014dt/fake_stars_f275w.list", "w+")
for j in range(x_pos.size):
    for i in range(f435w_mags.size):
        file.write('0 1 {x} {y} {f4} {f4} {f4} {f4}\n'.format(x=x_pos[j], y=y_pos[j], f4=f435w_mags[i]))
file.close()
# Loading in the false star test
catalog_false_f275w = np.loadtxt(Iax_dir + "/new_files/sn2014dt/sn2014dt_f275w.fake")
catalog_false_f438w = np.loadtxt(Iax_dir + "/new_files/sn2014dt/sn2014dt_f438w.fake")
recovered_mags = np.zeros([50000,14]) # 1. X position, 2. Y position, 3. Inserted Mag, 4. Recovered mag (f438w),
# 5. S/N (f438w), 6. Sharpness (f438w), 7. Roundness (f438w), 8. Crowding (f438w), 9.+ same, for other filters
for i in range(catalog_false_f438w.shape[0]):
    # X, Y, and Inserted Magnitude
    recovered_mags[i,0] = catalog_false_f438w[i,2]
    recovered_mags[i,1] = catalog_false_f438w[i,3]
    recovered_mags[i,2] = catalog_false_f438w[i,5]
    # Recovered mag in F275W
    if (catalog_false_f275w[i,14] == 0) and (catalog_false_f275w[i,15] == 0):
        recovered_mags[i,3] = 99.999
    else:
        recovered_mags[i,3] = catalog_false_f275w[i,27]
        recovered_mags[i,4:8] = catalog_false_f275w[i,17:21]
    # Recovered mag in F438W
    if (catalog_false_f438w[i,14] == 0) and (catalog_false_f438w[i,15] == 0):
        recovered_mags[i,8] = 99.999
    else:
        recovered_mags[i,8] = catalog_false_f438w[i,59]
        recovered_mags[i,9:13] = catalog_false_f438w[i,49:53]
plt.plot(recovered_mags[:,2], recovered_mags[:,3], '.', alpha=.05)
plt.plot(recovered_mags[:,2], recovered_mags[:,2], '--')
plt.ylim(29,22)
plt.xlim(29,22)
plt.xlabel("Inserted Magnitude")
plt.ylabel("Recovered Mag (f275w)")
plt.show()

plt.plot(recovered_mags[:,2], recovered_mags[:,8], '.', alpha=.05)
plt.plot(recovered_mags[:,2], recovered_mags[:,2], '--')
plt.ylim(29,22)
plt.xlim(29,22)
plt.xlabel("Inserted Magnitude")
plt.ylabel("Recovered Mag (f438w)")
plt.show()
# Make the same cuts we do for the real data
recovered_cut = np.zeros(recovered_mags.shape)
j = 0
for i in range(catalog_false_f438w.shape[0]):
    if ((recovered_mags[i,4] >= 2.5) and (abs(recovered_mags[i,5]) < .3) and (abs(recovered_mags[i,6]) < 1)
        and (recovered_mags[i,7] < .1)):
        recovered_cut[j,:8] = recovered_mags[i,:8]
    if ((recovered_mags[i,9] >= 2.5) and (abs(recovered_mags[i,10]) < .3) and (abs(recovered_mags[i,11]) < 1)
        and (recovered_mags[i,12] < .1)):
        recovered_cut[j,8:13] = recovered_mags[i,8:13]
    j += 1
recovered_cut = recovered_cut[~np.all(recovered_cut == 0, axis=1)]
plt.plot(recovered_cut[:,2], recovered_cut[:,3], '.', alpha=.5)
plt.plot(recovered_cut[:,2], recovered_cut[:,2], '--')
plt.ylim(29,22)
plt.xlim(29,22)
plt.xlabel("Inserted Magnitude")
plt.ylabel("Recovered Mag (f275w)")
plt.show()

plt.plot(recovered_cut[:,2], recovered_cut[:,8], '.', alpha=.5)
plt.plot(recovered_cut[:,2], recovered_cut[:,2], '--')
plt.ylim(29,22)
plt.xlim(29,22)
plt.xlabel("Inserted Magnitude")
plt.ylabel("Recovered Mag (f438w)")
plt.show()
f275w_mags = []
f275w_recover = []
f438w_mags = []
f438w_recover = []
# Calculate recovery percentage as a function of inserted mag
for i in range(catalog_false_f438w.shape[0]):
    f275w_mags.append(recovered_mags[i,2])
    f438w_mags.append(recovered_mags[i,2])
    if abs(recovered_mags[i,2] - recovered_mags[i,3]) < .5:
        f275w_recover.append(recovered_mags[i,2])
    if abs(recovered_mags[i,2] - recovered_mags[i,8]) < .5:
        f438w_recover.append(recovered_mags[i,2])
# Setting up the histogram bins w/ .1 mag resolution
b = []
for i in range(121):
    b.append(17 + i/10.)
plt.subplots(figsize=(15, 10))
#b = [17,17.5,18,18.5,19,19.5,20,20.5,21,21.5,22,22.5,23,23.5,24,24.5,25,25.5,26,26.5,27,27.5,28,28.5,29,29.5,30]
h2 = plt.hist(f275w_mags, bins=b, facecolor='navy', rwidth=.95, label='Inserted stars')
r2 = plt.hist(f275w_recover, bins=b, rwidth=.95, facecolor='khaki', label='Recovered stars')
plt.plot(h2[1][:-1]+.05, .9*h2[0], 'r_', markersize=6, label='90% inserted stars')
plt.xlabel("Inserted Magnitude (f275w)", fontsize=24)
plt.ylabel("Number", fontsize=24)
plt.legend(fontsize=20)
plt.show()

plt.subplots(figsize=(15, 10))
h4 = plt.hist(f438w_mags, bins=b, facecolor='navy', rwidth=.95, label='Inserted stars')
r4 = plt.hist(f438w_recover, bins=b, rwidth=.95, facecolor='khaki', label='Recovered stars')
plt.plot(h4[1][:-1]+.05, .9*h4[0], 'r_', markersize=6, label='90% inserted stars')
plt.xlabel("Inserted Magnitude (f438w)", fontsize=24)
plt.ylabel("Number", fontsize=24)
plt.legend(fontsize=20)
plt.show()

frac_f275w = r2[0]/h2[0]
frac_f438w = r4[0]/h4[0]
from scipy.interpolate import interp1d

plt.plot(h2[1][:-1]+.05, frac_f275w, 'm--')
plt.plot(h4[1][:-1]+.05, frac_f438w, 'b--')
plt.plot(h4[1][:-1]+.05, .9*np.ones(h4[1][:-1].shape), 'k-')

limit_mags = np.zeros(2)
limit_mags[0] = h2[1][:-1][np.argmin(abs(frac_f275w - .9))] + .05
limit_mags[1] = h4[1][60:-1][np.argmin(abs(frac_f438w[60:] - .9))] + .05
print(limit_mags)
kwargs = {"names": ["EEP", "log10_isochrone_age_yr", "initial_mass", "star_mass", "log_Teff", "log_g", "log_L",
                        "[Fe/H]_init", "[Fe/H]", "WFC3_UVIS_F200LP", "WFC3_UVIS_F218W", "WFC3_UVIS_F225W",
                        "WFC3_UVIS_F275W", "WFC3_UVIS_F280N", "WFC3_UVIS_F300X", "WFC3_UVIS_F336W",
                        "WFC3_UVIS_F343N", "WFC3_UVIS_F350LP", "WFC3_UVIS_F373N", "WFC3_UVIS_F390M",
                        "WFC3_UVIS_F390W", "WFC3_UVIS_F395N", "WFC3_UVIS_F410M", "WFC3_UVIS_F438W",
                        "WFC3_UVIS_F467M", "WFC3_UVIS_F469N", "WFC3_UVIS_F475W", "WFC3_UVIS_F475X",
                        "WFC3_UVIS_F487N", "WFC3_UVIS_F502N", "WFC3_UVIS_F547M", "WFC3_UVIS_F555W",
                        "WFC3_UVIS_F600LP", "WFC3_UVIS_F606W", "WFC3_UVIS_F621M", "WFC3_UVIS_F625W",
                        "WFC3_UVIS_F631N", "WFC3_UVIS_F645N", "WFC3_UVIS_F656N", "WFC3_UVIS_F657N",
                        "WFC3_UVIS_F658N", "WFC3_UVIS_F665N", "WFC3_UVIS_F673N", "WFC3_UVIS_F680N",
                        "WFC3_UVIS_F689M", "WFC3_UVIS_F763M", "WFC3_UVIS_F775W", "WFC3_UVIS_F814W",
                        "WFC3_UVIS_F845M", "WFC3_UVIS_F850LP", "WFC3_UVIS_F953N", "WFC3_IR_F098M",
                        "WFC3_IR_F105W", "WFC3_IR_F110W", "WFC3_IR_F125W", "WFC3_IR_F126N", "WFC3_IR_F127M",
                        "WFC3_IR_F128N", "WFC3_IR_F130N", "WFC3_IR_F132N", "WFC3_IR_F139M", "WFC3_IR_F140W",
                        "WFC3_IR_F153M", "WFC3_IR_F160W", "WFC3_IR_F164N", "WFC3_IR_F167N", "phase"],
            "delim_whitespace": True, "comment": "#"}
isochrones = {}
for filename in glob.glob(mist_dir + "/*.iso.cmd"):
    filename = filename.replace("\\", "/")
    feh_string = filename.split("/")[-1].split("=")[1].split(".iso")[0] # Pulls metalicity information
    feh = float(feh_string) # feh is [Fe/H]
    df = pandas.read_csv(filename, **kwargs)
    df['ages'] = 10 ** df.log10_isochrone_age_yr / 1e9
    isochrones[feh] = df # Creates dictionary accessible by entering a particular metalicity
df = isochrones[0.16]
ages = np.array(list(set(df.log10_isochrone_age_yr)))
ages.sort()

age_cmd = {}

ages = ages[(ages > 6.49) & (ages < 7.26)]
for i, age in enumerate(ages):
    idx = df.log10_isochrone_age_yr == age
    # Three adjustments: 1. Distance: added to mag 2. MW Extinction: added to both. 3. Galactic extinction
    dist_adjust = 5 * (np.log10(19.3) - 1)
    F275W_ext = 0.123 # extinction in F275W in M61 from NED
    F438W_ext = 0.081 # extinction in F555W in M61 from NED
    Gal_ext = 0
    # Multiply Gal_ext by A_lambda/E(B-V) values from Schlafly, Finkbeiner 2011
    mag = df[idx].WFC3_UVIS_F275W + dist_adjust + F275W_ext + 5.487*Gal_ext
    color = df[idx].WFC3_UVIS_F275W - df[idx].WFC3_UVIS_F438W  + F275W_ext - F438W_ext + (5.487 - 3.623) * Gal_ext
    mass = df[idx].initial_mass
    age_cmd[age] = (color, mag, mass)

norm = mpl.colors.Normalize(vmin=ages.min(), vmax=ages.max())
sm = mpl.cm.ScalarMappable(norm=norm, cmap='viridis')
sm.set_array(ages)
plt.subplots(figsize=(15, 10))

x_array = np.linspace(-1,2,100)
plt.fill_between(x_array, 99*np.ones(100), 24 - x_array, color='lightgrey', zorder=0)
for age, (color, mag, mass) in age_cmd.items():
    plt.plot(color, mag, alpha=0.5, color=sm.to_rgba(age), zorder=1, linewidth=4)
x_array = np.linspace(-1, 2, 100)
    
plt.xlabel("F275W - F438W")
plt.ylabel("F438W")
plt.ylim(28,22)
plt.xlim(-1, 2)
cb = plt.colorbar(sm)
cb.set_label('$\log_{10}$ Age [yr]')
plt.draw()
plt.show()
import scipy.integrate as integrate
def kroupa(m):
    norm = 0.32705769605379353 # integrated kroupa from .075 to 100
    if m < .08:
        alpha = .3
        return m**(-alpha)/norm
    elif m < .5:
        alpha = 1.3
        return .08 * m**(-alpha)/norm
    else:
        alpha = 2.3
        return .5 * .08 * m**(-alpha)/norm
def invertedKroupa(u):
    if u > kroupa(.08):
        return - np.log(u)/.3
    elif u > kroupa(.5):
        return - np.log(u)/(1.3 * .08)
    else:
        return - np.log(u)/(2.3 * .04)
    
print(integrate.quad(kroupa, 0.075, 100))
ages = np.array(list(set(df.log10_isochrone_age_yr)))
ages = ages[(ages > 6.49) & (ages < 8.51)]
ages.sort()
age_cmd = {}
dist_adjust = 5 * (np.log10(19.3e6) - 1)

ages_ruledout = []
for i, age in enumerate(ages):
    idx = df.log10_isochrone_age_yr == age
    Gal_ext = 0 
    mag2 = df[idx].WFC3_UVIS_F275W + dist_adjust + F275W_ext + 5.487*Gal_ext
    mag4 = df[idx].WFC3_UVIS_F438W + dist_adjust + F438W_ext + 3.623*Gal_ext
    mass = df[idx].initial_mass
    age_cmd[age] = (mag2, mag4, mass)

for age, (mag2, mag4, mass) in age_cmd.items():
    indices2 = []
    for i, mag in enumerate(mag2):
        if mag < limit_mags[0]:
            indices2.append(i)
    indices4 = []
    for i, mag in enumerate(mag4):
        if mag < limit_mags[1]:
            indices4.append(i)
            
    masses2 = []
    masses4 = []
    for i, mass in enumerate(mass):
        if np.isin(i, indices2):
            masses2.append(mass)
        if np.isin(i, indices4):
            masses4.append(mass)
            
    total_mass = 0
    n_inrange = 0
    while total_mass < 50:
        a = invertedKroupa(np.random.random())
        while a < .075 or a > 100: # Only uses stars with mass between .075 and 100 solar masses
            a = invertedKroupa(np.random.random())
        total_mass += a
        if not masses2 == []:
            if a > np.min(masses2) and a < np.max(masses2):
                n_inrange += 1
                continue
        if not masses4 == []:
            if a > np.min(masses4) and a < np.max(masses4):
                n_inrange += 1
                continue

    if n_inrange > 0:
        ages_ruledout.append(age)
print(masses4)
ages_ruledout.sort()
print(np.round(ages_ruledout,2))<jupyter_output>[]
{'fraction_non_alphanumeric': 0.17415199656504937, 'fraction_numerical': 0.09901245169600686, 'mean_word_length': 3.6068037974683542, 'pattern_counts': {'":': 3, '<': 22, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '570003', 'license': 'MIT', 'n_tokens_mistral': 5794, 'n_tokens_neox': 5086, 'n_words': 878, 'path': 'jupyter_notebooks/2014dt.ipynb', 'repo_name': 'TTakaro/Type-Iax-HST', 'size': '11645'}