Improved the routines for computing the interpolation of the fits files with simulation data (303e9971) · Commits · polarization / boundary_layer

build_model/make_model.py

+1 −3

Original line number	Diff line number	Diff line
		@@ -118,15 +118,13 @@ C_filtered=smooth_matrix_spectrum(specpol_array=C_matrix, energy_center=ene_cent
		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)


		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)

		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_matrix)
		create_fits_image(I_matrix=I_matrix, PD_matrix=PD_filtered, emin=emin, emax=emax)

build_model/model_utilities.py

+12 −3

Original line number	Diff line number	Diff line
		@@ -342,7 +342,8 @@ def call_modell(fitsfile=None, iobs=None, latmax=None):
		raise FileNotFoundError(f"❌ Il file '{dat_file}' non esiste.")

		#============================================================================
		def create_fits_image(filename="intensity_pd.fits", I_matrix=None, PD_matrix=None):
		def create_fits_image(filename="intensity_pd.fits", I_matrix=None,
		PD_matrix=None, emin=None, emax=None):

		if I_matrix.shape != PD_matrix.shape:
		print("Error the two arrays have different dimensions:")
		@@ -350,14 +351,22 @@ def create_fits_image(filename="intensity_pd.fits", I_matrix=None, PD_matrix=Non
		print("PD_matrix:", PD_matrix.shape)
		exit(1)

		nbin_u, nbin_ene=I_matrix.shape

		# Crea un PrimaryHDU e aggiungi keyword al suo header
		primary = fits.PrimaryHDU()
		primary.header['EMIN'] = (emin, 'Left boundary of the first energy bin')
		primary.header['EMAX'] = (emax, 'Right boundary of the last energy bin')
		primary.header['ENE_CHAN'] = (nbin_ene, 'Number of energy channels')
		primary.header['U_CHAN'] = (nbin_u, 'Number of view angles')

		# Crea due ImageHDU
		hdu1 = fits.ImageHDU(data=I_matrix, name='Stokes I')
		hdu2 = fits.ImageHDU(data=PD_matrix, name='Polarization')

		# PrimaryHDU vuoto (obbligatorio per file FITS validi)
		primary = fits.PrimaryHDU()

		# Costruisci l'HDUL (HDU List) e scrivi il file
		os.remove(filename)
		hdul = fits.HDUList([primary, hdu1, hdu2])
		hdul.writeto(filename, overwrite=True)
		print(f"✔️ File FITS '{filename}' creato con due estensioni di dimensione {I_matrix.shape}.")

headers.h

+10 −2

Original line number	Diff line number	Diff line
		@@ -45,8 +45,14 @@ extern gsl_interp_accel* yacc_2d;

		extern double emin;
		extern double emax;
		extern double emin_center;
		extern double emax_center;


		extern double blrad;
		extern double umin_center;
		extern double umax_center;

		double blackbody(double energy, double ktbb);
		double IntensityRestFrame(double E_rest, double CosAlphaPrime, double ktbb);
		double FunctERest(double E_obs, double delta, double u);
		@@ -93,10 +99,12 @@ double call_numerical_integration(double* function_params, double* xmin, double*
		polar_cap* add_signals (polar_cap* primary_cap, polar_cap* secundary_cap);
		float** read_fits_image(const char* filename, int ext_num, long* nx_out, long* ny_out);

		void read_energy_keywords(const char filename, double emin, double emax, long int ene_chan, long int* u_chan);

		void make_array_u(double array_u[], long nx, double start, double end);

		void make_array_energy(double array_energy[], double delta_ene[], double emin, double emax, int nstep);


		double Ivalue_interp(double E_rest, double CosAlphaPrime);
		double Qvalue_interp(double E_rest, double CosAlphaPrime);

		double IQvalue_interp(double E_rest, double CosAlphaPrime, int flag);

main_program.c

+59 −25

Original line number	Diff line number	Diff line
		@@ -22,6 +22,14 @@ double phase;
		double phase_obs;
		double emin;
		double emax;
		double emin_center;
		double emax_center;

		double blrad;
		double factor;
		double flux_energy;
		double umin_center;
		double umax_center;

		char* output;
		char* timefunct;
		@@ -117,6 +125,10 @@ int main(int argc, char* argv[])
		nsrad = numpar(argv[ii]);
		if (strstr(argv[ii], "compactness=") != NULL)
		compactness = numpar(argv[ii]);
		if (strstr(argv[ii], "compactness=") != NULL)
		compactness = numpar(argv[ii]);
		if (strstr(argv[ii], "blrad=") != NULL)
		blrad = numpar(argv[ii]);
		if (strstr(argv[ii], "nphases=") != NULL)
		nphases = numpar(argv[ii]);
		if (strstr(argv[ii], "simulfile=") != NULL)
		@@ -201,15 +213,34 @@ int main(int argc, char* argv[])
		}
		}

		read_energy_keywords(enefile, &emin, &emax, &nstep_ene, &nstep_u);

		array_ene = malloc(nstep_ene * sizeof(double));
		array_delta_ene = malloc(nstep_ene * sizeof(double));
		make_array_energy(array_ene, array_delta_ene, emin, emax, nstep_ene);
		emin_center = array_ene[0];
		emax_center = array_ene[nstep_ene - 1];

		array_u = malloc(nstep_u * sizeof(double));
		make_array_u(array_u, nstep_u, 0, 1);
		umin_center = array_u[0];
		umax_center = array_u[nstep_u - 1];

		/*

		make_array_u(array_u, nstep_u, -1, 1);
		make_array_energy(array_ene, array_delta_ene, emin, emax, nstep_ene);

		emin=array_ene[0];
		emax=array_ene[nstep_ene-1];
		double delta_u = array_u[1] - array_u[0];

		for (ii = 0; ii < nstep_u; ii++)
		{
		printf("%d u %lf [%5.4f - %5.4f]\n", ii, array_u[ii], array_u[ii] - delta_u / 2, array_u[ii] + delta_u / 2);
		}

		for (ii=0; ii< nstep_ene;ii++)
		{
		//printf("%d [%5.4f - %5.4f]\n", ii, array_ene[ii]-array_delta_ene[ii]/2, array_ene[ii]+array_delta_ene[ii]/2);
		printf("%d %5.4f\n", ii, array_ene[ii]);
		}*/

		gsl_spline2d_init(Spline_energy_angle_I, array_ene, array_u, I_upstream_za, nstep_ene, nstep_u);
		gsl_spline2d_init(Spline_energy_angle_Q, array_ene, array_u, Q_upstream_za, nstep_ene, nstep_u);
		@@ -228,7 +259,6 @@ int main(int argc, char* argv[])
		array_Qobs = malloc(nstep_ene * sizeof(double));
		array_Uobs = malloc(nstep_ene * sizeof(double));


		/*for (kk=0; kk< nstep_ene; kk++)
		{
		printf("Ene %lf +/ [%lf = %lf]\n", array_ene[kk], array_ene[kk]-array_delta_ene[kk]/2, array_ene[kk]+array_delta_ene[kk]/2);
		@@ -260,6 +290,8 @@ int main(int argc, char* argv[])
		nphases = 2;
		}

		flux_energy = 0;

		/*for (kk=0; kk< nstep_ene; kk++)
		{
		printf("LUCIO %lf\n", array_ene[kk]);
		@@ -358,6 +390,11 @@ int main(int argc, char* argv[])
		compactness = 1 / 2.5;
		}

		if (blrad == 0)
		{
		blrad = 1.2 * nsrad;
		}

		/=====================================================================/
		if (simulfile == NULL \|\| simulfile == 0x0)
		{
		@@ -471,7 +508,7 @@ int main(int argc, char* argv[])
		for (kk = 0; kk < nstep_ene; kk++)
		{
		function_params[3] = array_ene[kk];
		printf("Results computed for Ene %lf\n", function_params[3]);
		//printf("Results computed for Ene %lf\n", function_params[3]);

		for (ii = 0; ii < ncall; ii++)
		{
		@@ -510,7 +547,6 @@ int main(int argc, char* argv[])
		xmax_secondary[0] = PI - xmin[0];
		xmax_secondary[1] = xmax[1] + PI;


		secundary_cap->Ival[ii] =
		call_numerical_integration(function_params, xmin_secondary, xmax_secondary, STOKES_I);

		@@ -531,7 +567,6 @@ int main(int argc, char* argv[])
		secundary_cap->PD[ii] = sqrt(pow(secundary_cap->Qval[ii], 2) + pow(secundary_cap->Uval[ii], 2)) /
		secundary_cap->Ival[ii];


		secundary_cap->PA[ii] = PA_location(secundary_cap->Qval[ii], secundary_cap->Uval[ii]);
		}
		}
		@@ -597,13 +632,15 @@ int main(int argc, char* argv[])
		}
		else
		{
		fprintf(dat, "%5.4f %5.4f %5.4e %4.3e %4.3e\n", array_ene[kk]-array_delta_ene[kk]/2,
		array_ene[kk]+array_delta_ene[kk]/2,
		primary_cap->Ival[ii], primary_cap->Qval[ii], primary_cap->Uval[ii]);
		}

		factor = pow(blrad, 2);
		PA_obs=PA_location(primary_cap->Qval[ii], primary_cap->Uval[ii]);

		fprintf(dat, "%5.4f %5.4f %5.4e %4.3e %4.3e %5.4f\n", array_ene[kk] - array_delta_ene[kk] / 2,
		array_ene[kk] + array_delta_ene[kk] / 2, primary_cap->Ival[ii] * factor,
		primary_cap->Qval[ii] * factor, primary_cap->Uval[ii] * factor, PA_obs/PI*180);

		flux_energy = flux_energy + primary_cap->Ival[ii] * factor * array_delta_ene[kk];
		}

		} /End of loop over phases (cicle for)/

		@@ -611,10 +648,6 @@ int main(int argc, char* argv[])

		if (GeometryFlag == BOUNDARY_LAYER \|\| (GeometryFlag == POLAR_CAP && FlagMulti == FALSE))
		{
		/* printf("Emission properties\n");
		printf("I %5.4e PD %5.4f PA %4.3e\n", primary_cap->Ival[0], primary_cap->PD[0] * 100,
		primary_cap->PA[0] / PI * 180);
		*/


		if (GeometryFlag == POLAR_CAP)
		@@ -625,6 +658,8 @@ int main(int argc, char* argv[])
		}
		}
		}

		printf("Energy flux %5.4e\n", flux_energy);
		/===========================================================================/

		if (GeometryFlag == POLAR_CAP && FlagMulti == TRUE)
		@@ -894,8 +929,7 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do

		if (FlagEneDep == TRUE)
		{
		I_rest = Ivalue_interp(E_rest, cos_alpha_prime);
		Q_rest = Qvalue_interp(E_rest, cos_alpha_prime);
		I_rest = IQvalue_interp(E_rest, cos_alpha_prime, STOKES_I);
		}
		else
		{
		@@ -923,7 +957,7 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do

		if (FlagEneDep == TRUE)
		{
		PD = Qvalue_interp(E_rest, cos_alpha_prime);
		PD = IQvalue_interp(E_rest, cos_alpha_prime, STOKES_Q);
		}
		else
		{

read_spectropolarimetry.c

+61 −0

Original line number	Diff line number	Diff line
		@@ -61,3 +61,64 @@ float** read_fits_image(const char* filename, int ext_num, long* nx_out, long* n
		*ny_out = ny;
		return image;
		}


		void read_energy_keywords(const char filename, double emin, double emax, long int ene_chan, long int* u_chan)
		{
		fitsfile *fptr; // FITS file pointer
		int status = 0;
		int anynul;
		double val;

		// Apri il file FITS in sola lettura
		if (fits_open_file(&fptr, filename, READONLY, &status)) {
		fits_report_error(stderr, status);
		exit(status);
		}

		// Leggi la keyword EMIN
		if (fits_read_key(fptr, TDOUBLE, "EMIN", emin, NULL, &status)) {
		fprintf(stderr, "❌ Errore nel leggere EMIN dal file FITS\n");
		fits_report_error(stderr, status);
		exit(status);
		}

		// Leggi la keyword EMAX
		if (fits_read_key(fptr, TDOUBLE, "EMAX", emax, NULL, &status)) {
		fprintf(stderr, "❌ Errore nel leggere EMAX dal file FITS\n");
		fits_report_error(stderr, status);
		exit(status);
		}

		// Leggi la keyword EMAX
		if (fits_read_key(fptr, TINT, "ENE_CHAN", ene_chan, NULL, &status)) {
		fprintf(stderr, "❌ Errore nel leggere ENE_CHAN dal file FITS\n");
		fits_report_error(stderr, status);
		exit(status);
		}
		if (fits_read_key(fptr, TINT, "U_CHAN", u_chan, NULL, &status)) {
		fprintf(stderr, "❌ Errore nel leggere U_CHAN dal file FITS\n");
		fits_report_error(stderr, status);
		exit(status);
		}







		// Chiudi il file FITS
		fits_close_file(fptr, &status);
		if (status) {
		fits_report_error(stderr, status);
		exit(status);
		}

		/*printf("✔️ Letti dal file FITS:\n");
		printf(" EMIN = %.3f keV\n", *emin);
		printf(" EMAX = %.3f keV\n", emax);/
		}