Starting update where the case of two polar caps is considered (0ab1c3ec) · Commits · polarization / boundary_layer

compile.txt

+1 −1

gcc -O3 -Wall -o bl_pol main_program.c hcubature.c restframe_quantities.c polarization_angle.c numpar.c stringpar.c strrev.c interpolate_functions.c trig_funct.c read_file.c -lm

gcc -O3  -o bl_pol main_program.c hcubature.c restframe_quantities.c polarization_angle.c numpar.c stringpar.c strrev.c interpolate_functions.c trig_funct.c read_file.c sum_values.c -I $HEADERS_GSL -lm -L $LIB_GSL -lgsl -lgslcblas

headers.h

+16 −0

Original line number	Diff line number	Diff line
		@@ -5,6 +5,10 @@
		#include <stdint.h>
		#include <string.h>
		#include "cubature.h"
		#include "structures.h"
		#include <gsl/gsl_interp.h>
		#include <gsl/gsl_spline.h>


		#define TRUE 1
		#define FALSE 0
		@@ -16,6 +20,14 @@
		#define PI 3.14159
		#define enep 2.73

		#define BOUNDARY_LAYER 100
		#define POLAR_CAP 101

		#define STOKES_I 1
		#define STOKES_Q 2
		#define STOKES_U 3


		extern int tot_elements;
		extern int FlagFile;
		extern char* poltype;
		@@ -61,9 +73,13 @@ double reduce_azim(double phi);
		void Calculate_InterpolatedIntensity(double target_costheta, double* interpolated_intensity, double* interpolated_pdeg);
		void read_file(void);

		double PA_location(double Qval, double Uval);
		void interpolate_values(double target_costheta,
		double* Array_costheta,
		double* Array_pdeg,
		double* Array_Intensity,
		double* interpolated_pdeg,
		double* interpolated_intensity);

		double call_numerical_integration(double* function_params, double* xmin, double* xmax, int stokes_flag);
		polar_cap* add_signals (polar_cap* primary_cap, polar_cap* secundary_cap);

interpolate_functions.c

+10 −15

Original line number	Diff line number	Diff line
		@@ -12,8 +12,6 @@ void interpolate_values(double target_costheta,
		int i;
		for (i = 0; i < tot_elements - 1; i++)
		{


		if (target_costheta >= Array_costheta[i] && target_costheta <= Array_costheta[i + 1])
		{
		// printf("%lf [%lf - %lf]\n", target_costheta, Array_costheta[i], Array_costheta[i+1]);
		@@ -31,12 +29,9 @@ void interpolate_values(double target_costheta,
		}
		}

		void Calculate_InterpolatedIntensity(double target_costheta, double* interpolated_intensity,
		double* interpolated_pdeg)
		void Calculate_InterpolatedIntensity(double target_costheta, double* interpolated_intensity, double* interpolated_pdeg)
		{

		// printf("Target %lf %d\n", target_costheta, tot_elements);
		interpolate_values(target_costheta, Array_costheta, Array_pdeg, Array_Intensity,
		interpolated_pdeg, interpolated_intensity);

		}

main_program.c

+328 −55

Original line number	Diff line number	Diff line
		@@ -6,34 +6,58 @@ double phimin;
		double phimax;
		double i_obs;
		double spin;
		double nsrad;
		double compactness;

		char* poltype;
		char* output;
		char* timefunct;
		char* geom;
		char* simulfile;
		int FlagFile;
		int npt;
		int nphases;
		int GeometryFlag;

		int main(int argc, char* argv[])
		{
		static double xmin[2];
		static double xmax[2];
		FILE* dat;
		double xmin[2];
		double xmax[2];

		double xmin_secondary[2];
		double xmax_secondary[2];

		double function_params[15];

		double latmin_rad;
		double latmax_rad;
		double lat_center;

		double compactness;
		double E_obs;
		double ktbb;
		double PD_obs;
		double PA_obs;
		double ratio_UQ;
		double PD_obs_secundary;
		double PA_obs_secundary;

		double stokes1 = 0, err_stokes1 = 0;
		double Qval = 0, Qval_err = 0;
		double Uval = 0, Uval_err = 0;
		double deltaphase_primary;
		double deltaphase_secundary;

		double phi_step;

		double Ival;
		double Qval;
		double Uval;

		double Ival_secundary;
		double Qval_secundary;
		double Uval_secundary;

		FlagFile = FALSE;

		int ii;
		int ncall;
		int FlagMulti = FALSE;

		for (ii = 1; ii < argc; ii++)
		{
		@@ -49,10 +73,44 @@ int main(int argc, char* argv[])
		i_obs = numpar(argv[ii]);
		if (strstr(argv[ii], "spin=") != NULL)
		spin = numpar(argv[ii]);
		if (strstr(argv[ii], "nsrad=") != NULL)
		nsrad = numpar(argv[ii]);
		if (strstr(argv[ii], "compactness=") != NULL)
		compactness = numpar(argv[ii]);
		if (strstr(argv[ii], "nphases=") != NULL)
		nphases = numpar(argv[ii]);
		if (strstr(argv[ii], "simulfile=") != NULL)
		simulfile = stringpar(argv[ii]);
		if (strstr(argv[ii], "poltype=") != NULL)
		poltype = stringpar(argv[ii]);
		if (strstr(argv[ii], "geom=") != NULL)
		geom = stringpar(argv[ii]);
		if (strstr(argv[ii], "timefunct=") != NULL)
		timefunct = stringpar(argv[ii]);
		if (strstr(argv[ii], "npt=") != NULL)
		npt = numpar(argv[ii]);
		if (strstr(argv[ii], "output=") != NULL)
		output = stringpar(argv[ii]);
		}

		/=============================================================================================/

		if (geom == NULL)
		{
		printf(
		"Plese select the geometry of the region: it can be geom=bl\|xrp for boundary layer or "
		"X-ray pulsar (or FRB), respectively \n");
		return 1;
		}

		if (output == NULL)
		{
		output = "polarimetry.dat";
		}

		dat = fopen(output, "w"); // Apri in modalità scrittura
		if (dat == NULL)
		{
		printf("Errore while opening file %s\n", output);
		return 1;
		}

		if (!latmax \|\| latmax == 0x0)
		@@ -67,37 +125,92 @@ int main(int argc, char* argv[])
		return 1;
		}

		if (!latmin \|\| latmin == 0x0)
		if (nphases == 0)
		{
		latmin = 0;
		nphases = 2;
		}

		if (!phimin \|\| phimin == 0x0)
		/=======================================================================================/

		if (strcmp(geom, "bl") == 0)
		{
		GeometryFlag = BOUNDARY_LAYER;

		latmin = 0;
		phimin = 0;
		phimax = 360;
		}
		else if (strcmp(geom, "xrp") == 0)
		{
		GeometryFlag = POLAR_CAP;

		if (!phimax \|\| phimax == 0x0)
		if (latmin == 0)
		{
		phimax = 360;
		printf(
		"For geometry with polarcap please select also the minimum spot latitude with "
		"latmin=latmin\n");
		return (1);
		}

		if (!spin)
		spin = 500;
		if (timefunct == NULL)
		{
		printf(
		"For emission from a spot, select the case for a single pulse at given phase or results "
		"from the whole phase\n");
		printf("In the first case use timefunct=single, otherwise timefunct=all\n");
		return (1);
		}

		/=====================================================================/
		/Define the boundaries of integration/
		/=====================================================================/
		else if (strcmp(timefunct, "all") == 0)
		{
		if (npt == 0)
		{
		printf(
		"Please select the number of points for phase for which computing results with "
		"npt=npt\n");
		return (1);
		}
		FlagMulti = TRUE;
		}
		else if (strcmp(timefunct, "single") == 0)
		{
		if (!phimin \|\| phimin == 0x0)
		{
		printf("\nPlease select minimum azimuthal angle for the spot with phimin=phimin\n\n");
		}

		latmin_rad = latmin / 180 * PI;
		latmax_rad = latmax / 180 * PI;
		if (!phimax \|\| phimax == 0x0)
		{
		printf("\nPlease select minimum azimuthal angle for the spot with phimax=phimax\n\n");
		}
		}
		else
		{
		printf("Unknown value for parameter timefunct\n");
		return (1);
		}
		}
		else
		{
		printf("Unknown value for parameter geom\n");
		exit(1);
		}

		xmin[0] = PI / 2 - latmax_rad;
		xmax[0] = PI / 2 - latmin_rad;
		/=======================================================================================/
		if (spin == 0)
		{
		spin = 500;
		}

		xmin[1] = phimin / 180 * PI;
		xmax[1] = phimax / 180 * PI;
		if (nsrad == 0)
		{
		nsrad = 1e6;
		}

		if (compactness == 0)
		{
		compactness = 1 / 2.5;
		}

		/=====================================================================/
		if (simulfile == NULL \|\| simulfile == 0x0)
		@@ -108,15 +221,23 @@ int main(int argc, char* argv[])
		{
		FlagFile = TRUE;
		read_file();

		}

		/=====================================================================/
		/Define the boundaries of integration/
		/=====================================================================/

		latmin_rad = latmin / 180 * PI;
		latmax_rad = latmax / 180 * PI;

		xmin[0] = PI / 2 - latmax_rad;
		xmax[0] = PI / 2 - latmin_rad;

		/=====================================================================/
		/Here call the function which performs the integration/
		/=====================================================================/

		i_obs = i_obs / 180 * PI;
		compactness = 1 / 2.5;
		E_obs = 10;
		ktbb = 3;

		@@ -127,38 +248,187 @@ int main(int argc, char* argv[])
		function_params[4] = ktbb;

		/=====================================================================/
		/Computes components of the Stokes vector/

		/Intensity/
		int nfold = 4;

		if (FlagMulti == TRUE)
		{
		ncall = nfold * npt;
		phi_step = 360. * nfold / ncall;
		}
		else
		{
		ncall = 1;
		}

		/=====================================================================/
		/Allocate space for structures/
		/=====================================================================/

		polar_cap* primary_cap = malloc(sizeof(polar_cap));

		primary_cap->Ival = (double)malloc(sizeof(double) ncall);
		primary_cap->Qval = (double)malloc(sizeof(double) ncall);
		primary_cap->Uval = (double)malloc(sizeof(double) ncall);
		primary_cap->PD = (double)malloc(sizeof(double) ncall);
		primary_cap->PA = (double)malloc(sizeof(double) ncall);
		primary_cap->delta_phase = (double)malloc(sizeof(double) ncall);
		primary_cap->nval = ncall;

		polar_cap* secundary_cap = malloc(sizeof(polar_cap));

		secundary_cap->Ival = (double)malloc(sizeof(double) ncall);
		secundary_cap->Qval = (double)malloc(sizeof(double) ncall);
		secundary_cap->Uval = (double)malloc(sizeof(double) ncall);
		secundary_cap->PD = (double)malloc(sizeof(double) ncall);
		secundary_cap->PA = (double)malloc(sizeof(double) ncall);
		secundary_cap->delta_phase = (double)malloc(sizeof(double) ncall);
		secundary_cap->nval = ncall;

		/=====================================================================/

		function_params[5] = 1;
		for (ii = 0; ii < ncall; ii++)
		{
		if (FlagMulti == TRUE)
		{
		phimin = ii * phi_step;
		phimax = (ii + 1) * phi_step;
		}

		hcubature(1, observed_flux, function_params, 2, xmin, xmax, 2500, 0, 1e-40, ERROR_INDIVIDUAL,
		&stokes1, &err_stokes1);
		xmin[1] = phimin / 180 * PI;
		xmax[1] = phimax / 180 * PI;

		/====================================================/
		/Q/
		/====================================================/
		// printf("%d %lf %lf\n", ii, phimin, phimax);

		function_params[5] = 2;
		/=====================================================================/
		/Computes components of the Stokes vector for primary cap/
		/=====================================================================/

		hcubature(1, observed_flux, function_params, 2, xmin, xmax, 2500, 0, 1e-40, ERROR_INDIVIDUAL, &Qval, &Qval_err);
		primary_cap->Ival[ii] = call_numerical_integration(function_params, xmin, xmax, STOKES_I);
		primary_cap->Qval[ii] = call_numerical_integration(function_params, xmin, xmax, STOKES_Q);
		primary_cap->Uval[ii] = call_numerical_integration(function_params, xmin, xmax, STOKES_U);

		/=====================================================================/
		/U/
		/For polar cap geometry compute results also from the secondary spot/
		/=====================================================================/

		function_params[5] = 3;
		if (GeometryFlag == POLAR_CAP)
		{
		xmin_secondary[0] = PI - xmax[0];
		xmin_secondary[1] = xmin[1] + PI;

		xmax_secondary[0] = PI - xmin[0];
		xmax_secondary[1] = xmax[1] + PI;

		hcubature(1, observed_flux, function_params, 2, xmin, xmax, 2500, 0, 1e-40, ERROR_INDIVIDUAL, &Uval, &Uval_err);
		// printf("%lf %lf\n", xmin_secondary[0]/PI180, xmax_secondary[0]/PI180);

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

		if (secundary_cap->Ival[ii] == 0)
		{
		secundary_cap->Qval[ii] = 0;
		secundary_cap->Uval[ii] = 0;
		secundary_cap->PD[ii] = 0;
		secundary_cap->PA[ii] = 0;
		}
		else
		{
		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]);
		}
		}

		/=====================================================================/
		/In the Q-U plane, the angle is here defined bewteen -PI/2 and PI/2/
		/=====================================================================/

		printf("I %lf Q %5.4e U %5.4e\n", stokes1, Qval, Uval);
		primary_cap->PD[ii] =
		sqrt(pow(primary_cap->Qval[ii], 2) + pow(primary_cap->Uval[ii], 2)) / primary_cap->Ival[ii];

		primary_cap->PA[ii] = PA_location(primary_cap->Qval[ii], primary_cap->Uval[ii]);

		if (isnan(primary_cap->PD[ii]))
		primary_cap->PD[ii] = 0;

		if (isnan(primary_cap->PA[ii]))
		primary_cap->PA[ii] = 0;

		/=====================================================================/
		/Use equations (14) and (15) for delta phase/
		/=====================================================================/

		lat_center = (latmin_rad + latmax_rad) / 2;

		PD_obs = sqrt(pow(Qval, 2) + pow(Uval, 2)) / stokes1;
		deltaphase_primary =
		2 * PI * spin * nsrad / C_LIGHT * sin(i_obs) * sin(lat_center) * (1 - cos(phimin / 180 * PI));

		primary_cap->delta_phase[ii] = phimin / 360 + deltaphase_primary;

		deltaphase_secundary =
		2 * PI * spin * nsrad / C_LIGHT *
		(2 * cos(i_obs) * cos(lat_center) + sin(i_obs) * sin(lat_center) * (1 + cos(phimin / 180 * PI)));

		secundary_cap->delta_phase[ii] = phimin / 360 + deltaphase_secundary;

		/=====================================================================/

		if (GeometryFlag == POLAR_CAP && FlagMulti == TRUE)
		{
		fprintf(dat, "%5.4f %5.4e %4.3e %4.3e %5.4f %5.4e %4.3e %4.3e\n",
		primary_cap->delta_phase[ii], primary_cap->Ival[ii], primary_cap->PD[ii] * 100,
		primary_cap->PA[ii] / PI * 180, secundary_cap->delta_phase[ii],
		secundary_cap->Ival[ii], secundary_cap->PD[ii] * 100, secundary_cap->PA[ii] / PI * 180);
		}

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

		/===========================================================================/

		if (GeometryFlag == BOUNDARY_LAYER)
		{
		printf("Spreading layer properties\n");
		printf("I %5.4e PD %5.4f PA %4.3f\n", primary_cap->Ival[0], primary_cap->PD[0] * 100,
		primary_cap->PA[0] / PI * 180);
		}

		/===========================================================================/

		if (GeometryFlag == POLAR_CAP && FlagMulti == TRUE)
		{
		printf("Output results written to: %s\n", output);

		polar_cap* total_signal = add_signals(primary_cap, secundary_cap);
		}

		return 0;
		}

		/================================================================================/

		double call_numerical_integration(double* function_params, double* xmin, double* xmax, int stokes_flag)
		{
		double res = 0;
		double res_err = 0;

		function_params[5] = stokes_flag;

		hcubature(1, observed_flux, function_params, 2, xmin, xmax, 2500, 0, 1e-20, ERROR_INDIVIDUAL, &res, &res_err);

		return res;
		}

		/====================================================================/

		double PA_location(double Qval, double Uval)
		{
		double ratio_UQ;
		double PA_obs;

		ratio_UQ = Uval / Qval;

		@@ -185,9 +455,7 @@ int main(int argc, char* argv[])
		}
		}

		printf("\nPD %5.3f PA %5.3f \n", PD_obs * 100, PA_obs / PI * 180);

		return 0;
		return PA_obs;
		}

		/====================================================================/
		@@ -398,6 +666,11 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do
		retval = flag value * PD * sin(2 * PA);
		}

		if (isnan(PA) \|\| isnan(PD))
		{
		printf("Warning: PD %5.4f PA %5.4f\n", PD, PA);
		}

		return EXIT_SUCCESS;
		}

restframe_quantities.c

+4 −6

Original line number	Diff line number	Diff line
		@@ -40,11 +40,9 @@ double PolarizationRestFrame(double E_rest, double CosAlphaPrime)
		Pdeg = P_interp;
		}


		/To define the observed PA following the definition of Poutanen (2020) /
		/* change sign to the polarization of the slab*/


		Pdeg = -Pdeg;

		return Pdeg;