Commit 59fdfb0f authored by Ruben Farinelli's avatar Ruben Farinelli
Browse files

Changed the computation method for the smoothed polarization degree as a function of energy, 

and fixed a bug in the main script where ktbb was incorrectly kept constant instead of varying with the grid values
parent 303e9971

array_binning.c

+2 −2
Original line number Diff line number Diff line
@@ -8,11 +8,11 @@ void make_array_u(double array_u[], long nx, double start, double end) 
        if (nx == 1) array_u[0] = start;

        return;

    }

    double step = (end - start) / (nx - 1);

    double step = (end - start) / nx;

    

    for (long i = 0; i < nx; i++) 

    {

        array_u[i] = start + i * step;

        array_u[i] = start + step/2+ i * step;

    }

}

build_model/make_model.py

+218 −56
Original line number Diff line number Diff line 
import numpy as np



from model_utilities import *

import os

import os,sys

from scipy.optimize import curve_fit



#=======================================================================
@@ -16,11 +16,13 @@ except ImportError: 
    # Your code to handle the failure goes here





TableXspec = table()



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

#Set table descriptors and the energy array

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



TableXspec = table()



TableXspec.setModelName("Polcomp")

TableXspec.setModelUnits(" ")

TableXspec.setisRedshift(True)
@@ -32,12 +34,11 @@ filters = ["Stokes:0", "Stokes:1", "Stokes:2"] 
TableXspec.setFiltExps(filters)

FlagKeywords=True





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

#Here set the number of model parameters

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



TableXspec.setNumIntParams(4)

TableXspec.setNumIntParams(5)

TableXspec.setNumAddParams(0)



#====================================================
@@ -47,7 +48,7 @@ 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):



    min_val=array_values[0]
@@ -56,7 +57,6 @@ def table_init(init_val=None, array_values=None, name=None, table=None): 
    testpar = tableParameter()

    testpar.setName(name)





    testpar.setInterpolationMethod(0)

    testpar.setInitialValue(init_val)

    testpar.setDelta(0.1)
@@ -74,6 +74,85 @@ 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):



    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)



    Q_matrix=-Q_matrix



    #Q_filtered = smooth_matrix_Q(specpol_array=Q_matrix, energy_center=ene_center,

    #                             delta_energy=delta_energy)



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

    # Write the I,Q,U spectra

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





    PD_filtered = smooth_matrix_pd(Q_matrix=Q_matrix, C_matrix=C_matrix,

                                   energy_center=ene_center, emin=0.3, emax=30)



    Norm_flux = normalize_counts(counts_array=C_matrix, energy_center=ene_center, delta_energy=delta_energy,

                                 u_array=u_centers, delta_u=delta_u, ktbb=float(ktbb))



    I_matrix = intensity_matrix(counts_array=C_filtered, energy_center=ene_center, delta_energy=delta_energy,

                                delta_u=delta_u, norm_value=Norm_flux)



    #print('Python: I %5.4e Q %5.4e' % (I_matrix[10][10], Q_filtered[10][10]))



    if fitsfile is None:

        fitsfile = 'intensity_pd.fits'



    try:

        create_fits_image(I_matrix=I_matrix, PD_matrix=PD_filtered, emin=emin, emax=emax, filename=fitsfile)

    except:

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

        os.sys.exit(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))





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



        filename = 'polarimetry.dat'



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

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



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

    if xspec is True:



        for tt in range(3):



            testspec = tableSpectrum()

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



            if tt == 0:

                testspec.setFlux(I_flux)

            elif tt == 1:

                testspec.setFlux(Q_flux)

            else:

                testspec.setFlux(U_flux)



            TableXspec.pushSpectrum(testspec)





    return C_filtered, PD_filtered



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

#Angles with respect to the normal of the slab

#============================================================
@@ -105,51 +184,125 @@ if module_heasp==True: 
    TableXspec.setEnergies(logbins)



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

for inputval in sys.argv:



input_file='../grid0/test_BL_slab_Z0seed_grid__1.5_3.0_0.5__IQ_ADMP.dat'

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

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



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

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

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



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

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



PD_matrix=Q_matrix/C_matrix

PD_matrix = np.nan_to_num(PD_matrix, nan=0.0)

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

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



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

Q_filtered=smooth_matrix_Q(specpol_array=Q_matrix, energy_center=ene_center, delta_energy=delta_energy)

PD_filtered=smooth_matrix_Q(specpol_array=PD_matrix, energy_center=ene_center, delta_energy=delta_energy)

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

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



Norm_flux=normalize_counts(counts_array=C_matrix, energy_center=ene_center, delta_energy=delta_energy,

                 u_array=u_centers, delta_u=delta_u, ktbb=1.5)

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

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



I_matrix=intensity_matrix(counts_array=C_filtered, energy_center=ene_center, delta_energy=delta_energy,

                          delta_u=delta_u, norm_value=Norm_flux)



create_fits_image(I_matrix=I_matrix, PD_matrix=PD_filtered, emin=emin, emax=emax)

try:

    runtable

except:

    runtable=None



try:

    mcfile

except:

    mcfile=None



try:

    ktbb

except:

    ktbb=None



try:

    kte

except:

    kte=None



try:

    tau

except:

    tau=None



try:

    fitsfile

except:

    fitsfile='intensity_pd.fits'



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

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=='n':



    if mcfile is None:



        print('Please select the ASCII file derived from MC simulations used to  build the fits image')

        print('with command mcfile=filename')

        exit()



    if ktbb is None:

        print('\nError: ktbb must be defined in order to compute the normalized specific intensity of the 2D fits file')

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

        os.sys.exit(1)





    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)



exit(1)



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



for ii in range(0,20):

    data = 'I'



    for ii in range(10, 15):



    pdeg=Q_matrix[ii]/C_matrix[ii]

        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)



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

        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)



    #plot_data(x_value=ene_center, y_mc=PD_matrix[ii]*100, y_smooth=PD_filtered[ii]*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)



        else:

            a = 1



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



exit()

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

#Here define the grid values for kTbb, kTe and tau



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



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



if runtable=='n':

    exit(1)



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

#Here define the grid values for kTbb, kTe, tau and i_obs

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



#Latmax

latmax_array = np.array([10., 30., 50., 70., 90.])



#latmax_array = np.array([10.,50])



lat_init=45

table_init(init_val=lat_init, array_values=latmax_array, name='latmax', table=TableXspec)
@@ -170,62 +323,71 @@ table_init(init_val=kte_init, array_values=kte_array, name='kTe', table=TableXsp 
#tau



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





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)



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



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

#Start the loop over the parameters of the model

#cos(i_obs)

uobs_array=np.linspace(0,1, 10)

uobs_init=0.5



table_init(init_val=uobs_init, array_values=uobs_array, name='cos_iobs', table=TableXspec)



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

#NS Spin



for latmax in latmax_array:

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

spin_init=200



    for kte in formatted_kte_val:

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



        for tau in formatted_tau_val:



            for ktbb in formatted_ktbb_val:

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



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

                file_path='../grid0/'+filename

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

#Seed photons flag

#0 is for photons at the base of the BL

#1 is for uniformly distributed photons

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

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

flag_init=0



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



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



                I_flux,Q_flux, U_flux=(

                    read_spectra(ktbb=ktbb, kte=kte, tau=tau, latmax=latmax, array_ene=ene_center))

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

#Start the loop over the parameters of the model

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



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

                #Write the I,Q,U spectra

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

for latmax in latmax_array:



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

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

    for ktbb in formatted_ktbb_val:



                call_modell(fitsfile='random_image.fits', iobs=40, latmax=latmax)

        for kte in formatted_kte_val:



                for tt in range(3):

            for tau in formatted_tau_val:



                    testspec = tableSpectrum()

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

                for u_obs in uobs_array:



                    if tt == 0:

                        testspec.setFlux(I_flux)

                    elif tt == 1:

                        testspec.setFlux(Q_flux)

                    else:

                        testspec.setFlux(U_flux)



                    TableXspec.pushSpectrum(testspec)

                    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)







# now write out the table.

try:

    TableXspec.write(TableName)

except:

    print('Error writing table data')

    exit(1)



filter_fits(filename=TableName, stokes='I')

filter_fits(filename=TableName, stokes='Q')

filter_fits(filename=TableName, stokes='U')

build_model/model_utilities.py

+206 −61

File changed.

Preview size limit exceeded, changes collapsed.

grid1/test_BL_slab_ISOsed_grid__0.5_1.5_0.5__IQ_ADMP.dat

0 → 100644
+5120 −0

File added.

Preview size limit exceeded, changes collapsed.

grid1/test_BL_slab_ISOsed_grid__0.5_1.5_1.0__IQ_ADMP.dat

0 → 100644
+5120 −0

File added.

Preview size limit exceeded, changes collapsed.