Modified the Makefile in order to compile source files in new folder...

Modified the Makefile in order to compile source files in new folder ${MC_INSTDIR}/../../gr_code/synchrotron

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:

compila.txt

+
+$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

main_program.c

@@ -287,7 +287,10 @@ 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);
   }

metropolis_gibbs.c

+#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;
+}