|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
|
|
#define NPTS 5000
|
|
#define NTB 2
|
|
#define NTBEG 128
|
|
#define NUMNT 13
|
|
|
|
double ptheta(double theta, double lambda);
|
|
void trova_massimo(double lambda, double *Pmax, double *thmax);
|
|
double theta_prob(double lambda, double Pmax);
|
|
void salva_energia(double Uarr[NUMNT + 1]);
|
|
|
|
int main() {
|
|
int i, n, nt;
|
|
double pmax, theta, lambda, thmax;
|
|
double Uarr[NUMNT + 1] = {0};
|
|
srand48(42);
|
|
|
|
do {
|
|
printf("Inserisci il valore di lambda (in [0.1, 10]):");
|
|
scanf("%lf", &lambda);
|
|
}
|
|
while(lambda < 0.1 || lambda > 10.0);
|
|
|
|
trova_massimo(lambda, &pmax, &thmax);
|
|
printf("Pmax = %.4f, thetamax = %.4f\n", pmax, thmax);
|
|
for(n = 0; n <= NUMNT; n++) {
|
|
nt = NTBEG * pow(NTB, n);
|
|
for(i = 0; i < nt; i++) {
|
|
theta = theta_prob(lambda, pmax);
|
|
Uarr[n] += -cos(theta);
|
|
}
|
|
Uarr[n] /= nt;
|
|
}
|
|
salva_energia(Uarr);
|
|
return 0;
|
|
}
|
|
|
|
double ptheta(double theta, double lambda) {
|
|
double W = lambda / (exp(lambda) - exp(-lambda));
|
|
return sin(theta) * exp(lambda * cos(theta)) * W;
|
|
}
|
|
|
|
void trova_massimo(double lambda, double *Pmax, double *thmax) {
|
|
int i;
|
|
*Pmax = 0;
|
|
double dth = M_PI / NPTS;
|
|
for(i = 0; i < NPTS; i++) {
|
|
double theta = dth * i;
|
|
double Ptheta = ptheta(theta, lambda);
|
|
if(Ptheta > *Pmax || i == 0) {
|
|
*Pmax = Ptheta;
|
|
*thmax = theta;
|
|
}
|
|
}
|
|
}
|
|
|
|
double theta_prob(double lambda, double Pmax) {
|
|
double x, y;
|
|
do {
|
|
x = drand48() * M_PI;
|
|
y = drand48() * Pmax;
|
|
} while (y > ptheta(x, lambda));
|
|
|
|
return x;
|
|
}
|
|
|
|
void salva_energia(double Uarr[NUMNT + 1]) {
|
|
int n, nt;
|
|
FILE *f = fopen("energia.dat", "w");
|
|
if(f == NULL) {
|
|
printf("Errore nell'apertura del file\n");
|
|
exit(1);
|
|
}
|
|
for(n = 0; n <= NUMNT; n++) {
|
|
nt = NTBEG * pow(NTB,n);
|
|
fprintf(f, "%d %.4f\n", nt, Uarr[n]);
|
|
}
|
|
fclose(f);
|
|
}
|