New functions implemented to fit the MC spectra but still to be optimized

import numpy as np

from scipy.special import gamma as gammln, zeta

import matplotlib.pyplot as plt

# Seed spectrum: modified blackbody

def bbody(R, z, bb_par):

    return z**bb_par / (np.exp(R * z) - 1)

# Interpolation function (linear)

def interp_funct(x, y, N, val):

    return np.interp(val, x[:N], y[:N])

# Main solver

def xscomptb(energy,  ktbb, kte, alpha, logA, Norm):

    # Unpack parameters

    delta=0

    bb_par=3

    N=4000

    #ktbb, alpha, kte, logA = parameters

    A = 10**logA

    fb = 1.0  # neglecting (vb/c)^2 terms

    R = kte / ktbb

    cn = 8.0525 * (kte / ktbb)**bb_par * ktbb**(bb_par - 4)

    gamma_val = alpha * (alpha + delta + 3) + 3 * delta

    # Grid

    tmin, tmax = -10, 10

    h = (tmax - tmin) / N

    t = tmin + h * np.arange(N+1)

    xnum = np.exp(t)

    # Coefficient arrays

    p = np.full(N+1, fb)

    g = xnum - 3 * fb - delta

    r = xnum - gamma_val + 3 * delta

    s = -bbody(R, xnum, bb_par)

    s[s > -1e-15] = 0.0

    a = p / h**2

    b = -(2 * p / h**2 + g / h - r)

    c = p / h**2 + g / h

    # Boundary conditions

    a[0], b[0], c[0], s[0] = 0, 1, 0, 0

    a[N], b[N], c[N], s[N] = 0, 1, 0, 0

    # Thomas algorithm

    L = np.zeros(N+1)

    K = np.zeros(N+1)

    L[0] = -c[0] / b[0]

    K[0] = 0

    for n in range(1, N+1):

        denom = a[n] * L[n-1] + b[n]

        L[n] = -c[n] / denom

        K[n] = (s[n] - a[n] * K[n-1]) / denom

    # Backward substitution

    spec = np.zeros(N+1)

    integ = 0

    for n in reversed(range(N)):

        spec[n] = L[n] * spec[n+1] + K[n]

        integ += spec[n] * h

        if spec[n] < 0:

            print("Warning: negative flux")

    # Evaluate output spectrum

    Nflux = len(energy)

    flux = np.zeros(Nflux)

    flux_comp = np.zeros(Nflux)

    flux_bb = np.zeros(Nflux)

    unnorm_flux = np.zeros(Nflux)

    bb = np.zeros(Nflux)

    x = energy / kte

    flux_comp = interp_funct(xnum, spec, N, x) / energy

    bb = bbody(R, x, bb_par) / energy

    delta_ene = np.diff(energy)

    unnorm_flux[1:] = 0.5 * (flux_comp[1:] + flux_comp[:-1])

    flux_bb[1:] = 0.5 * (bb[1:] + bb[:-1])

    # Photon number conservation normalization

    phtot = (ktbb**bb_par) * gammln(bb_par) * zeta(bb_par) / (integ * kte**bb_par)

    flux = Norm / (A + 1) * (flux_bb + A * phtot * unnorm_flux)

    #plt.xscale('log')

    #plt.yscale('log')

    #plt.plot(energy, flux, marker='o', markersize=1, label='', linewidth=2)

    #plt.show()

    return flux

import numpy as np

#========================================================================================

def piecewise_bb_cpl_smooth_deriv(E, A, kTbb, Gamma, Ecut, Esplit):

    """

    Funzione spezzata BB+CPL con continuità di valore e derivata in Esplit

    e una sola normalizzazione libera A.

    """

    # Corpo nero (non normalizzato)

    BB = (E**2) / (np.exp(E / kTbb) - 1)

    # Punto di raccordo

    x = Esplit

    BB_split = (x**2) / (np.exp(x / kTbb) - 1)

    # Derivata del BB in Esplit

    expx = np.exp(x / kTbb)

    dBB_dx = (

        2 * x / (expx - 1)

        - (x**2) * (expx / (kTbb * (expx - 1)**2))

    )

    # CPL e derivata in Esplit (senza normalizzazione)def fit_log_poly_cpl(energy_center=None, delta_energy=None, array_y=None, emin=None, emax=None, deg=5,

                     split_energy=20, ktbb=None, kte=None, tau=None):

    #result = minimize_scalar(

    #    lambda split: compute_chisq_for_split(split, energy_center, array_y, emin, deg),

    #    bounds=(4, 30), method='bounded'

    # split_energy=result.x

    #)

    y_err = np.sqrt(np.abs(array_y))

    array_y=array_y/delta_energy

    safe_y_err = np.where(y_err == 0, 1e8, y_err)/delta_energy

    energy_at_max = energy_center[np.argmax(array_y)]

    split_energy=2*energy_at_max

    split_energy=10

    # ==========================================================

    mask_low = (energy_center > emin) & (energy_center <= split_energy) & (array_y > 0)

    log_x = np.log10(energy_center[mask_low])

    log_y = np.log10(array_y[mask_low])

    # Fit polinomiale in scala log-log

    poly_coeffs = np.polyfit(log_x, log_y, deg=deg)

    poly_fit = np.poly1d(poly_coeffs)

    # Valore e derivata a split_energy

    log_split = np.log10(split_energy)

    log_y_split = poly_fit(log_split)

    # Converti in valore reale

    y_split = 10 ** log_y_split

    #==========================================================

    # Mask per regione CPL

    # =========================================================

    mask_high = (energy_center > split_energy) & (array_y > 0) & (array_y/safe_y_err > 2) & (energy_center < 40)

    while (np.sum(mask_high) < 5):

        split_energy=split_energy-1

        mask_high = (energy_center > split_energy) & (array_y > 0) & (array_y / safe_y_err > 2)

    x_high = energy_center[mask_high]

    y_high = array_y[mask_high]

    # Stima iniziale: usa continuità per inizializzare A

    gamma0 = 2

    Ecut0 = split_energy

    if kte is not None:

        Ecut0=3*kte

    A0 = y_split / (split_energy ** (-gamma0) * np.exp(-split_energy / Ecut0))

    bounds=([0, 0.5, 0.1], [np.inf, 5, 50])

    try:

        with warnings.catch_warnings():

            warnings.simplefilter("ignore", OptimizeWarning)

            popt, _ = curve_fit(

                cutoff_powerlaw,

                x_high,

                y_high,

                p0=[A0, gamma0, Ecut0],

                sigma=y_err[mask_high],

                absolute_sigma=True,

                bounds=bounds,

                maxfev=2000

            )

            if not np.all(np.isfinite(popt)):

                return np.zeros(len(array_y), dtype=bool)

            # Costruisci modello combinato

            result_full = np.zeros_like(energy_center)

            print("A %5.3e gamma %5.3f Ecut %5.4f" % (popt[0], popt[1], popt[2]))

            mask_below=(energy_center < split_energy)

            mask_above = (energy_center > split_energy)

            # Parte polinomiale

            result_full[mask_below] = 10 ** poly_fit(np.log10(energy_center[mask_below]))

            #result_full[mask_low]=np.ones(len(energy_center[mask_low]))

            # Parte CPL

            result_full[mask_above] = cutoff_powerlaw(energy_center[mask_above], *popt)

    except RuntimeError:

        # Fit failed (did not converge)

        result_full=np.zeros(len(array_y), dtype=bool)

    # Calcolo del chi²

    chi_squared = np.sum(

        ((array_y - result_full) / safe_y_err) ** 2

    )/(len(array_y)-3)

    #print(f"Chi² del fit: {chi_squared:.2f}")

    return result_full

    CPL_split = x**(-Gamma) * np.exp(-x / Ecut)

    dCPL_dx = CPL_split * (-Gamma / x - 1 / Ecut)

    # Fattore di raccordo

    R = BB_split / CPL_split

    R_prime = dBB_dx / dCPL_dx

    # Per derivabilità: imposta R = R_prime (sistema omogeneo)

    # Sostituisci R con media geometrica (più stabile numericamente)

    R = np.sqrt(np.abs(BB_split / CPL_split * dBB_dx / dCPL_dx))

    CPL = R * E**(-Gamma) * np.exp(-E / Ecut)

    # Spezzata

    model = np.where(E <= Esplit, BB, CPL)

    return A * model

#!/usr/bin/env python3

import numpy as np

from model_utilities import *

import os,sys

from scipy.optimize import curve_fit

from comptb import *

#=======================================================================

try:

#Here set the number of model parameters

#====================================================

FreeParams=5

FreeParams=7

TableXspec.setNumIntParams(FreeParams)

TableXspec.setNumAddParams(0)

#====================================================

TableName='blcomp.fits'

if TableName is not None:

    if os.path.exists(TableName):

        os.remove(TableName)

#========================================================================

def table_init(init_val=None,  array_values=None, name=None, table=None):

    return array_values

#================================================================================================

def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, mc_file=None,

def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, spin=None, mc_file=None,

                      xspec=None, fitsfile=None, paramvalues=None):

    print('Reading file: %s' % (mc_file))

    C_matrix, Q_matrix = read_mc_simulations(mc_file=mc_file, n_i=20, n_j=256)

    C_filtered = smooth_matrix_spectrum(specpol_array=C_matrix, energy_center=ene_center,

                                        delta_energy=delta_energy)

                                        delta_energy=delta_energy, ktbb=ktbb, kte=kte, tau=tau)

    #for ii in range(len(ene_center)):

    #    print("[%5.4f - %5.4f]"% (ene_center[ii]-delta_energy[ii]/2, ene_center[ii]+delta_energy[ii]/2))

    Q_matrix=-Q_matrix

        print('Error while creating the fits file %s ' % (fitsfile))

        os.sys.exit(1)

    if spin is None:

        spin=500

    if spin==0:

        spin=0.1

        print('Set spin to 0.1')

    if u_obs is not None and latmax is not None:

        iobs = np.arccos(u_obs) / np.pi * 180

        print('Calling model with latmax=%4.2f ktbb=%4.2f kte=%4.2f tau=%4.2f iobs=%4.3f\n' %

              (float(latmax), float(ktbb), float(kte), float(tau), iobs))

        print('Calling model with latmax=%4.2f ktbb=%4.2f kte=%4.2f tau=%4.2f iobs=%4.3f spin=%4.3f\n' %

              (float(latmax), float(ktbb), float(kte), float(tau), iobs, spin))

        call_modell(fitsfile=fitsfile, iobs=iobs, latmax=latmax)

        call_modell(fitsfile=fitsfile, iobs=iobs, latmax=latmax, spin=spin)

        filename = 'polarimetry.dat'

    if inputval.startswith('tau') == True:

        tau = float(inputval[inputval.rfind('=') + 1:])

    if inputval.startswith('param') == True:

        param = str(inputval[inputval.rfind('=') + 1:])

    if inputval.startswith('ymin') == True:

        ymin = float(inputval[inputval.rfind('=') + 1:])

    if inputval.startswith('tabname') == True:

        tabname = str(inputval[inputval.rfind('=') + 1:])

try:

    runtable

except:

    fitsfile='intensity_pd.fits'

try:

    param

except:

    param=None

try:

    ymin

except:

    ymin=-5

#====================================================================

if runtable is None:

    print('Please select if the XSPEC table model must bu built or not with command runtable=y|n')

    exit()

if runtable=='y':

    try:

        tabname

    except:

        print('Please select the name of the FITS model name with parameter tabname=<tabname>')

        exit(1)

    TableName = tabname

    if TableName is not None:

        if os.path.exists(TableName):

            os.remove(TableName)

#====================================================================

if runtable=='n':

    if mcfile is None:

        print('Set the temperature with ktbb=value')

        os.sys.exit(1)

    if param is None:

        print('Please select the parameter to be plotted with param=name')

        print('Allowed values are  C, Q or pd')

        exit(1)

    allowed_params=['C', 'pd']

    if param not in allowed_params:

        print('Unknow param value')

        print('Allowed values are ', allowed_params)

        exit(1)

    #=========================================================================

    #Get parametes of the simulation from filename

    # =========================================================================

    ktbb,kte,tau=filename2params(filename=mcfile)

    pdf_name = os.path.basename(mcfile.replace(".dat", ".pdf"))

    pdf_name = pdf_name.replace("test_BL", param)

    C_filtered, PD_filtered = write_table_value(latmax=None, ktbb=ktbb, kte=kte,

                                                    tau=tau, u_obs=None, mc_file=mcfile, xspec=False, fitsfile=fitsfile)

    C_matrix, Q_matrix = read_mc_simulations(mc_file=mcfile, n_i=20, n_j=256)

    plt.figure(figsize=(10, 4))

    plt.figure(figsize=(10, 5))

    data = 'PD'

    for ii in range(1, 2):

    for ii in range(10, 12):

        if data == 'Q':

        if param == 'Q':

            plot_data(x_value=ene_center, y_mc=Q_matrix[ii], y_smooth=Q_filtered[ii],

                      title=input_file, ymin=-5500, ymax=5000, logscale=True)

                      title=input_file, ymin=-5500, ymax=5000, logscale=True, pdf_name=pdf_name)

        elif data == 'PD':

        elif param == 'pd':

            safe_Cmatrix = np.where(C_matrix == 0, 1e-20, C_matrix)

            error = np.sqrt(2 / safe_Cmatrix[ii])

            plot_data(x_value=ene_center, y_mc=PD_matrix[ii] * 100, y_smooth=PD_filtered[ii] * 100, ymin=-5, ymax=5,

                        emin=0.3, emax=30, y_err=error*100)

            plot_data(x_value=ene_center, y_mc=PD_matrix[ii] * 100, y_smooth=PD_filtered[ii] * 100, ymin=ymin, ymax=5,

                        emin=0.3, emax=30, y_err=error*100, y_label='Polarization degree (%)', pdf_name=pdf_name)

        elif data == 'I':

        elif param == 'C':

            plot_data(x_value=ene_center, y_mc=C_matrix[ii], y_smooth=C_filtered[ii],

                      logscale=True, ymin=0.1, ymax=1e6, y_err=np.sqrt(C_matrix[ii]))

                      logscale=True, ymin=0.1, ymax=1e7, y_err=np.sqrt(C_matrix[ii]),

                      delta_energy=delta_energy, y_label='Counts/keV', pdf_name=pdf_name)

            I_flux, Q_flux, U_flux=read_spectra(filename='polarimetry.dat')

            #I_flux, Q_flux, U_flux=read_spectra(filename='polarimetry.dat')

        else:

            a = 1

            print('Unknown parameter name to be plotted')

            exit(1)

    #plt.pause(2)  # Mostra il plot per 0.5 secondi

    #plt.close()  # Chiude la figura

    #sys.exit()

    plt.show()

    print('=================================================================================')

    print('\nIf you want to calculate the I,Q,U spectra for the BL around NS now from the prompt you can type:\n')

    print('bl_pol geom=bl enedep=y enefile=%s iobs=value latmax=value spin=value' % (fitsfile))

mcfile='../grid0/test_BL_slab_Z0seed_grid__2.5_4.0_4.0__IQ_ADMP.dat'

#====================================================================

if runtable=='n':

table_init(init_val=ktbb_init, array_values=ktbb_array, name='kTbb', table=TableXspec)

#kTe

kte_array = np.linspace(1.5,5,8)

#kte_array = np.linspace(1.5,5,8)

kte_array=np.array([1.5,2,2.5,3,3.5,4,4.5,5,6,7])

formatted_kte_val=[f"{val:.1f}" for val in kte_array]

kte_init=3

table_init(init_val=tau_init, array_values=tau_array, name='tau', table=TableXspec)

#print(formatted_tau_val)

#====================================================================

#cos(i_obs)

#====================================================================

#NS Spin

#====================================================================

#spin_array=np.linspace(0,500, 5)

#spin_init=200

#table_init(init_val=spin_init, array_values=spin_array, name='spin', table=TableXspec)

spin_array=np.linspace(0,700, 5)

spin_init=300

spin=500

array_parnames=['latmax', 'ktbb', 'kte', 'tau', 'cos_iobs']

table_init(init_val=spin_init, array_values=spin_array, name='spin', table=TableXspec)

#====================================================================

#Seed photons flag

#0 is for photons at the base of the BL

#1 is for uniformly distributed photons

#Flag for polarized or unpolarized seed photos

#====================================================================

flag_array=np.array([0,1])

flag_array=np.array([0., 1.])

flag_init=0

table_init(init_val=flag_init, array_values=flag_array, name='flag', table=TableXspec)

#array_parnames=['latmax', 'ktbb', 'kte', 'tau', 'cos_iobs']

#====================================================================

#Start the loop over the parameters of the model

                for u_obs in uobs_array:

                    for spin in spin_array:

                        for flag in flag_array:

                            filename='test_BL_slab_Z0seed_grid__'+ktbb+'_'+kte+'_'+tau+'__IQ_ADMP.dat'

                    file_path='../grid0/'+filename

                    paramvalues = np.array([float(latmax), float(ktbb), float(kte), float(tau), u_obs])

                            if flag==0:

                                file_path = '../unpolarized_seed/'+filename

                            else:

                                file_path = '../polarized_seed/'+ filename

                            #paramvalues = np.array([float(latmax), float(ktbb), float(kte), float(tau), u_obs])

                            paramvalues = np.array([float(latmax), float(ktbb), float(kte), float(tau), u_obs, spin, flag])

                            if len(paramvalues) != FreeParams:

                                print('WARNING: the number of free parameters determined by the loop is %d' % (len(paramvalues)))

                                print('while that defined for the command TableXspec.setNumIntParams(N) is %d' % (FreeParams))

                                os.sys.exit(1)

                    write_table_value(latmax=latmax, ktbb=ktbb, kte=kte, tau=tau, u_obs=u_obs, mc_file=file_path,

                            write_table_value(latmax=latmax, ktbb=ktbb, kte=kte, tau=tau, u_obs=u_obs, spin=spin, mc_file=file_path,

                                              xspec=True, paramvalues=paramvalues)

# now write out the table.

try:

    TableXspec.write(TableName)

#!/bin/bash

# Directory contenente i file .dat

DATADIR="../polarized_seed"

# Seleziona casualmente 10 file .dat

FILES=$(find "$DATADIR" -maxdepth 1 -type f -name "*.dat"  | shuf | head -n 10)

# Loop su ciascun file selezionato

for FILE in $FILES; do

    echo "Processing $FILE"

    make_model.py runtable=n mcfile="$FILE" ktbb=1.5 param=C

done