<commit_before>import csv
import os
import sys
import numpy as np
from fitelp.label_tools import line_label
from fitelp.read_spectra import GalaxyRegion
from fitelp.fit_line_profiles import FittingProfile, plot_profiles, vel_to_wave, wave_to_vel
from fitelp.make_latex_tables import average_velocities_table_to_latex, halpha_regions_table_to_latex, comp_table_to_latex
from fitelp.bpt_plotting import calc_bpt_points, bpt_plot
import fitelp.constants as constants

constants.init()
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '../scripts'))
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '../Input_Galaxy_Region_Information'))
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '../Input_Data_Files'))


def calc_vel_dispersion(sigmaObs, sigmaObsError, sigmaTemp2, filter, rp, correctCopiedSigma=False):
    # Assuming negligible error in temp or instrument
    if filter == 'blue':
        sigmaInstr = rp.sigmaInstrBlue
    elif filter == 'red':
        sigmaInstr = rp.sigmaInstrRed

    #Old version (added by Vero)
    # totalSigmaSquared = sigmaObs ** 2 - sigmaInstr ** 2 - sigmaTemp2
    # totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
    # try:
    #     intrinsic = np.sqrt(totalSigmaSquared)
    #     intrinsicError = 0.5 * totalSigmaSquared ** (-0.5) * totalSigmaSquaredError
    # except ValueError:
    #     "ERROR: INVALID SIGMA"
    #     intrinsic = 0
    #     intrinsicError = 0
    #
    # Last version (added by Vero)
    #if sigmaObs**2 > sigmaTemp2:
        # totalSigmaSquared = sigmaObs ** 2 - sigmaInstr ** 2 - sigmaTemp2
        # totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
    # else:
    #     totalSigmaSquared = (sigmaInstr + np.sqrt(sigmaTemp2))**2
    #     totalSigmaSquaredError = 0
    # intrinsic = np.sqrt(totalSigmaSquared)
    # intrinsicError = 0.5 * totalSigmaSquared ** (-0.5) * totalSigmaSquaredError
    #
    if correctCopiedSigma:
        # SigmaObs here is the SigmaIntrinsic of the copied line
        totalSigmaSquared = sigmaObs**2 + sigmaInstr**2 + sigmaTemp2
        totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
    else:
        if sigmaObs**2 > (sigmaTemp2 + sigmaInstr**2):
            totalSigmaSquared = sigmaObs**2 - sigmaInstr**2 - sigmaTemp2
            totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
        else:
            print("ERROR: INVALID SIGMA")
            totalSigmaSquared = 1e-99
            totalSigmaSquaredError = 0
    intrinsic = np.sqrt(totalSigmaSquared)
    intrinsicError = 0.5 * totalSigmaSquared**(-0.5) * totalSigmaSquaredError

    return intrinsic, intrinsicError


def calc_flux(height, sigmaObs, heightError, sigmaObsError):
    insideSqrt = 2 * np.pi * sigmaObs**2
    insideSqrtError = 2 * sigmaObs * sigmaObsError
    sqrt = np.sqrt(insideSqrt)
    sqrtError = 0.5 * insideSqrt**(-0.5) * insideSqrtError

    calcFlux = height * sqrt
    calcFluxError = calcFlux * np.sqrt((heightError/height)**2 + (sqrtError/sqrt)**2)

    return calcFlux, calcFluxError


def calc_emf(model, numComponents):
    componentFluxes = []
    componentFluxErrors = []
    for i in range(numComponents):
        height = model.params['g%d_height' % (i + 1)].value
        sigmaObs = model.params['g%d_sigma' % (i + 1)].value
        heightError = model.params['g%d_height' % (i + 1)].stderr
        sigmaObsError = model.params['g%d_sigma' % (i + 1)].stderr
        amplitude = model.params['g%d_amplitude' % (i + 1)].value
        amplitudeError = model.params['g%d_amplitude' % (i + 1)].stderr
        calcFlux, calcFluxError = amplitude, amplitudeError

        componentFluxes.append(calcFlux)
        componentFluxErrors.append(calcFluxError)
    componentFluxes = np.array(componentFluxes)
    componentFluxErrors = np.array(componentFluxErrors)
    totalFlux = sum(componentFluxes)
    totalFluxError = np.sqrt(sum([s**2 for s in componentFluxErrors]))

    calcEMF = componentFluxes/sum(componentFluxes) * 100

    return calcEMF, componentFluxes, componentFluxErrors, totalFlux, totalFluxError


def calc_continuum(model, emName, rp):
    continuumList = model.best_values['lin_slope'] * np.array(rp.emProfiles[emName]['centerList']) + model.best_values['lin_intercept']
    globalContinuum = model.best_values['lin_slope'] * np.mean(rp.emProfiles[emName]['centerList']) + model.best_values['lin_intercept']

    return continuumList, globalContinuum


def calc_luminosity(rp):
    if 'H-Alpha' in rp.emProfiles:
        calcLuminosity = (4 * np.pi * rp.emProfiles['H-Alpha']['globalFlux'] * (
                    rp.distance * 3.086) ** 2 / rp.scaleFlux) * 1e48
        calcLuminosityError = (4 * np.pi * (rp.distance * 3.086) ** 2 * rp.emProfiles['H-Alpha'][
            'globalFluxErr'] / rp.scaleFlux) * 1e48
    else:
        calcLuminosity = (4 * np.pi * rp.emProfiles['H1r_6563A']['globalFlux'] * (
                    rp.distance * 3.086) ** 2 / rp.scaleFlux) * 1e48
        calcLuminosityError = (4 * np.pi * (rp.distance * 3.086) ** 2 * rp.emProfiles['H1r_6563A'][
            'globalFluxErr'] / rp.scaleFlux) * 1e48

    starFormRate = 7.9e-42 * calcLuminosity  # (*1e6 M_sun/year) by Kennicutt, 1998, 0.1-100 M_sun, Te=10^4k, ne=100 cm^-3
    starFormRateError = 7.9e-42 * calcLuminosityError

    #starFormRate = 5.5e-42 * calcLuminosity  # (*1e6 M_sun/year) by Starburst99, Leitherer et al. 1999, 0.1-100 M_sun, Te=10^4k, ne=100 cm^-3
    #starFormRateError = 5.5e-42 * calcLuminosityError  #(*1e6 M_sun/year)

    logLuminosity = np.log10(calcLuminosity)
    logLuminosityError = 0.434 * calcLuminosityError / calcLuminosity

    return logLuminosity, logLuminosityError, starFormRate, starFormRateError


def save_measurements(measurementInfo, rp):<commit_msg>f deleted of the print(f"example_region.add_em_line(name='{emName}'<commit_after>import csv
import os
import sys
import numpy as np
from fitelp.label_tools import line_label
from fitelp.read_spectra import GalaxyRegion
from fitelp.fit_line_profiles import FittingProfile, plot_profiles, vel_to_wave, wave_to_vel
from fitelp.make_latex_tables import average_velocities_table_to_latex, halpha_regions_table_to_latex, comp_table_to_latex
from fitelp.bpt_plotting import calc_bpt_points, bpt_plot
import fitelp.constants as constants

constants.init()
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '../scripts'))
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '../Input_Galaxy_Region_Information'))
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '../Input_Data_Files'))


def calc_vel_dispersion(sigmaObs, sigmaObsError, sigmaTemp2, filter, rp, correctCopiedSigma=False):
    # Assuming negligible error in temp or instrument
    if filter == 'blue':
        sigmaInstr = rp.sigmaInstrBlue
    elif filter == 'red':
        sigmaInstr = rp.sigmaInstrRed

    #Old version (added by Vero)
    # totalSigmaSquared = sigmaObs ** 2 - sigmaInstr ** 2 - sigmaTemp2
    # totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
    # try:
    #     intrinsic = np.sqrt(totalSigmaSquared)
    #     intrinsicError = 0.5 * totalSigmaSquared ** (-0.5) * totalSigmaSquaredError
    # except ValueError:
    #     "ERROR: INVALID SIGMA"
    #     intrinsic = 0
    #     intrinsicError = 0
    #
    # Last version (added by Vero)
    #if sigmaObs**2 > sigmaTemp2:
        # totalSigmaSquared = sigmaObs ** 2 - sigmaInstr ** 2 - sigmaTemp2
        # totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
    # else:
    #     totalSigmaSquared = (sigmaInstr + np.sqrt(sigmaTemp2))**2
    #     totalSigmaSquaredError = 0
    # intrinsic = np.sqrt(totalSigmaSquared)
    # intrinsicError = 0.5 * totalSigmaSquared ** (-0.5) * totalSigmaSquaredError
    #
    if correctCopiedSigma:
        # SigmaObs here is the SigmaIntrinsic of the copied line
        totalSigmaSquared = sigmaObs**2 + sigmaInstr**2 + sigmaTemp2
        totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
    else:
        if sigmaObs**2 > (sigmaTemp2 + sigmaInstr**2):
            totalSigmaSquared = sigmaObs**2 - sigmaInstr**2 - sigmaTemp2
            totalSigmaSquaredError = 2 * sigmaObs * sigmaObsError
        else:
            print("ERROR: INVALID SIGMA")
            totalSigmaSquared = 1e-99
            totalSigmaSquaredError = 0
    intrinsic = np.sqrt(totalSigmaSquared)
    intrinsicError = 0.5 * totalSigmaSquared**(-0.5) * totalSigmaSquaredError

    return intrinsic, intrinsicError


def calc_flux(height, sigmaObs, heightError, sigmaObsError):
    insideSqrt = 2 * np.pi * sigmaObs**2
    insideSqrtError = 2 * sigmaObs * sigmaObsError
    sqrt = np.sqrt(insideSqrt)
    sqrtError = 0.5 * insideSqrt**(-0.5) * insideSqrtError

    calcFlux = height * sqrt
    calcFluxError = calcFlux * np.sqrt((heightError/height)**2 + (sqrtError/sqrt)**2)

    return calcFlux, calcFluxError


def calc_emf(model, numComponents):
    componentFluxes = []
    componentFluxErrors = []
    for i in range(numComponents):
        height = model.params['g%d_height' % (i + 1)].value
        sigmaObs = model.params['g%d_sigma' % (i + 1)].value
        heightError = model.params['g%d_height' % (i + 1)].stderr
        sigmaObsError = model.params['g%d_sigma' % (i + 1)].stderr
        amplitude = model.params['g%d_amplitude' % (i + 1)].value
        amplitudeError = model.params['g%d_amplitude' % (i + 1)].stderr
        calcFlux, calcFluxError = amplitude, amplitudeError

        componentFluxes.append(calcFlux)
        componentFluxErrors.append(calcFluxError)
    componentFluxes = np.array(componentFluxes)
    componentFluxErrors = np.array(componentFluxErrors)
    totalFlux = sum(componentFluxes)
    totalFluxError = np.sqrt(sum([s**2 for s in componentFluxErrors]))

    calcEMF = componentFluxes/sum(componentFluxes) * 100

    return calcEMF, componentFluxes, componentFluxErrors, totalFlux, totalFluxError


def calc_continuum(model, emName, rp):
    continuumList = model.best_values['lin_slope'] * np.array(rp.emProfiles[emName]['centerList']) + model.best_values['lin_intercept']
    globalContinuum = model.best_values['lin_slope'] * np.mean(rp.emProfiles[emName]['centerList']) + model.best_values['lin_intercept']

    return continuumList, globalContinuum


def calc_luminosity(rp):
    if 'H-Alpha' in rp.emProfiles:
        calcLuminosity = (4 * np.pi * rp.emProfiles['H-Alpha']['globalFlux'] * (
                    rp.distance * 3.086) ** 2 / rp.scaleFlux) * 1e48
        calcLuminosityError = (4 * np.pi * (rp.distance * 3.086) ** 2 * rp.emProfiles['H-Alpha'][
            'globalFluxErr'] / rp.scaleFlux) * 1e48
    else:
        calcLuminosity = (4 * np.pi * rp.emProfiles['H1r_6563A']['globalFlux'] * (
                    rp.distance * 3.086) ** 2 / rp.scaleFlux) * 1e48
        calcLuminosityError = (4 * np.pi * (rp.distance * 3.086) ** 2 * rp.emProfiles['H1r_6563A'][
            'globalFluxErr'] / rp.scaleFlux) * 1e48

    starFormRate = 7.9e-42 * calcLuminosity  # (*1e6 M_sun/year) by Kennicutt, 1998, 0.1-100 M_sun, Te=10^4k, ne=100 cm^-3
    starFormRateError = 7.9e-42 * calcLuminosityError
    #starFormRate = 5.5e-42 * calcLuminosity  # (*1e6 M_sun/year) by Starburst99, Leitherer et al. 1999, 0.1-100 M_sun, Te=10^4k, ne=100 cm^-3
    #starFormRateError = 5.5e-42 * calcLuminosityError  #(*1e6 M_sun/year)

    logLuminosity = np.log10(calcLuminosity)
    logLuminosityError = 0.434 * calcLuminosityError / calcLuminosity

    return logLuminosity, logLuminosityError, starFormRate, starFormRateError


def save_measurements(measurementInfo, rp):
{'fraction_non_alphanumeric': 0.11646724627700783, 'fraction_numerical': 0.02978393733321856, 'mean_word_length': 3.525905726529022, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2210639', 'n_tokens_mistral': 4393, 'n_tokens_neox': 3891, 'n_words': 775}