This Python script is designed to analyze synthetic stellar populations and their fiducial lines in the context of astronomy. It involves generating synthetic populations based on specific parameters, such as metallicity, age, and binary fraction, and then measuring the goodness-of-fit between observed and synthetic fiducial lines using a Chi-squared method. The script imports modules for data manipulation (pandas), mathematical operations (NumPy, SciPy), and file management, indicating a comprehensive approach to processing astronomical data. The code includes functionalities for loading pre-computed empirical fiducials, handling bolometric correction tables, interpolating synthetic fiducial lines, and performing optimizations to minimize the difference between observed and modeled data. This script is likely part of a larger research project aimed at understanding the properties of globular clusters or stellar populations within galaxies.

from fidanka.population.synthesize import population
from fidanka.fiducial.fiducial import fiducial_line
from fidanka.isofit.fit import pfD
 
import pickle as pkl
import os
import pathlib
import pandas as pd
from scipy.interpolate import interp1d
from scipy.optimize import minimize
from itertools import permutations, product
import re
import numpy as np
 
from tqdm import tqdm
 
ISOPATH = "../../../outputs.denseAlpha/"
ARTSTARTEST = "/home/tboudreaux/d/Astronomy/GraduateSchool/Thesis/GCConsistency/NGC2808/ArtificialStarTests/ArtificialStarCalibratedErrorFunctions.pkl"
RESPONSEFUNCTIONS = "../../../inputs/responseFunctions/"
BOLTABLES = "../../../bolTables/HSTWFC3/"
 
 
def obs_synth_chi2(obs, synth):
    oMag, oColor, oWidth = obs[:,0], obs[:,1], obs[:,2]
    sMag, sColor, sWidth = synth[:,0], synth[:,1], synth[:,2]
 
    # Interpolate the synthetic fiducial to the observed magnitudes
    sColorFunc = interp1d(sMag, sColor, kind='linear', fill_value='extrapolate', bounds_error=False)
    sWidthFunc = interp1d(sMag, sWidth, kind='linear', fill_value='extrapolate', bounds_error=False)
 
    sColor = sColorFunc(oMag)
    sWidth = sWidthFunc(oMag)
 
    oColorFunc = interp1d(oMag, oColor, kind='linear', fill_value='extrapolate', bounds_error=False)
 
    colorChi2 = np.empty(len(oMag))
    colorChi2[:] = np.nan
    widthChi2 = np.empty(len(oMag))
 
    for idx, (mag, sC, sW, oW) in enumerate(zip(oMag, sColor, sWidth, oWidth)):
        d = pfD((sC,mag), oColorFunc)
        nearest = minimize(d, mag, method='Nelder-Mead')
        if not nearest.success:
            print("Warning: minimization failed")
        else:
            colorChi2[idx] = d(nearest.x)**2
        widthChi2[idx] = ((oW - sW)**2/sW)
 
    colorChi2 = np.sum(colorChi2)
    widthChi2 = np.sum(widthChi2)
    return (colorChi2 + widthChi2)/(2*len(oMag))
 
def build_and_measure_population(mu, E, fb, artStar, targetAge, iso, alpha, bolTables, n=10000):
    pop = population(
            iso,
            alpha,
            fb, lambda x: x-x + targetAge,
            n,
            targetAge,
            targetAge,
            0.25,
            2,
            artStar,
            10**((mu+5)/5),
            E,
            "F606W",
            bolTables
            )
 
    df = pop.to_pandas()
    fid = fiducial_line(
            df["F606W"],
            df["F814W"],
            df["F606W_err"],
            df["F814W_err"],
            mcruns=1,
            percHigh=97,
            percLow=4,
            method='spline'
            )
 
    # only retain mag bin, color loc, and width
    return df, fid[:,(0,1,4)]
 
 
if __name__ == "__main__":
    assert 'HUGSFiducials.pkl' in os.listdir('.'), "Cannot find pre computed empirical fiducials"
 
    with open(ARTSTARTEST, 'rb') as f:
        artStar = pkl.load(f)
 
    with open('./HUGSFiducials.pkl', 'rb') as f:
        hugFids = pkl.load(f)
 
    bolTables = list(filter(lambda x: re.match(r'feh[mp]\d{3}\.HST_WFC3', x), os.listdir(BOLTABLES)))
    bolTablePaths = [os.path.join(BOLTABLES, x) for x in bolTables]
 
    isos = list()
    isoGlob = pathlib.Path(ISOPATH).rglob("isochrones.txt")
    for idx, iso in tqdm(enumerate(isoGlob), desc="Scanning for bolometrically corrected isochrones"):
        isos.append(str(iso))
    AIsos = list(filter(lambda x: 'PopA' in x, isos))
    EIsos = list(filter(lambda x: 'PopE' in x, isos))
 
    uCombs = list(product(AIsos, EIsos))
    uCombsValid = filter(lambda x: float(re.search('Pop[AE]\+(\d+\.\d+)', x[0]).groups()[0]) < float(re.search('Pop[AE]\+(\d+\.\d+)', x[1]).groups()[0]), uCombs)
    uCombsValid = filter(lambda x: re.search('alpha-(\d+\.\d+)', x[0]).groups()[0] == re.search('alpha-(\d+\.\d+)', x[1]).groups()[0], uCombsValid)
 
    mu = 14.89
    E = 0.17
    fb = 0.02
    targetAge = 12.5e9
    for idx, isoComb in enumerate(uCombsValid):
        # if idx < 200:
        #     continue
        print(f"Working on {idx} with {isoComb[0]}")
        print(artStar.keys())
        moveOn = False
        while not moveOn:
            synthPopFiducial = build_and_measure_population(mu, E, fb, artStar.copy(), targetAge, list(isoComb), -0.84, bolTablePaths)
            if len(synthPopFiducial[1][:,0]) == 0:
                print("Warning: synthetic fiducial failure, retrying")
            else:
                moveOn = True
        chi2 = obs_synth_chi2(hugFids, synthPopFiducial[1])
        print(chi2)
        with open(f"DumpTestResults_{idx}.pkl", "wb") as f:
            pkl.dump(synthPopFiducial, f)
 
        if idx == 0:
            exit()