Added new structure datatype for handling the case of emission from polar caps (954fc6c9) · Commits · polarization / boundary_layer

USERGUIDE.txt

+62 −27

Original line number	Diff line number	Diff line
		INSTRUCTION FOR COMPILING AND EXECUTING THE CODE


		==============================================================
		COMPILING THE CODE
		==============================================================
		@@ -8,7 +7,7 @@ COMPILING THE CODE
		The command line for creating the executable is defined in file compile.txt
		From shell of the user just type

		$chmod 755
		$chmod 755 compile.txt
		$./compile.txt

		If everything works well you find the binary executable bp_pol
		@@ -17,34 +16,47 @@ EXECUTING THE CODE

		The code is run by calling the executable and then the parameter list in format par1=value par2=value etc

		MANDATORY PARAMETERS
		==============================================================
		GEOMETRY OF THE EMISSION
		==============================================================

		latmax=maximum latitude of the boundary layer with angle computed from the equator
		1) BOUNDARY LAYER

		iobs=observer’s viewing angle
		This geometry is defined by the parameters

		geom=bl latmax=value

		==============================================================
		OPTIONAL PARAMETERS FOR THE SURFACE OF THE BOUNDARY LAYER
		==============================================================
		where latmax is latitude of the boundary layer with angle computed from the equator

		2) POLAR CAP

		The second case is emission from two polar caps and is defined by

		1) Surface area
		geom=xrp

		By default the polarization degree (PD) and polarization angle (PA) are computed by integrating over an equatorial belt
		so that the azimuthal angle ranges from 0 to 2PI and the minimum polar angle is zero.
		However, one may wish to calculate PD and PA from a spot on the neutron star surface at given orbital phase. In this case
		also the minimum latitude and the interval of phi can be defined
		The size of the two caps is set by the parameter (units of km)

		capsize=value

		latmin=maximum latitude of the boundary layer with angle computed from the equator
		phimin=minimum azimuthal angle
		phimax=maximum azimuthal angle
		while the cap latitude is given by

		thetacap=value

		2) Rest frame polarization degree and intensity as a function of the normal to the surface
		The are two possibility then: either make computations across the whole phase (XRP) or at a given
		phase only (e.g., FRB). In the first case the following parameters must be defined

		By default if no input file is provided, the local P(u) and I(u) laws are computed assuming Chandrasekhar’s profile. If one wants to
		use pre-defined P(u) and I(u) values, they must be defined in a ASCII file with three columns having the format
		timefunct=all

		while for emission only at given azimuthal angle

		timefunct=single phicap=value

		==============================================================
		Rest frame polarization degree and intensity
		as a function of the normal to the surface
		==============================================================

		By default if no input file is provided, the local P(u) and I(u) laws are computed assuming Chandrasekhar’s profile. If one wants to use pre-defined P(u) and I(u) values, they must be defined in a ASCII file with three columns having the format

		u P(u) I(u)

		@@ -52,23 +64,46 @@ Then, the file is passed from shell as

		simulfile=/path_to_file/user_file

		Note that P(u) must not be in percentage (i.e. x 100) but just as a fraction of 1, while
		the intensity I(u) must be defined as normalizaed to the value I(u=1)

		==============================================================
		WORKING EXAMPLES
		MANDATORY PARAMETERS
		==============================================================

		1) Calculate PD and PA for a BL with latitude 45 degrees and observer’s viewing angle 70 degrees
		iobs=observer’s viewing angle (degrees)

		==============================================================
		OPTIONAL PARAMETERS
		==============================================================

		NS radius in Km (default 10)

		nsrad=value

		$ ./bl_pol latmax=45 iobs=70
		Compactness, defined as Rg/Rns (default 1/2.5)

		compactness=value

		2) Calculate PD and PA for a BL with latitude 45 degrees, observer’s viewing angle 70 degrees and rest frame P(u) and I(u)
		defined in a ASCII file
		NS rotation frequency (default 500 Hz)

		$ ./bl_pol latmax=45 iobs=70 simul1=/Users/Shared/astro/slab_polarization/results/intensity_ode_tau5_seed3.dat
		spin=value

		ASCII file with results, valid for cap geometry and phase-dependent case (default name polarimetry.dat)

		output=filename

		==============================================================
		WORKING EXAMPLES
		==============================================================

		1) Calculate PD and PA for a BL with latitude 45 degrees and observer’s viewing angle 70 degrees
		with PD and angular distribution read from a file

		$ ./bl_pol latmax=45 iobs=70 geom=bl simulfile=/Users/Shared/astro/slab_polarization/results/intensity_ode_tau5_seed3.dat


		2) Polar cap(s) with size 100 m located at latitude 45 degrees with angular distribution of intensity and polarization defined in a file

		bl_pol geom=xrp timefunct=all iobs=60 capsize=0.1 thetacap=45 spin=500 timefunct=all simulfile=/astro/slab_polarization/results/intensity_ode_tau1_seed3.qdp output=newfile.dat

main_program.c

+198 −62

Original line number	Diff line number	Diff line
		@@ -7,7 +7,19 @@ double phimax;
		double i_obs;
		double spin;
		double nsrad;
		double nsrad_km;
		double compactness;
		double capsize;
		double thetacap;
		double thetacap_rad;
		double phicap;
		double phicap_rad;
		double beta;
		double gamma_val;
		double delta_phi;
		double delta_theta;
		double phase;
		double phase_obs;

		char* output;
		char* timefunct;
		@@ -89,6 +101,12 @@ int main(int argc, char* argv[])
		npt = numpar(argv[ii]);
		if (strstr(argv[ii], "output=") != NULL)
		output = stringpar(argv[ii]);
		if (strstr(argv[ii], "capsize=") != NULL)
		capsize = numpar(argv[ii]);
		if (strstr(argv[ii], "thetacap=") != NULL)
		thetacap = numpar(argv[ii]);
		if (strstr(argv[ii], "phicap=") != NULL)
		phicap = numpar(argv[ii]);
		}

		/=============================================================================================/
		@@ -113,12 +131,6 @@ int main(int argc, char* argv[])
		return 1;
		}

		if (!latmax \|\| latmax == 0x0)
		{
		printf("\nPlease select latitude of the BL with parameter latmax=latmax\n\n");
		return 1;
		}

		if (!i_obs)
		{
		printf("\nPlease select the observer's view angle with parameter iobs=iobs\n\n");
		@@ -139,19 +151,33 @@ int main(int argc, char* argv[])
		latmin = 0;
		phimin = 0;
		phimax = 360;

		if (!latmax \|\| latmax == 0x0)
		{
		printf("\nPlease select latitude of the BL with parameter latmax=latmax\n\n");
		return 1;
		}
		}
		else if (strcmp(geom, "xrp") == 0)
		{
		GeometryFlag = POLAR_CAP;

		if (latmin == 0)
		if (capsize == 0)
		{
		printf(
		"For geometry with polarcap please select also the minimum spot latitude with "
		"latmin=latmin\n");
		printf("Please select the polar caps size in km with parameter capsize=value\n");
		return (1);
		}

		if (thetacap == 0)
		{
		printf("Please select the primary cap latitude in degrees with parameter thetacap=value\n");
		return (1);
		}
		else
		{
		thetacap_rad = thetacap / 180 * PI;
		}

		if (timefunct == NULL)
		{
		printf(
		@@ -163,26 +189,16 @@ int main(int argc, char* argv[])

		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");
		}

		if (!phimax \|\| phimax == 0x0)
		if (phicap == 0)
		{
		printf("\nPlease select minimum azimuthal angle for the spot with phimax=phimax\n\n");
		printf("Please select the azimuthal angle of the primary cap with phicap=value\n");
		return (1);
		}
		phicap_rad = phicap / 180 * PI;
		}
		else
		{
		@@ -204,8 +220,14 @@ int main(int argc, char* argv[])

		if (nsrad == 0)
		{
		nsrad = 1e6;
		nsrad = NS_RADIUS;
		}
		else
		{
		nsrad = 1e5 * nsrad;
		}

		nsrad_km = nsrad / 1e5;

		if (compactness == 0)
		{
		@@ -224,12 +246,48 @@ int main(int argc, char* argv[])
		}

		/=====================================================================/
		/Define the boundaries of integration/
		/=====================================================================/
		if (GeometryFlag == POLAR_CAP)
		{
		if (thetacap < 0 \|\| thetacap >= 90)
		{
		printf("Error: latitude of the polar cap must be between 0 and 90\n");
		return (1);
		}

		beta = FunctBetaValue(spin, compactness, PI / 2 - thetacap_rad);
		gamma_val = 1 / sqrt(1 - beta * beta);

		delta_theta = capsize / (nsrad_km);
		delta_phi = capsize / (gamma_val * nsrad_km * sin(PI / 2 - thetacap_rad));

		npt = (int)ceil((2 * PI) / delta_phi);

		// printf("Number of points %d NS rad %4.3f\n", npt, nsrad_km);

		// printf("Delta_theta (rad) %lf delta_phi %lf\n", delta_theta, delta_phi);

		latmin_rad = thetacap_rad - delta_theta / 2;
		latmax_rad = thetacap_rad + delta_theta / 2;

		if (FlagMulti == FALSE)
		{
		phimin = phicap_rad - delta_phi / 2;
		phimax = phicap_rad + delta_phi / 2;
		}
		}
		else
		{
		latmin_rad = latmin / 180 * PI;
		latmax_rad = latmax / 180 * PI;

		phimin = 0;
		phimax = 2 * PI;
		}

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

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

		@@ -249,7 +307,7 @@ int main(int argc, char* argv[])

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

		int nfold = 4;
		int nfold = 2;

		if (FlagMulti == TRUE)
		{
		@@ -291,14 +349,14 @@ int main(int argc, char* argv[])
		{
		if (FlagMulti == TRUE)
		{
		phimin = ii * phi_step;
		phimax = (ii + 1) * phi_step;
		phimin = -delta_phi / 2 + ii * delta_phi;
		phimax = -delta_phi / 2 + (ii + 1) * delta_phi;
		}

		xmin[1] = phimin / 180 * PI;
		xmax[1] = phimax / 180 * PI;
		xmin[1] = phimin;
		xmax[1] = phimax;

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

		/=====================================================================/
		/Computes components of the Stokes vector for primary cap/
		@@ -320,10 +378,15 @@ int main(int argc, char* argv[])
		xmax_secondary[0] = PI - xmin[0];
		xmax_secondary[1] = xmax[1] + PI;

		// printf("%lf %lf\n", xmin_secondary[0]/PI180, xmax_secondary[0]/PI180);
		//printf("Call secundary cap\n");

		/* printf("Theta bound %lf %lf\n", xmin_secondary[0]/PI180, xmax_secondary[0]/PI180);
		printf("Phi bound %lf %lf\n\n", xmin_secondary[1]/PI180, xmax_secondary[1]/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] =
		@@ -340,10 +403,37 @@ 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];


		/if ( secundary_cap->PD[ii]100 > 15)
		{

		printf("ii %d Q %5.4e U %5.4e I %5.4e\n", ii, secundary_cap->Qval[ii], secundary_cap->Uval[ii],secundary_cap->Ival[ii]);
		printf("%lf\n", secundary_cap->Qval[ii]/secundary_cap->Qval[ii-1]);

		} */

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

		/=============================================================================/
		/Here just compute differential flux for coordinated dthetadphi*/
		/=============================================================================/

		int zio;
		double res = 0;

		function_params[5] = STOKES_I;

		zio = observed_flux(2, xmin, function_params, 1, &res);

		/*
		printf("======================================================\n");
		printf("Approximate fluxdthetadphi\n");
		printf(" %5.4e\n", res * delta_phi * delta_theta);
		printf("======================================================\n");
		*/
		/=====================================================================/
		/In the Q-U plane, the angle is here defined bewteen -PI/2 and PI/2/
		/=====================================================================/
		@@ -365,16 +455,17 @@ int main(int argc, char* argv[])

		lat_center = (latmin_rad + latmax_rad) / 2;

		deltaphase_primary =
		2 * PI * spin * nsrad / C_LIGHT * sin(i_obs) * sin(lat_center) * (1 - cos(phimin / 180 * PI));
		phase = (phimin + phimax) / 2 * 1. / (2 * PI);

		deltaphase_primary = 2 * PI * spin * nsrad / C_LIGHT * sin(i_obs) * sin(lat_center) * (1 - cos(phimin));

		primary_cap->delta_phase[ii] = phimin / 360 + deltaphase_primary;
		primary_cap->delta_phase[ii] = phase + 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)));
		(2 * cos(i_obs) * cos(lat_center) + sin(i_obs) * sin(lat_center) * (1 + cos(phimin)));

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

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

		@@ -390,11 +481,18 @@ int main(int argc, char* argv[])

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

		if (GeometryFlag == BOUNDARY_LAYER)
		if (GeometryFlag == BOUNDARY_LAYER \|\| (GeometryFlag == POLAR_CAP && FlagMulti == FALSE))
		{
		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,
		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)
		{
		printf("Secundary spot\n");
		printf("I %5.4e PD %5.4f PA %4.3e\n", secundary_cap->Ival[0], secundary_cap->PD[0] * 100,
		secundary_cap->PA[0] / PI * 180);
		}
		}

		/===========================================================================/
		@@ -418,7 +516,7 @@ double call_numerical_integration(double* function_params, double* xmin, double*

		function_params[5] = stokes_flag;

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

		return res;
		}
		@@ -474,10 +572,17 @@ double doppler_factor(double spin, double compact, double cos_alpha, double Psi,

		beta = 2 * PI * spin * sin(theta) * NS_RADIUS / (C_LIGHT * sqrt(1 - compact));

		cos_csi = -sin_alpha / sin(Psi) * sin(i_obs) * sin(phi);

		gamma = 1 / sqrt(1 - beta * beta);

		if (sin_alpha == 0)
		{
		cos_csi = 0;
		}
		else
		{
		cos_csi = -sin_alpha / sin(Psi) * sin(i_obs) * sin(phi);
		}

		delta = 1 / (gamma * (1 - beta * cos_csi));

		return delta;
		@@ -559,15 +664,15 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do
		double I_rest;
		double E_rest;
		double stokes_idx;
		double PA;
		double PD;
		double PA=0;
		double PD=0;
		double alpha;
		double beta;
		double gamma_val;

		theta = x[0];
		phi = x[1];

		// printf("theta %lf phi %lf\n", theta/PI180, phi/PI180);

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

		i_obs = ((double*)params)[0];
		@@ -588,6 +693,13 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do

		Psi = acos(CosPsi);

		if (isnan(Psi))
		{
		printf("Error: Psi is Nan\n");
		printf("Line %d\n", __LINE__);
		exit(1);
		}

		/===========================================================/
		/Call function which determines Cos(alpha)/
		/Equation (5)/
		@@ -595,6 +707,13 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do

		cos_alpha = FunctCosAlpha(u, CosPsi);

		if (isnan(cos_alpha) \|\| (cos_alpha > 1))
		{
		printf("Error: cos_alpha is %5.4f\n", cos_alpha);
		printf("Line %d\n", __LINE__);
		exit(1);
		}

		alpha = acos(cos_alpha);

		if (cos_alpha >= 0)
		@@ -613,8 +732,14 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do

		delta = doppler_factor(spin, u, cos_alpha, Psi, i_obs, theta, phi);

		/*printf("Delta %lf\n", delta);
		*/
		if (isnan(delta))
		{
		printf("Error: delta is Nan\n");
		printf("Line %d\n", __LINE__);
		printf("cos_alpha %5.4f theta %5.4f phi %5.4f\n", cos_alpha, theta / PI * 180, phi);
		exit(1);
		}

		/===========================================================/
		/Compute the lensing factor/
		/===========================================================/
		@@ -638,7 +763,7 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do

		I_rest = IntensityRestFrame(E_rest, cos_alpha_prime, ktbb);

		/printf("I_rest %4.3e\n\n", I_rest);/
		// printf("I_rest %4.3e\n\n", I_rest);

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

		@@ -647,6 +772,8 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do

		value = pow(1 - u, 1.5) * pow(delta, 4) * I_rest * cos_alpha * Lensing * sin(theta) * gamma_val;

		//printf("stokes %3.1f VALUE %5.4e\n", stokes_idx, value);

		if (stokes_idx == 1)
		{
		retval = flag value;
		@@ -666,9 +793,18 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do
		retval = flag value * PD * sin(2 * PA);
		}

		if (PD*100 > 10 && flag==1)
		{
		printf("WARNING %3.1f %lf %lf\n", stokes_idx, cos_alpha_prime, PD);


		}

		if (isnan(PA) \|\| isnan(PD))
		{
		printf("Warning: PD %5.4f PA %5.4f\n", PD, PA);
		printf("Line %d\n", __LINE__);
		// exit(1);
		}

		return EXIT_SUCCESS;

restframe_quantities.c

+1 −0

Original line number	Diff line number	Diff line
		@@ -44,6 +44,7 @@ double PolarizationRestFrame(double E_rest, double CosAlphaPrime)
		/* change sign to the polarization of the slab*/

		Pdeg = -Pdeg;
		//printf("Return for u %lf PD %lf\n", CosAlphaPrime,Pdeg);

		return Pdeg;
		}

sum_values.c

+15 −55

Original line number	Diff line number	Diff line
		@@ -4,83 +4,43 @@ polar_cap* add_signals(polar_cap* primary_cap, polar_cap* secundary_cap)
		{
		int ii;
		int count;
		double min_val;
		double max_val;

		char* file_results = "total_result.dat";
		FILE* dat;


		dat = fopen(file_results, "w");


		min_val = primary_cap->delta_phase[0];
		max_val = primary_cap->delta_phase[secundary_cap->nval - 1];

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

		printf("%lf %lf\n", min_val, max_val);
		count=secundary_cap->nval;

		/=====================================================================================/
		/Find the number of points of the primary cap phase which overlap to the secundary/
		/so that just in that interval the Stokes parameters are summed/
		/=====================================================================================/

		for (ii = 0; ii < secundary_cap->nval; ii++)
		{
		if (secundary_cap->delta_phase[ii] >= min_val && secundary_cap->delta_phase[ii] < max_val)
		{
		count++;
		}
		}

		gsl_interp_accel* acc_I = gsl_interp_accel_alloc();
		gsl_spline* spline_I = gsl_spline_alloc(gsl_interp_linear, primary_cap->nval);
		gsl_spline_init(spline_I, primary_cap->delta_phase, primary_cap->Ival, primary_cap->nval);

		gsl_interp_accel* acc_Q = gsl_interp_accel_alloc();
		gsl_spline* spline_Q = gsl_spline_alloc(gsl_interp_linear, primary_cap->nval);
		gsl_spline_init(spline_Q, primary_cap->delta_phase, primary_cap->Qval, primary_cap->nval);

		gsl_interp_accel* acc_U = gsl_interp_accel_alloc();
		gsl_spline* spline_U = gsl_spline_alloc(gsl_interp_linear, primary_cap->nval);
		gsl_spline_init(spline_U, primary_cap->delta_phase, primary_cap->Uval, primary_cap->nval);

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

		int idx = 0;
		double interp_val_I;
		double interp_val_Q;
		double interp_val_U;

		for (ii = 0; ii < secundary_cap->nval; ii++)
		{
		if (secundary_cap->delta_phase[ii] >= min_val && secundary_cap->delta_phase[ii] <= max_val)
		{
		interp_val_I = gsl_spline_eval(spline_I, secundary_cap->delta_phase[ii], acc_I);
		interp_val_Q = gsl_spline_eval(spline_Q, secundary_cap->delta_phase[ii], acc_Q);
		interp_val_U = gsl_spline_eval(spline_U, secundary_cap->delta_phase[ii], acc_U);

		total->Ival[ii] = interp_val_I + secundary_cap->Ival[ii];
		total->Qval[ii] = interp_val_Q + secundary_cap->Qval[ii];
		total->Uval[ii] = interp_val_U + secundary_cap->Uval[ii];
		total->Ival[ii] = primary_cap->Ival[ii] + secundary_cap->Ival[ii];
		total->Qval[ii] = primary_cap->Qval[ii] + secundary_cap->Qval[ii];
		total->Uval[ii] = primary_cap->Uval[ii] + secundary_cap->Uval[ii];

		total->PD[ii] = sqrt(pow(total->Qval[ii], 2) + pow(total->Uval[ii], 2)) / total->Ival[ii];

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

		fprintf(dat, "%lf %lf %lf %lf\n", secundary_cap->delta_phase[ii], total->Ival[ii],
		fprintf(dat, "%lf %lf %lf %lf\n", primary_cap->delta_phase[ii], total->Ival[ii],
		total->PD[ii] * 100, total->PA[ii] / PI * 180);


		printf("ALCE %lf %lf %lf\n", secundary_cap->delta_phase[ii], interp_val_I, total->Ival[ii]);
		}
		// printf("ALCE %lf %lf %lf\n", secundary_cap->delta_phase[ii], interp_val_I, total->Ival[ii]);
		}

		fclose(dat);