Giù Marcer
093faf3fe8
move all the functions to lib.c and lib.h; add this two librearies to makefile; make plot.py more easy to use for passing to one mode to the other (show or save figure); create fit_plot and create the figure 5-fit.pdf.
87 lines
1.7 KiB
C
87 lines
1.7 KiB
C
#include <gsl/gsl_fit.h>
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include "lib.h"
|
|
|
|
|
|
// Wrapper for the gsl_function structure
|
|
//
|
|
double function (double * x, size_t dim, void * params)
|
|
{
|
|
return exp(x[0]);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////////////
|
|
|
|
// Results printer.
|
|
//
|
|
void results (size_t calls, double result, double error, double chi)
|
|
{
|
|
if (calls != 0)
|
|
{
|
|
printf ("%6.0e | ", (double)calls);
|
|
}
|
|
printf ("%5f | ", result);
|
|
printf ("%5f | ", error);
|
|
if (chi != 0)
|
|
{
|
|
printf ("%5f", chi);
|
|
}
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////////////
|
|
|
|
// Perform a fit in order to compare the data with the expected function:
|
|
//
|
|
// y = a⋅x^b → ln(y) = ln(a) + b⋅ln(x) → Y = A + B⋅X
|
|
//
|
|
// with:
|
|
//
|
|
// A = ln(a) → a = e^A
|
|
// B = b → b = B
|
|
//
|
|
// using a linear regression. For b, the results is hence compared with the
|
|
// expected one:
|
|
//
|
|
// b_exp = - 0.5
|
|
//
|
|
void fit (struct bag full_bag, double* p,
|
|
double* a, double* a_err,
|
|
double* b, double* b_err)
|
|
{
|
|
// Expected value.
|
|
//
|
|
double b_exp = -0.5;
|
|
|
|
// Parse arguments.
|
|
//
|
|
double size = full_bag.size;
|
|
double* x = full_bag.pokets.x;
|
|
double* y = full_bag.pokets.y;
|
|
for (size_t i = 0; i < size; i++)
|
|
{
|
|
x[i] = log(x[i]);
|
|
y[i] = log(y[i]);
|
|
}
|
|
|
|
// Do fit.
|
|
//
|
|
double A, B, A_err, B_err, AB_cov, sum2;
|
|
gsl_fit_linear(x, 1, y, 1, size, &A, &B, &A_err, &AB_cov, &B_err, &sum2);
|
|
|
|
// Parse results.
|
|
//
|
|
A_err = sqrt(A_err);
|
|
B_err = sqrt(B_err);
|
|
*a = exp(A);
|
|
*b = B;
|
|
*a_err = *a * A_err;
|
|
*b_err = B_err;
|
|
|
|
// Check compatibility with expected values.
|
|
//
|
|
double t = fabs(*b - b_exp)/ *b_err;
|
|
*p = 1 - erf(t/sqrt(2));
|
|
}
|
|
|