/**************************************************************/
// Compilare con:
//
// gcc -Wall -lgsl -lgslcblas -lm -o best_fit best_fit.c
//
/**************************************************************/


#include <stdio.h>
#include <string.h>
#include <gsl/gsl_fit.h>

int main (void) {

  // la x=Vm e' la tensione del Multimetro V
  // e la y=Vo e' la tensione dell'oscilloscopio (su cui abbiamo l'incertezza y_err)
  
  int i;

  double Vo[20] = {-4.10,
-3.80,
-3.60,
-3.40,
-3.20,
-3.05,
-2.85,
-2.65,
-2.45,
-2.20,
-2.00,
-1.80,
-1.60,
-1.40,
-1.20,
-1.00,
-0.80,
-0.60,
-0.40,
-0.20};
	
  double Vm[20] = {-4.0,
-3.8,
-3.6,
-3.4,
-3.2,
-3.0,
-2.8,
-2.6,
-2.4,
-2.2,
-2.0,
-1.8,
-1.6,
-1.4,
-1.2,
-1.0,
-0.8,
-0.6,
-0.4,
-0.2};

  double So[20] = {0.1,
0.1,
0.1,
0.1,
0.1,
0.05,
0.05,
0.05,
0.05,
0.05,
0.05,
0.05,
0.05,
0.05,
0.05,
0.05,
0.05,
0.03,
0.01,
0.005};
double Wo[20];
for (i=0; i<20; i++) {
	Wo[i]=(1/pow(So[i], 2));
	}


double Ic100[4] = {-19.40,
-19.00,
-18.50,
-17.50
};

double Vc100[4] = {-3.7,
-2.964,
-1.976,
-0.988,
};

double Sc100[4] = {0.1,
0.1,
0.1,
0.02
};
double Wc100[4];

for (i=0; i<4; i++) {
	Wc100[i]=(1/pow(Sc100[i], 2));
	}

double Ic50[9] = {-5.70,
-5.66,
-5.55,
-5.48,
-5.40,
-5.36,
-5.32,
-5.29,
-5.27
};

double Vc50[9] = {-3.7,
-2.964,
-1.976,
-1.482,
-0.988,
-0.790,
-0.593,
-0.494,
-0.395
};

double Sc50[9] = {0.1,
0.1,
0.1,
0.05,
0.02,
0.02,
0.02,
0.01,
0.01
};
double Wc50[9];

for (i=0; i<9; i++) {
	Wc50[i]=(1/pow(Sc50[i], 2));
	}


	
  double c0, c1, cov00, cov01, cov11, chisq, sigma0, sigma1;
  
  gsl_fit_wlinear (Vm, 1, Wo, 1, Vo, 1, 20, 
                   &c0, &c1, &cov00, &cov01, &cov11, 
                   &chisq);
  sigma0 = sqrt(cov00);
  sigma1 = sqrt(cov11);
  printf ("# 20 campioni\n");
  printf ("# Retta di best fit: Vm = (%g +/- %g) Vo + (%g +/- %g)\n", 1/c1, sigma1/pow(c1,2), -c0/c1, sqrt(pow(sigma0/c1, 2)+pow(c0*sigma1/pow(c1,2),2)));
  printf ("# sigma0 = %g\n", sigma0);
  printf ("# sigma1 = %g\n", sigma1);
  printf ("# Matrice di covarianza:\n");
  printf ("# [ %g, %g\n#   %g, %g]\n", cov00, cov01, cov01, cov11);
  printf ("# chisq = %g\n", chisq);
  printf ("# Retta per GNUplot: %g + %g*x\n\n\n", 1/c1, -c0/c1);

  gsl_fit_wlinear (Ic100, 1, Wc100, 1, Vc100, 1, 4, 
                   &c0, &c1, &cov00, &cov01, &cov11, 
                   &chisq);
  sigma0 = sqrt(cov00);
  sigma1 = sqrt(cov11);
  printf ("# 4 campioni\n");
  printf ("# Retta di best fit: Ic100 = (%g +/- %g)Vc100 + (%g +/- %g)\n", 1/c1, sigma1/pow(c1,2), -c0/c1, sqrt(pow(sigma0/c1, 2)+pow(c0*sigma1/pow(c1,2),2)));
  printf ("# sigma0 = %g\n", sigma0);
  printf ("# sigma1 = %g\n", sigma1);
  printf ("# Matrice di covarianza:\n");
  printf ("# [ %g, %g\n#   %g, %g]\n", cov00, cov01, cov01, cov11);
  printf ("# chisq = %g\n", chisq);
  printf ("# Retta per GNUplot: %g + %g*x\n\n\n", 1/c1, -c0/c1);


  gsl_fit_wlinear (Ic50, 1, Wc50, 1, Vc50, 1, 9, 
                   &c0, &c1, &cov00, &cov01, &cov11, 
                   &chisq);
  sigma0 = sqrt(cov00);
  sigma1 = sqrt(cov11);
  printf ("# 9 campioni\n");
  printf ("# Retta di best fit: Ic50 = (%g +/- %g)Vc50 + (%g +/- %g)\n", 1/c1, sigma1/pow(c1,2), -c0/c1, sqrt(pow(sigma0/c1, 2)+pow(c0*sigma1/pow(c1,2),2)));
  printf ("# sigma0 = %g\n", sigma0);
  printf ("# sigma1 = %g\n", sigma1);
  printf ("# Matrice di covarianza:\n");
  printf ("# [ %g, %g\n#   %g, %g]\n", cov00, cov01, cov01, cov11);
  printf ("# chisq = %g\n", chisq);
  printf ("# Retta per GNUplot: %g + %g*x\n\n\n", 1/c1, -c0/c1);

/*   for (i = 0; i < n; i++) {
    gsl_fit_linear_est (x[i], c0, c1, cov00, cov01, cov11, &yf, &yf_err);
    printf ("%g\t\t%g\t\t%g\t\t%g\n", y[i], x[i], yf_err, Voerr[i] );  
  }
 */
  return 0;
}