From 4f592803fb2ab7616fef4bbef75a275f93a2ee51 Mon Sep 17 00:00:00 2001
From: Ruben Farinelli <ruben.farinelli@inaf.it>
Date: Wed, 18 Oct 2023 09:22:16 +0200
Subject: [PATCH] Moved the routine for sampling photons from synchrotron
distribution to main MC folder
---
Makefile | 6 +--
compila.txt | 2 +-
main_program.c | 23 ++++-----
main_test_sync.c | 24 +++++++++
mc_slab.c | 24 ++++-----
metropolis_gibbs.c | 122 --------------------------------------------
simul_sync.c | 72 ++++++++++++++++++++++++++
slab_polarization.c | 14 ++---
8 files changed, 125 insertions(+), 162 deletions(-)
create mode 100644 main_test_sync.c
delete mode 100644 metropolis_gibbs.c
create mode 100644 simul_sync.c
diff --git a/Makefile b/Makefile
index 57ff3d6..94d23c8 100644
--- a/Makefile
+++ b/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
diff --git a/compila.txt b/compila.txt
index 0546478..cd35a5f 100755
--- a/compila.txt
+++ b/compila.txt
@@ -1,3 +1,3 @@
-$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
diff --git a/main_program.c b/main_program.c
index 177828a..41549a5 100644
--- a/main_program.c
+++ b/main_program.c
@@ -1,35 +1,30 @@
#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);
diff --git a/main_test_sync.c b/main_test_sync.c
new file mode 100644
index 0000000..f2b72db
--- /dev/null
+++ b/main_test_sync.c
@@ -0,0 +1,24 @@
+#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;
+}
diff --git a/mc_slab.c b/mc_slab.c
index a4936af..2a658c6 100644
--- a/mc_slab.c
+++ b/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;
diff --git a/metropolis_gibbs.c b/metropolis_gibbs.c
deleted file mode 100644
index 5891ad0..0000000
--- a/metropolis_gibbs.c
+++ /dev/null
@@ -1,122 +0,0 @@
-#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;
-}
diff --git a/simul_sync.c b/simul_sync.c
new file mode 100644
index 0000000..1f0e78c
--- /dev/null
+++ b/simul_sync.c
@@ -0,0 +1,72 @@
+#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);
+
+ }
+
+
+
+}
+
diff --git a/slab_polarization.c b/slab_polarization.c
index 160b7a3..3f01bb3 100644
--- a/slab_polarization.c
+++ b/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;
}
}
--
GitLab