diff --git a/Makefile b/Makefile
index ccff84bd3111a249f94085a28871b5c68a14274e..57ff3d676843f2502b64b8a37775237598db7074 100644
--- a/Makefile
+++ b/Makefile
@@ -12,6 +12,7 @@ MC_GENERAL = ${MC_INSTDIR}/general_routines/
MC_COORD = ${MC_INSTDIR}/coordinates_systems
MC_SEED = ${MC_INSTDIR}/seed_photons/
MC_HYDRO = ${MC_INSTDIR}/../../gr_code/hydrodynamics
+MC_SYNC = ${MC_INSTDIR}/../../gr_code/synchrotron
#========================================================================
@@ -31,6 +32,7 @@ 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_SLAB_OBJ=${MC_SLAB}/electricfield_fromfile.o ${MC_SLAB}/read_resultsfile.o ${MC_SLAB}/setup_gsl_objects.o ${MC_SLAB}/compute_cumulative_angdist.o
@@ -40,12 +42,12 @@ all: start_iteration
%.o: %.c $(DEPS)
$(CC) -c -o $@ $< $(CFLAGS)
-start_iteration: ${LOCAL_DIR_OBJ} ${MC_GENERAL_OBJ} ${MC_POLARIZATION_OBJ} ${MC_ELE_DISTRIBUTION_OBJ} ${MC_COORD_OBJ} ${MC_SEED_OBJ} ${MC_HYDRO_OBJ} ${MC_SLAB_OBJ}
+start_iteration: ${LOCAL_DIR_OBJ} ${MC_GENERAL_OBJ} ${MC_POLARIZATION_OBJ} ${MC_ELE_DISTRIBUTION_OBJ} ${MC_COORD_OBJ} ${MC_SEED_OBJ} ${MC_HYDRO_OBJ} ${MC_SLAB_OBJ} ${MC_SYNC_OBJ}
${CC} $^ -o $@ ${LDFLAGS}
clean:
- rm -f main_program.o ${LOCAL_DIR_OBJ} ${MC_GENERAL_OBJ} ${MC_POLARIZATION_OBJ} ${MC_ELE_DISTRIBUTION_OBJ} ${MC_COORD_OBJ} ${MC_SEED_OBJ} ${MC_HYDRO_OBJ} ${MC_SLAB_OBJ}
+ rm -f main_program.o ${LOCAL_DIR_OBJ} ${MC_GENERAL_OBJ} ${MC_POLARIZATION_OBJ} ${MC_ELE_DISTRIBUTION_OBJ} ${MC_COORD_OBJ} ${MC_SEED_OBJ} ${MC_HYDRO_OBJ} ${MC_SLAB_OBJ} ${MC_SYNC_OBJ}
ciao:
diff --git a/compila.txt b/compila.txt
new file mode 100755
index 0000000000000000000000000000000000000000..60a214f21878b92c5a6f89285e306dbab0d35fb1
--- /dev/null
+++ b/compila.txt
@@ -0,0 +1,2 @@
+
+$CC -o zio metropolis_gibbs.c ${MC_INSTDIR}/../../gr_code/clock_random.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 f585dc189d584ccc4fea6c6e3ea16f33a8c9a188..177828a2e01303d34cea35db8fd2052c0eb62766 100644
--- a/main_program.c
+++ b/main_program.c
@@ -287,8 +287,11 @@ int main(int argc, char* argv[])
integral = set_filename("integralpolar_mc_tau", tau_char, seed_char, albedo_char, "mc");
polarfile = set_filename("energypolar_mc_tau", tau_char, seed_char, albedo_char, "mc");
- status = slab_mc(nph, seed);
+ //status = slab_mc(nph, seed);
+ optimize_envelop();
+
+
// check_results(Spline_sample_Iu, Spline_eval_pdeg, Spline_eval_limbdark);
}
diff --git a/metropolis_gibbs.c b/metropolis_gibbs.c
new file mode 100644
index 0000000000000000000000000000000000000000..5f2f731311fd47bf0f5ec1fda10d9924a1cb2c72
--- /dev/null
+++ b/metropolis_gibbs.c
@@ -0,0 +1,131 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include <local_header.h>
+
+#define NUM_SAMPLES 400000
+#define BURN_IN 1000 // Numero di iterazioni di "bruciatura"
+#define INITIAL_X 0.1 // Valore iniziale di x (deve essere positivo)
+#define INITIAL_Y 0.1 // Valore iniziale di y (deve essere positivo)
+
+#define MIN_X 0.01
+#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 z_function(double x, double y);
+double acceptance_ratio(double x_new, double x_old, double y_new, double y_old);
+void metropolis_gibbs_sampling(double samples[NUM_SAMPLES][2]);
+
+
+double z_function(double x, double y)
+{
+
+ double f1;
+ double f2;
+ double bessel_arg;
+ double value;
+
+ bessel_arg=x / 2.0 * pow(1 + y * y, 1.5);
+
+
+ f1=x * pow(1 + y * y, 2) * pow(gsl_sf_bessel_Knu(2.0/3.0, bessel_arg),2);
+ f2=x * pow(y,2) *(1 +y*y) * pow(gsl_sf_bessel_Knu(1.0/3.0, bessel_arg),2);
+
+
+ return f1+f2; // Z(x, y) = x * BesselK^(2/3)(...)
+}
+
+double acceptance_ratio(double x_new, double x_old, double y_new, double y_old)
+{
+ double z_new;
+ double z_old;
+
+ if (x_new <= MIN_X || x_new > MAX_X || x_old <= MIN_X || x_old > MAX_X
+ || y_new > MAX_Y || y_old > MAX_Y)
+ {
+ return 0;
+ }
+ else
+ {
+
+ z_new = z_function(x_new, y_new);
+ z_old = z_function(x_old, y_old);
+
+ if (z_new > z_old || clock_random() < z_new / z_old)
+
+ return 1.0;
+
+ else
+ return 0;
+ }
+
+}
+
+
+void metropolis_gibbs_sampling(double samples[NUM_SAMPLES][2]) {
+ // Inizializza il generatore di numeri casuali
+ srand(time(NULL));
+
+ // Valori iniziali
+ double x = INITIAL_X; // Inizializza con un valore positivo
+ double y = INITIAL_Y; // Inizializza con un valore positivo
+
+ for (int i = 0; i < NUM_SAMPLES; i++)
+ {
+ // Campiona x utilizzando il Metropolis step
+ double x_new;
+ do
+ {
+ x_new = x + (clock_random() - 0.5); // Propone un nuovo valore di x
+ } while (x_new < MIN_X || x_new > MAX_X); // Assicura che x sia positivo e rientri nel dominio
+
+ double acceptance_x = acceptance_ratio(x_new, x, y, y);
+
+ if (clock_random() < acceptance_x)
+ x = x_new;
+
+ // Campiona y utilizzando il Metropolis step
+ double y_new;
+ do
+ {
+ y_new = y + (clock_random() - 0.5); // Propose un nuovo valore di y
+ } while (y_new < MIN_Y || y_new > MAX_Y); // Assicura che y sia positivo e rientri nel dominio
+
+ double acceptance_y = acceptance_ratio(x, x_new, y_new, y);
+
+ if (clock_random() < acceptance_y)
+ y = y_new;
+
+ // Salva il campione dopo il periodo di "bruciatura"
+ if (i >= BURN_IN) {
+ samples[i - BURN_IN][0] = x;
+ samples[i - BURN_IN][1] = y;
+ printf("%lf %lf %5.4e\n", x, y, z_function(x, y));
+
+ }
+ }
+}
+
+
+
+int main() {
+ double samples[NUM_SAMPLES][2];
+
+ // Esegui il campionamento utilizzando Metropolis-Gibbs
+ metropolis_gibbs_sampling(samples);
+
+ // Stampa i campioni
+ //printf("Campioni dalla distribuzione Z(x, y) = x * BesselK^(2/3)(...):\n");
+
+ /* for (int i = 0; i < NUM_SAMPLES; i++)
+ {
+ printf("(%5.4e, %5.4e)\n", samples[i][0], samples[i][1]);
+ }
+*/
+ return 0;
+}