Commit 62d83f0e authored by Ruben Farinelli's avatar Ruben Farinelli
Browse files

Added modification to the C program so that when the rest frame energy is above the maximum 

value determined by the Montecarlo grid it returns zero for the specific intensity and polarization
parent 59fdfb0f

build_model/make_model.py

+54 −29
Original line number Diff line number Diff line
@@ -16,6 +16,18 @@ except ImportError: 
    # Your code to handle the failure goes here





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

#Energy bins

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



energy_bins=257

emin=0.3

emax=70



logbins = np.logspace(np.log10(emin), np.log10(emax), num=energy_bins)



ene_center=np.asarray([round((logbins[ii]+logbins[ii+1])/2,4) for ii in range(energy_bins-1)])

delta_energy=np.asarray([round(logbins[ii+1]-logbins[ii],4) for ii in range(energy_bins-1)])



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

#Set table descriptors and the energy array
@@ -23,7 +35,7 @@ except ImportError: 


TableXspec = table()



TableXspec.setModelName("Polcomp")

TableXspec.setModelName("blcomp")

TableXspec.setModelUnits(" ")

TableXspec.setisRedshift(True)

TableXspec.setisAdditive(True)
@@ -38,7 +50,9 @@ FlagKeywords=True 
#Here set the number of model parameters

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



TableXspec.setNumIntParams(5)

FreeParams=5



TableXspec.setNumIntParams(FreeParams)

TableXspec.setNumAddParams(0)



#====================================================
@@ -75,7 +89,8 @@ 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, xspec=None, fitsfile=None):

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

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



    print('Reading file: %s' % (mc_file))
@@ -130,8 +145,7 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, mc 
        #Here read the I,Q,U spectra of the C program

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



        I_flux, Q_flux, U_flux = (read_spectra(ktbb=ktbb, kte=kte, tau=tau, latmax=latmax,

                     array_ene=ene_center, filename=filename))

        I_flux, Q_flux, U_flux = read_spectra(array_ene=ene_center, filename=filename)



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

    if xspec is True:
@@ -139,7 +153,10 @@ def write_table_value(latmax=None, ktbb=None, kte=None, tau=None, u_obs=None, mc 
        for tt in range(3):



            testspec = tableSpectrum()

            testspec.setParameterValues(np.array([float(latmax), float(ktbb), float(kte), float(tau), u_obs]))

            #testspec.setParameterValues(np.array([float(latmax), float(ktbb), float(kte), float(tau), u_obs]))



            testspec.setParameterValues(paramvalues)





            if tt == 0:

                testspec.setFlux(I_flux)
@@ -164,18 +181,6 @@ u_centers = 0.5 * (u_bin_edges[:-1] + u_bin_edges[1:]) 
#for ii in range(len(u_bin_edges)-1):

#    print('Bin %d [%4.2f : %4.2f]' % (ii, u_bin_edges[ii], u_bin_edges[ii+1]))



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

#Energy bins

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



energy_bins=257

emin=0.3

emax=70



logbins = np.logspace(np.log10(emin), np.log10(emax), num=energy_bins)



ene_center=np.asarray([round((logbins[ii]+logbins[ii+1])/2,4) for ii in range(energy_bins-1)])

delta_energy=np.asarray([round(logbins[ii+1]-logbins[ii],4) for ii in range(energy_bins-1)])



#for ii in range(energy_bins-1):

#    print('Channel %d [%5.3f - %5.3f]' % (ii, ene_center[ii]-delta[ii]/2, ene_center[ii]+delta[ii]/2 ))
@@ -262,20 +267,31 @@ if runtable=='n': 


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



    data = 'I'

    data = 'PD'



    for ii in range(10, 15):

    for ii in range(10, 12):



        if data == '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)



        elif data == 'PD':

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

                        emin=0.3, emax=30)



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



            PD_matrix = -Q_matrix / safe_Cmatrix



            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)



        elif data == 'I':

            plot_data(x_value=ene_center, y_mc=C_matrix[ii], y_smooth=C_filtered[ii], logscale=True, ymin=1, ymax=1e6)

            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]))



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





        else:

            a = 1
@@ -324,12 +340,14 @@ table_init(init_val=kte_init, array_values=kte_array, name='kTe', table=TableXsp 


tau_array = np.linspace(0.5,4, 8)



tau_array=np.array([0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6])



formatted_tau_val=[f"{val:.1f}" for val in tau_array]

tau_init=4



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



#print(formatted_tau_val)



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

#cos(i_obs)
@@ -341,10 +359,10 @@ table_init(init_val=uobs_init, array_values=uobs_array, name='cos_iobs', table=T 
#====================================================================

#NS Spin



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

spin_init=200

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

#spin_init=200



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

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





array_parnames=['latmax', 'ktbb', 'kte', 'tau', 'cos_iobs']
@@ -357,8 +375,6 @@ array_parnames=['latmax', 'ktbb', 'kte', 'tau', 'cos_iobs'] 
flag_array=np.array([0,1])

flag_init=0







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

#Start the loop over the parameters of the model

#====================================================================
@@ -377,7 +393,16 @@ for latmax in latmax_array: 
                    filename='test_BL_slab_Z0seed_grid__'+ktbb+'_'+kte+'_'+tau+'__IQ_ADMP.dat'

                    file_path='../grid0/'+filename



                    write_table_value(latmax=latmax, ktbb=ktbb, kte=kte, tau=tau, u_obs=u_obs, mc_file=file_path, xspec=True)

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



                    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,

                                      xspec=True, paramvalues=paramvalues)

build_model/model_utilities.py

+11 −3
Original line number Diff line number Diff line
@@ -173,7 +173,7 @@ def angular_integration(u_array=None, Y_array=None): 


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

def plot_data(x_value=None, y_mc=None, y_smooth=None, title='', ymin=None, ymax=None, emin=None, emax=None,

              logscale=None):

              logscale=None, y_err=None):



    if emin is not None:

        plt.xlim(left=emin)
@@ -192,13 +192,21 @@ def plot_data(x_value=None, y_mc=None, y_smooth=None, title='', ymin=None, ymax= 


    plt.xlabel("Energy (keV)")



    plt.xscale('log')



    if logscale is True:

        plt.xscale('log')

        plt.yscale('log')



    if y_mc is not  None:





        if y_err is None:

            plt.plot(x_value, y_mc, marker='s', markersize=2, label='', linewidth=0.5)



        else:

            plt.errorbar(x_value, y_mc, yerr=y_err, fmt='s', markersize=2, linewidth=1.5, capsize=2, label='Monte Carlo')



    if y_smooth is not None:

        plt.plot(x_value, y_smooth, marker='o', markersize=1, label='')
@@ -213,7 +221,7 @@ def model(x, N, alpha, xc): 


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



def read_spectra(ktbb=None, kte=None, tau=None, latmax=None, array_ene=None, filename=None):

def read_spectra(array_ene=None, filename=None):



    if filename is not  None:

create_random_image.py

0 → 100755
+22 −0
Original line number Diff line number Diff line 
import numpy as np

from astropy.io import fits



def create_random_fits(filename="random_image.fits", shape=(40, 100)):

    # Crea due array random 40x100

    data1 = np.random.rand(*shape)

    data2 = -1*np.random.rand(*shape)



    # Crea due ImageHDU

    hdu1 = fits.ImageHDU(data=data1, name='RANDOM1')

    hdu2 = fits.ImageHDU(data=data2, name='RANDOM2')



    # PrimaryHDU vuoto (obbligatorio per file FITS validi)

    primary = fits.PrimaryHDU()



    # Costruisci l'HDUL (HDU List) e scrivi il file

    hdul = fits.HDUList([primary, hdu1, hdu2])

    hdul.writeto(filename, overwrite=True)

    print(f"✔️ File FITS '{filename}' creato con due estensioni di dimensione {shape}.")



# Esegui

create_random_fits()

grid0/test_BL_slab_Z0seed_grid__0.5_1.5_4.5__IQ_ADMP.dat

0 → 100644
+5120 −0

File added.

Preview size limit exceeded, changes collapsed.

grid0/test_BL_slab_Z0seed_grid__0.5_1.5_5.0__IQ_ADMP.dat

0 → 100644
+5120 −0

File added.

Preview size limit exceeded, changes collapsed.