diff --git a/compila.txt b/compila.txt
index 60a214f21878b92c5a6f89285e306dbab0d35fb1..054647896d2a60e6134cee346ab60bfd0eeaa566 100755
--- a/compila.txt
+++ b/compila.txt
@@ -1,2 +1,3 @@
-$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
+$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
+
diff --git a/metropolis_gibbs.c b/metropolis_gibbs.c
index 5f2f731311fd47bf0f5ec1fda10d9924a1cb2c72..5891ad07d510754395a7f32f999069e69a3b4ed4 100644
--- a/metropolis_gibbs.c
+++ b/metropolis_gibbs.c
@@ -3,129 +3,120 @@
#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.1 // Valore iniziale di x (deve essere positivo)
-#define INITIAL_Y 0.1 // Valore iniziale di y (deve essere positivo)
+#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.01
+#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 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]);
+void metropolis_gibbs_sampling(double* x, double* y, int Nburns);
-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;
+ double FlagValue;
- 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)
+ 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)
{
- return 0;
+ FlagValue=1;
}
- else
+ else
{
-
- z_new = z_function(x_new, y_new);
- z_old = z_function(x_old, y_old);
+ FlagValue=0;
+ }
+
+ return FlagValue;
- if (z_new > z_old || clock_random() < z_new / z_old)
-
- return 1.0;
- else
- return 0;
- }
-
}
+/*====================================================================================*/
+void metropolis_gibbs_sampling(double* x_sample, double* y_sample, int Nburns) {
-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
+ int ii;
+ double x = INITIAL_X;
+ double y = INITIAL_Y;
- for (int i = 0; i < NUM_SAMPLES; i++)
+ double x_new;
+ double y_new;
+ double acceptance_x;
+ double acceptance_y;
+
+ for (ii = 0; ii <= Nburns; ii++)
{
+ /*===================================================*/
// 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
+ }
+ while (x_new < MIN_X || x_new > MAX_X);
- double acceptance_x = acceptance_ratio(x_new, x, y, y);
+ acceptance_x = acceptance_ratio(x_new, x, y, y);
- if (clock_random() < acceptance_x)
+ if (acceptance_x==1)
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
+ }
+ while (y_new < MIN_Y || y_new > MAX_Y);
- double acceptance_y = acceptance_ratio(x, x_new, y_new, y);
+ acceptance_y = acceptance_ratio(x, x_new, y_new, y);
- if (clock_random() < acceptance_y)
+ if (acceptance_y==1)
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));
-
+ /*===================================================*/
+
+ if (ii == Nburns)
+ {
+ *x_sample = x;
+ *y_sample = y;
}
}
}
-int main() {
- double samples[NUM_SAMPLES][2];
+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));
+ }
- // 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;
}