|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <time.h>
|
|
#include <math.h>
|
|
|
|
#define NPTS 10000
|
|
#define NTRIALS 100000
|
|
#define ALPHAMIN 1
|
|
#define ALPHAMAX 10
|
|
|
|
double onsager(double theta, double alpha);
|
|
void findmax(double alpha, double *f_max, double *theta_max);
|
|
double onsager_prob(double alpha, double f_max);
|
|
void salva_Sarr(double *Sarr);
|
|
|
|
int main() {
|
|
srand48(time(NULL));
|
|
|
|
int i;
|
|
double f_max, theta_max;
|
|
double Sarr[ALPHAMAX - ALPHAMIN + 1] = {0};
|
|
|
|
int alpha;
|
|
for(alpha = ALPHAMIN; alpha <= ALPHAMAX; alpha++) {
|
|
findmax(alpha, &f_max, &theta_max);
|
|
if(alpha == ALPHAMAX) {
|
|
fprintf(stderr, "Valori di f_ons massima, e theta corrispondente: %.5lf %.5lf\n", f_max, theta_max);
|
|
}
|
|
|
|
for(i = 0; i < NTRIALS; i++) {
|
|
double theta = onsager_prob(alpha, f_max);
|
|
double P2_theta = (3 * cos(theta) * cos(theta) - 1.0) / 2.0;
|
|
Sarr[alpha - ALPHAMIN] += P2_theta;
|
|
}
|
|
Sarr[alpha - ALPHAMIN] /= NTRIALS;
|
|
}
|
|
|
|
salva_Sarr(Sarr);
|
|
|
|
return 0;
|
|
}
|
|
|
|
double onsager(double theta, double alpha) {
|
|
return sin(theta) * alpha * cosh(alpha * cos(theta)) / (2.0 * sinh(alpha));
|
|
}
|
|
|
|
void findmax(double alpha, double *f_max, double *theta_max) {
|
|
double dtheta = M_PI / NPTS;
|
|
double theta;
|
|
for(*f_max = 0.0, theta = 0.0; theta <= M_PI; theta += dtheta) {
|
|
double val = onsager(theta, alpha);
|
|
if(val > *f_max) {
|
|
*f_max = val;
|
|
*theta_max = theta;
|
|
}
|
|
}
|
|
}
|
|
|
|
double onsager_prob(double alpha, double f_max) {
|
|
double x, y;
|
|
do {
|
|
x = drand48() * M_PI;
|
|
y = drand48() * f_max;
|
|
} while(y > onsager(x, alpha));
|
|
|
|
return x;
|
|
}
|
|
|
|
void salva_Sarr(double *Sarr) {
|
|
FILE *output = fopen("S_alpha.dat", "w");
|
|
int alpha;
|
|
for(alpha = ALPHAMIN; alpha <= ALPHAMAX; alpha += 1) {
|
|
fprintf(output, "%d %.5lf\n", alpha, Sarr[alpha - ALPHAMIN]);
|
|
}
|
|
fclose(output);
|
|
}
|