Moved the routine for sampling photons from synchrotron distribution to main MC folder

Makefile

@@ -18,7 +18,7 @@ MC_SYNC                 = ${MC_INSTDIR}/../../gr_code/synchrotron
 #========================================================================
 
 
-LOCAL_DIR_OBJ=  ${LOCAL_DIR}/main_program.o   ${LOCAL_DIR}/check_results.o  ${LOCAL_DIR}/compute_results.o ${LOCAL_DIR}/function_integrals.o ${LOCAL_DIR}/legendre_integration.o ${LOCAL_DIR}/set_deltatau.o ${LOCAL_DIR}/slab_polarization.o ${LOCAL_DIR}/solve_rte.o ${LOCAL_DIR}/mc_slab.o ${LOCAL_DIR}/compute_stokes.o ${LOCAL_DIR}/collect_photons.o  
+LOCAL_DIR_OBJ=  ${LOCAL_DIR}/main_program.o   ${LOCAL_DIR}/check_results.o  ${LOCAL_DIR}/compute_results.o ${LOCAL_DIR}/function_integrals.o ${LOCAL_DIR}/legendre_integration.o ${LOCAL_DIR}/set_deltatau.o ${LOCAL_DIR}/slab_polarization.o ${LOCAL_DIR}/solve_rte.o ${LOCAL_DIR}/mc_slab.o ${LOCAL_DIR}/compute_stokes.o ${LOCAL_DIR}/collect_photons.o  ${LOCAL_DIR}/simul_sync.o
 
 MC_POLARIZATION_OBJ=${MC_POLARIZATION}/ortonormal_axes.o ${MC_POLARIZATION}/unpolarized_seedphot.o ${MC_POLARIZATION}/polarization_probability.o ${MC_POLARIZATION}/depolarization_probability.o ${MC_POLARIZATION}/electfield_as_new.o ${MC_POLARIZATION}/sample_azim_angle.o ${MC_POLARIZATION}/allocate_stokes_pointer.o ${MC_POLARIZATION}/disk_polarizationelectricfield.o ${MC_POLARIZATION}/pdeg_atmosphere.o ${MC_POLARIZATION}/write_polardata.o ${MC_POLARIZATION}/select_quadrant.o $ ${MC_POLARIZATION}/detector_axis.o
 
@@ -26,13 +26,13 @@ MC_ELE_DISTRIBUTION_OBJ=${MC_ELE_DISTRIBUTION}/hybrid/eledist_hybrid.o ${MC_ELE_
 
 MC_COORD_OBJ=${MC_COORD}/versor_boost.o ${MC_COORD}/elmag_lorentztransform.o 
 
-MC_GENERAL_OBJ=${MC_INSTDIR}/../../gr_code/clock_random.o ${MC_INSTDIR}/../../gr_code/dmatrix.o ${MC_INSTDIR}/../../gr_code/round_to_digits.o ${MC_INSTDIR}/general_routines/gauleg.o ${MC_INSTDIR}/general_routines/numpar.o ${MC_INSTDIR}/general_routines/stringpar.o ${MC_INSTDIR}/general_routines/strrev.o ${MC_INSTDIR}/general_routines/stringconcat.o ${MC_INSTDIR}/general_routines/interpolate_kn.o ${MC_INSTDIR}/../../gr_code/dot_prod.o  ${MC_GENERAL}/read_array.o ${MC_INSTDIR}/../../gr_code/vect_prod.o ${MC_GENERAL}/read_matrix.o ${MC_GENERAL}/energy_lorentz_boost.o ${MC_GENERAL}/new_energy_erf.o ${MC_GENERAL}/make_array_angles.o ${MC_INSTDIR}/../../gr_code/init_mpi.o ${MC_INSTDIR}/../../gr_code/abs_val.o ${MC_GENERAL}/make_array_energy.o
+MC_GENERAL_OBJ=${MC_INSTDIR}/../../gr_code/clock_random.o ${MC_INSTDIR}/../../gr_code/dmatrix.o ${MC_INSTDIR}/../../gr_code/round_to_digits.o ${MC_INSTDIR}/general_routines/gauleg.o ${MC_INSTDIR}/general_routines/numpar.o ${MC_INSTDIR}/general_routines/stringpar.o ${MC_INSTDIR}/general_routines/strrev.o ${MC_INSTDIR}/general_routines/stringconcat.o ${MC_INSTDIR}/general_routines/interpolate_kn.o ${MC_INSTDIR}/../../gr_code/dot_prod.o  ${MC_GENERAL}/read_array.o ${MC_INSTDIR}/../../gr_code/vect_prod.o ${MC_GENERAL}/read_matrix.o ${MC_GENERAL}/energy_lorentz_boost.o ${MC_GENERAL}/new_energy_erf.o ${MC_GENERAL}/make_array_angles.o ${MC_INSTDIR}/../../gr_code/init_mpi.o ${MC_INSTDIR}/../../gr_code/abs_val.o ${MC_GENERAL}/make_array_energy.o  ${MC_INSTDIR}/../../gr_code/trig_funct.o
 
 MC_SEED_OBJ=${MC_SEED}/bb_spec.o
 
 MC_HYDRO_OBJ=${MC_HYDRO}/reverse_array.o
 
-MC_SYNC_OBJ=${MC_SYNC}/envelop_gaussian.o ${MC_SYNC}/optimize_envelop.o ${MC_SYNC}/sync_distribution.o ${MC_SYNC}/setup_syncvariables.o
+MC_SYNC_OBJ=${MC_SYNC}/envelop_gaussian.o ${MC_SYNC}/optimize_envelop.o ${MC_SYNC}/sync_distribution.o ${MC_SYNC}/setup_syncvariables.o ${MC_SYNC}/sync_seedphotons.o ${MC_SYNC}/metropolis_gibbs.o
 
 MC_SLAB_OBJ=${MC_SLAB}/electricfield_fromfile.o ${MC_SLAB}/read_resultsfile.o  ${MC_SLAB}/setup_gsl_objects.o ${MC_SLAB}/compute_cumulative_angdist.o 
 

compila.txt

 
-$CC -o zio metropolis_gibbs.c  ${MC_INSTDIR}/../../gr_code/clock_random.c ${MC_INSTDIR}/../../gr_code/synchrotron/sync_distribution.c -I ${HEADERS_GSL}   -I ${MC_INSTDIR}/../../include/ -I ${HEADERS_MPI} -I ${HEADERS_CFITSIO}  -L ${LIB_GSL}  -lm -lgsl -lgslcblas
+$CC -o zio main_test_sync.c  ${MC_INSTDIR}/../../gr_code/synchrotron/metropolis_gibbs.c  ${MC_INSTDIR}/../../gr_code/clock_random.c ${MC_INSTDIR}/../../gr_code/synchrotron/sync_distribution.c -I ${HEADERS_GSL}   -I ${MC_INSTDIR}/../../include/ -I ${HEADERS_MPI} -I ${HEADERS_CFITSIO}  -L ${LIB_GSL}  -lm -lgsl -lgslcblas
 

main_program.c

 #include "iteration_functions.h"
+#include "globals.h"
 #include <mpi.h>
 #include <stdio.h>
 
 char* set_filename(char* root, char* tau_char, char* seed_char, char* albedo, char* method);
 
 int seed;
-double kte;
-double ktbb;
 double tau_c;
-double disktau;
 double albedobase;
-double gmin;
-double gmax;
-double p1;
+
 
 char* polarization_file;
 char* diffusion_file;
-char* polardisk;
 
-int FlagEleDist;
-int nph;
-int rank;
+
+//int nph;
+//int rank;
 
 int kiter;
 int DiskPolarizationFlag;
 int Albedoflag;
 
-char* outspec;
-char* polarfile;
+//char* outspec;
+//char* polarfile;
 char* integral;
-char* diffusion;
+//char* diffusion;
 char* polarinput;
 
 FILE* dat_spectrum;
@@ -289,7 +284,7 @@ int main(int argc, char* argv[])
 
     //status = slab_mc(nph, seed);
 
-    optimize_envelop();
+    simul_sync();
     
     
     // check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);

main_test_sync.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <local_header.h>
+#include <synchrotron.h>
+
+
+int main()
+{
+
+	int ii;
+	double energy=0;
+	double csi=0;
+
+	for (ii=0; ii< 150000; ii++)
+	{
+		metropolis_gibbs_sampling(&energy, &csi, 50);
+		printf("%lf %lf  %5.4e\n", energy, csi, Psync_total(energy, csi));
+	}
+
+
+    return 0;
+}

mc_slab.c

@@ -25,28 +25,22 @@ double distance_top(double z0, double kz, double H);
 /*Datatypes defined in the header globals.h*/
 /*==================================================================*/
 
-double obsmindeg;
-double obsmaxdeg;
+//double obsmindeg;
+//double obsmaxdeg;
 
-double* col1;
-double* col2;
-double* array_angles;
-double* array_ubounds;
 
-double** array_gb;
-double** ptr_gb;
-double** ptr_2colfile;
-double** array_gb;
+//double* array_angles;
+//double* array_ubounds;
 
-int nstepangles;
+//double** array_gb;
+//double** ptr_2colfile;
 
-char* matrix_gb;
 
 const gsl_root_fsolver_type* T_phi;
 gsl_root_fsolver* s_phi;
 
-const gsl_root_fsolver_type* T_disk;
-gsl_root_fsolver* s_disk;
+//const gsl_root_fsolver_type* T_disk;
+//gsl_root_fsolver* s_disk;
 
 const gsl_root_fsolver_type* T_maxw;
 gsl_root_fsolver* s_maxw;
@@ -83,7 +77,7 @@ int slab_mc(int nphot, int seed)
   double Hphot;
   double L_top;
   double tau;
-  double eps;
+  //double eps;
   double l_int;
   double csi;
   double Pdeg;

metropolis_gibbs.c

-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <time.h>
-#include <local_header.h>
-#include <synchrotron.h>
-
-#define NUM_SAMPLES 400000
-#define BURN_IN 1000  // Numero di iterazioni di "bruciatura"
-#define INITIAL_X 0.29  // Valore iniziale di x (deve essere positivo)
-#define INITIAL_Y 0  // Valore iniziale di y (deve essere positivo)
-
-#define MIN_X 0.1
-#define MAX_X 10.0  // Massimo valore di x
-
-#define MIN_Y -5 // Massimo valore di y
-#define MAX_Y 5  // Massimo valore di y
-
-
-double acceptance_ratio(double x_new, double x_old, double y_new, double y_old);
-void metropolis_gibbs_sampling(double* x, double* y, int Nburns);
-
-
-/*====================================================================================*/
-
-double acceptance_ratio(double x_new, double x_old, double y_new, double y_old) 
-{
-    double z_new;
-    double z_old;
-    double FlagValue;
-
-    z_new = Psync_total(x_new, y_new);
-    z_old = Psync_total(x_old, y_old);
-
-    if (z_new > z_old || clock_random() < z_new / z_old)
-    {
-    	FlagValue=1;
-    }
-    else
-    {
-    	FlagValue=0;
-    }
-
-    return FlagValue;
-
-
-}
-
-/*====================================================================================*/
-void metropolis_gibbs_sampling(double* x_sample, double* y_sample, int Nburns) {
-
-
-	int ii;
-    double x = INITIAL_X;
-    double y = INITIAL_Y;
-
-    double x_new;
-    double y_new;
-    double acceptance_x;
-    double acceptance_y;
-
-    for (ii = 0; ii <= Nburns; ii++)
-    {
-    	/*===================================================*/
-        // Campiona x utilizzando il Metropolis step
-    	/*===================================================*/
-        do 
-        {
-            x_new = x + (clock_random() - 0.5);  // Propone un nuovo valore di x
-        }
-        while (x_new < MIN_X || x_new > MAX_X);
-
-        acceptance_x = acceptance_ratio(x_new, x, y, y);
-
-        if (acceptance_x==1)
-            x = x_new;
-
-        /*===================================================*/
-        // Campiona y utilizzando il Metropolis step
-        /*===================================================*/
-
-        do 
-        {
-            y_new = y + (clock_random() - 0.5);  // Propose un nuovo valore di y
-        }
-        while (y_new < MIN_Y || y_new > MAX_Y);
-
-        acceptance_y = acceptance_ratio(x, x_new, y_new, y);
-
-        if (acceptance_y==1)
-            y = y_new;
-
-        /*===================================================*/
-        // Salva il campione dopo il periodo di "bruciatura"
-        /*===================================================*/
-
-        if (ii == Nburns)
-        {
-        	*x_sample = x;
-        	*y_sample = y;
-        }
-    }
-}
-
-
-
-int main()
-{
-
-	int ii;
-	double energy=0;
-	double csi=0;
-
-	for (ii=0; ii< 150000; ii++)
-	{
-		metropolis_gibbs_sampling(&energy, &csi, 50);
-		printf("%lf %lf  %5.4e\n", energy, csi, Psync_total(energy, csi));
-	}
-
-
-    return 0;
-}

simul_sync.c

+#include <local_header.h>
+#include "common_header.h"
+#include "functions.h"
+#include "globals.h"
+#include "hydrodynamics.h"
+#include "variables.h"
+
+void simul_sync()
+{
+	
+	int ii;
+	double energy=0;
+	double csi=0;
+	
+	kerr_data* ptr_data;
+	ptr_data = (kerr_data*)malloc(sizeof(kerr_data));
+	
+	stokes_parameters ptr_stokes[POLAR_NBOUNDS - 1];
+	stokes_parameters ptr_stokes_average[POLAR_NBOUNDS - 1];
+	
+	/*============================================================================*/
+	
+	 matrix_gb =
+	      stringconcat(getenv("MC_INSTDIR"), "/electron_distributions/maxwellian/gammabeta_lab.fits");
+	  array_gb = read_matrix(matrix_gb);
+
+	  ptr_data->array_gb = array_gb;
+	  if (array_gb == NULL)
+	  {
+	    printf("\nError in assigning pointer memory to pointer array_gb\n");
+	    printf("Line %d of file %s\n", __LINE__, __FILE__);
+	    exit(1);
+	  }
+
+	 ptr_gb = read_array(
+	      stringconcat(getenv("MC_INSTDIR"), "/electron_distributions/maxwellian/temp_random.txt"));
+	  ptr_data->ptr_gb = ptr_gb;
+	  
+	  if (ptr_gb == NULL)
+	  {
+	    printf("\nError in assigning pointer memory to pointer ptr_gb\n");
+	    printf("Line %d of file %s\n", __LINE__, __FILE__);
+	    exit(1);
+	  }
+	  
+	  /*============================================================================*/
+
+	// init_globals(ptr_data, ptr_stokes, ptr_stokes_average);
+	  
+	 double r0=10;
+	 double theta0=PI/2;
+	 double phi=0;
+	 double source_fact=1;
+	 double ph_ene=0;
+	 
+	 sync_seedphotons(ptr_data, r0, &theta0, &phi0, &source_fact, &ph_ene);
+	
+	
+	        for (ii=0; ii< 300; ii++)
+	        {
+	                //metropolis_gibbs_sampling(&energy, &csi, 50);
+	                //printf("%lf %lf  %5.4e\n", energy, csi, Psync_total(energy, csi));
+	                
+	                //printf("%lf  %lf\n", ParalleSyncMode, PerpendicularSyncMode);
+	                sync_seedphotons(ptr_data, r0, &theta0, &phi0, &source_fact, &ph_ene);
+	                
+	        }
+
+	
+	
+}
+

slab_polarization.c

@@ -18,7 +18,7 @@ gsl_interp_accel* yacc_1d;
 
 int Nstep_tau;
 int Nstep_mu;
-double eps;
+double epsilon;
 double u_min;
 double u_max;
 double tau_s;
@@ -53,7 +53,7 @@ void slab_polarization(int kiter, int seed_location)
   Bk_array = malloc(Nstep_tau * sizeof(double));
   Ck_array = malloc(Nstep_tau * sizeof(double));
 
-  eps = 1e-5;
+  epsilon = 1e-5;
 
   tau_s = 2 * tau0 - tau_c;
 
@@ -237,8 +237,8 @@ void slab_polarization(int kiter, int seed_location)
         }
         else
         {
-          gsl_integration_qags(&I_l, array_tau[ii], 2 * tau0, eps, eps, 1000, w, &result_Il, &error);
-          gsl_integration_qags(&I_r, array_tau[ii], 2 * tau0, eps, eps, 1000, w, &result_Ir, &error);
+          gsl_integration_qags(&I_l, array_tau[ii], 2 * tau0, epsilon, epsilon, 1000, w, &result_Il, &error);
+          gsl_integration_qags(&I_r, array_tau[ii], 2 * tau0, epsilon, epsilon, 1000, w, &result_Ir, &error);
 
           if (isnan(result_Il))
           {
@@ -423,7 +423,7 @@ void Compute_Ak_array(AngularFunctionParams* ptr_data,
 
     // printf("Integration between %lf - %lf\n", u_min, u_max);
 
-    gsl_integration_qags(&F_Ak, u_min, u_max, eps, eps, 1000, w, &result, &error);
+    gsl_integration_qags(&F_Ak, u_min, u_max, epsilon, epsilon, 1000, w, &result, &error);
     Ak_array[ii] = result;
   }
 }
@@ -465,7 +465,7 @@ void Compute_Bk_array(AngularFunctionParams* ptr_data,
     ptr_data->k = k;
     F_Bk.params = ptr_data;
 
-    gsl_integration_qags(&F_Bk, u_min, u_max, eps, eps, 1000, w, &result, &error);
+    gsl_integration_qags(&F_Bk, u_min, u_max, epsilon, epsilon, 1000, w, &result, &error);
     Bk_array[ii] = result;
   }
 }
@@ -507,7 +507,7 @@ void Compute_Ck_array(AngularFunctionParams* ptr_data,
     ptr_data->k = k;
     F_Ck.params = ptr_data;
 
-    gsl_integration_qags(&F_Ck, u_min, u_max, eps, eps, 1000, w, &result, &error);
+    gsl_integration_qags(&F_Ck, u_min, u_max, epsilon, epsilon, 1000, w, &result, &error);
     Ck_array[ii] = result;
   }
 }