Commit fdeddc5e authored by Ruben Farinelli's avatar Ruben Farinelli
Browse files

Added new routine which reads just once the file with results for plane-parallel geometry

parent 7c3a9cfe

USERGUIDE.txt

0 → 100644
+74 −0
Original line number Diff line number Diff line 
INSTRUCTION FOR  COMPILING AND EXECUTING THE CODE





==============================================================

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 

$./compile.txt



If everything works well you find the binary executable bp_pol



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



latmax=maximum latitude of the boundary layer with angle computed from the equator



iobs=observer’s viewing angle





==============================================================

OPTIONAL PARAMETERS FOR THE SURFACE OF THE BOUNDARY LAYER

==============================================================



1) Surface area



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





latmin=maximum latitude of the boundary layer with angle computed from the equator

phimin=minimum azimuthal angle

phimax=maximum azimuthal angle





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



Then, the file is passed from shell as



simulfile=/path_to_file/user_file





==============================================================

WORKING EXAMPLES

==============================================================



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



$ ./bl_pol latmax=45 iobs=70





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



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

compile.txt

+1 −1
Original line number Diff line number Diff line 
gcc -O3 -Wall -o bl_pol main_program.c hcubature.c restframe_quantities.c polarization_angle.c numpar.c stringpar.c strrev.c interpolate_trial.c trig_funct.c -lm 
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

headers.h

+14 −1
Original line number Diff line number Diff line
@@ -16,9 +16,14 @@ 
#define PI 3.14159

#define enep 2.73



extern int tot_elements;

extern int FlagFile;

extern char* poltype;

extern char* simul1;

extern char* simulfile;



extern double* Array_costheta;

extern double* Array_pdeg;

extern double* Array_Intensity;
@@ -54,3 +59,11 @@ double csc(double a); 
double arctan(double yq, double xq);

double reduce_azim(double phi);

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

void read_file(void);



void interpolate_values(double target_costheta,

                 double* Array_costheta,

                 double* Array_pdeg,

                 double* Array_Intensity,

                 double* interpolated_pdeg,

			double* interpolated_intensity);

interpolate_trial.c→interpolate_functions.c

+42 −0
Original line number Diff line number Diff line 
#include "headers.h"



#include <stdio.h>

#include <stdlib.h>

#include <string.h>



// Function for linear interpolation



void interpolate(double target_costheta,

void interpolate_values(double target_costheta,

                 double* Array_costheta,

                 double* Array_pdeg,

                 double* Array_Intensity,

                 int data_size,

                 double* interpolated_pdeg,

                 double* interpolated_intensity)

{

  int i;

  for (i = 0; i < data_size-1; 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]);



      *interpolated_pdeg = Array_pdeg[i - 1] + (target_costheta - Array_costheta[i - 1]) *
@@ -35,100 +31,12 @@ void interpolate(double target_costheta, 
  }

}



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

{

  double* Array_costheta;

  double* Array_pdeg;

  double* Array_Intensity;





  int tot_lines = 0;

  int ii;

  char line[50];



  FILE* file;



  file = fopen(simul1, "r");

  if (file == NULL)

void Calculate_InterpolatedIntensity(double target_costheta, double* interpolated_intensity, 

		double* interpolated_pdeg)

{

    printf("Error opening file '%s'. Please check the file name and location.\n", simul1);

    exit(1);

  }



  /*====================================================*/

  /*First determine the number of lines of the file*/

  /*====================================================*/



  while (!feof(file))

  {

    if (fgets(line, 50, file) == NULL)

    {

      break;

    }

    else if (strstr(line, "#") != NULL)

    {

      continue;

    }

    else if (strstr(line, "!") != NULL)

    {

      continue;

    }

    else

    {

      tot_lines++;

    }

  }



  rewind(file);





  //printf("Number of lines %d\n", tot_lines);

  Array_costheta = (double*) malloc(tot_lines * sizeof(double));

  Array_pdeg = (double*) malloc(tot_lines * sizeof(double));

  Array_Intensity = (double*)malloc(tot_lines * sizeof(double));



  if (Array_costheta == NULL || Array_pdeg == NULL || Array_Intensity == NULL)

  {

    printf("Memory allocation failed. Exiting...\n");

    free(Array_costheta);

    free(Array_pdeg);

    free(Array_Intensity);

    fclose(file);

    exit(1);

  }







  ii=0;

  while (!feof(file))

  {

	  if (fgets(line, 50, file) == NULL)

	  {

	      break;

	  }

	  else if (strstr(line, "!") != NULL)

	  {

		  continue;

	  }

	  else

	  {

		 sscanf(line, "%lg %lg %lg\n", &Array_costheta[ii], &Array_pdeg[ii], &Array_Intensity[ii]);

		 //printf("CIAO %lf %lf %lf\n", Array_costheta[ii], Array_pdeg[ii], Array_Intensity[ii]);



		  ii++;

	  }

  }



  //exit(1);





  // Close the file

  fclose(file);

	

  interpolate(target_costheta, Array_costheta, Array_pdeg, Array_Intensity, tot_lines,

	//printf("Target %lf  %d\n", target_costheta, tot_elements);

  interpolate_values(target_costheta, Array_costheta, Array_pdeg, Array_Intensity, 

              interpolated_pdeg, interpolated_intensity);

 

  free(Array_costheta);

  free(Array_pdeg);

  free(Array_Intensity);

}

main_program.c

+10 −12
Original line number Diff line number Diff line
@@ -8,7 +8,7 @@ double i_obs; 
double spin;



char* poltype;

char* simul1;

char* simulfile;

int FlagFile;



int main(int argc, char* argv[])
@@ -49,8 +49,8 @@ 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], "simul1=") != NULL)

      simul1 = stringpar(argv[ii]);

    if (strstr(argv[ii], "simulfile=") != NULL)

      simulfile = stringpar(argv[ii]);

    if (strstr(argv[ii], "poltype=") != NULL)

      poltype = stringpar(argv[ii]);

  }
@@ -96,16 +96,19 @@ int main(int argc, char* argv[]) 
  xmax[0] = PI / 2 - latmin_rad;



  xmin[1] = phimin / 180 * PI;

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

  xmax[1] = phimax / 180 * PI;

  ;



  /*=====================================================================*/

  if (simul1 == NULL || simul1 == 0x0)

  if (simulfile == NULL || simulfile == 0x0)

  {

    FlagFile = FALSE;

  }

  else

  {

    FlagFile = TRUE;

    read_file();

    

  }



  /*=====================================================================*/
@@ -183,8 +186,6 @@ int main(int argc, char* argv[]) 
  printf("\nPD %5.3f \n", PD_obs * 100);

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







  return 0;

}
@@ -383,16 +384,13 @@ int observed_flux(unsigned dim, const double* x, void* params, unsigned fdim, do 
  }

  else if (stokes_idx == 2)

  {



    PA = polarization_angle(alpha, Psi, theta, phi, i_obs, beta);

    PD = PolarizationRestFrame(E_rest, cos_alpha_prime);



    *retval = flag * value * PD * cos(2 * PA);



  }

  else

  {



    PA = polarization_angle(alpha, Psi, theta, phi, i_obs, beta);

    PD = PolarizationRestFrame(E_rest, cos_alpha_prime);