#include <stdio.h>
#include <math.h>
#include <limits.h>

#define RADIUS    25
#define N          5
#define NR       125

void fit(float x[], float y[], int n, float *A, float *B, float *dA, float *dB);
float integrate(float (*f)(float), float a, float b, int n);
float density(float r);

int main() {
  float M[N][NR];
  float x[N];
  float y[N];
  int i, r;
  for (i = 0; i < N; i++) {
    int n = pow(2, i);
    for (r = 0; r < NR; r++) {
      M[i][r] = integrate(density, 0., r, n);
    }
    x[i] = 1./n;
    y[i] = M[i][NR - 1];
  }
  FILE *fp = fopen("integrals.dat", "w");
  for (r = 0; r < NR; r++) {
    fprintf(fp, "%5d ", r);
    for (i = 0; i < N; i++) {
      fprintf(fp, "%f ", M[i][r]);
    }
    fprintf(fp,"\n");
  }
  fclose(fp);
  for (i = 0; i < N; i++) {
    printf("%f %f\n", x[i], y[i]);
  }
  float A, B, dA, dB;
  fit(x, y, N, &A, &B, &dA, &dB);
  printf("# %f +- %f\n", A, dA);
}

void fit(float x[], float y[], int n, float *A, float *B, float *dA, float *dB) {
  /* esegue la regressione lineare dei dati (x,y) */
  float Stt = 0.;
  float Sty = 0.;
  float Sx  = 0.;
  float Sy  = 0.;
  float S   = (float)n; // solo per chiarezza (S non servirebbe nel caso sigma = 1)
  int i;
  for (i = 0; i < n; i++) {
    Sx += x[i];
    Sy += y[i];
  }
  for (i = 0; i < n; i++) {
    float t = x[i] - Sx/S;
    Sty += t*y[i];
    Stt += t*t;
  }
  *B = Sty/Stt;
  *A = (Sy-Sx*(*B))/S;
  *dA = sqrt((1.+Sx*Sx/(S*Stt)))/S;
  *dB = sqrt(1./Stt);
}

float integrate(float (*f)(float), float a, float b, int n) {
  /* integra, con il metodo dei rettangoli, la funzione f moltiplicata per 4pi*r^2 */
  double S = 0.;
  double x = a;
  double dx = (b-a)/n;
  int i;
  for (i = 0; i < n; i++) {
    S += 4.*M_PI*x*x*f(x)*dx;
    x += dx;
  }
  return S;
}

float density(float r) {
  // restituisce la densita' dell'ammasso in funzione di r
  return exp(-r*r/(2.*RADIUS*RADIUS));
}