/**************************************************************/
// 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 scrivi (const char* elemento, double *x, double *y, double *w, int n);

int main (void) {
  // Silicio
  // la y e' la tensione in V (su cui abbiamo l'incertezza y_err)
  // e la x e' il LOGARITMO NATURALE della corrente in mA

  int i, nsi = 19, nge=21;
    double c0, c1, cov00, cov01, cov11, chisq, sigma0, sigma1;

    double Vo[17] = {800,
690,
595,
545,
495,
445,
395,
345,
298,
249,
198,
148,
99,
74,
50,
24,
8.5};
	
  double Vm[17] = {800,
701,
600,
550,
500,
451,
400,
350,
300,
250,
201,
150.3,
100.1,
74.7,
50.7,
24.4,
7.9};

  double So[17] = {10,
10,
5,
5,
5,
5,
5,
5,
3,
3,
3,
3,
3,
1,
5E-1,
1,
1};
double Wo[17];
for (i=0; i<17; i++) {
	Wo[i]=(1/pow(So[i], 2));
	}
  gsl_fit_wlinear (Vm, 1, Wo, 1, Vo, 1, 17, 
                   &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 ("# 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);

  
  double siV[19] = {0.76, 0.72, 0.68, 0.64, 0.62, 0.6, 0.58, 0.56,
  	0.54, 0.52, 0.5, 0.48, 0.46, 0.44, 0.42, 0.4, 0.38, 0.36, 0.34};
	
  double siI[19] = {2.6, 1.92, 1.203, 0.438, 0.095, -0.408, -0.664,
  	-1.143, -1.575, -1.981, -2.430, -2.847, -3.270, -3.730, -4.135,
	-4.510, -4.962, -5.298, -5.809};

  double siS[19] = {0.02, 0.02, 0.02,0.02,0.02,0.02,0.01,0.01,0.01,0.01,
  	0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.005,0.005};
	
	
	
	
  double geV[21] = {0.3, 0.29, 0.28, 0.275, 0.260, 0.250, 0.240, 0.230, 0.220, 0.210, 0.200, 0.190,
  	0.180, 0.170, 0.160, 0.150, 0.140, 0.130, 0.120, 0.110, 0.100};
	
  double geI[21] = {0.683, 0.405, 0.182, -0.002, -0.286, -0.594, -0.918, -1.207, -1.609, -1.845, -2.180,
    -2.513, -2.882, -3.270, -3.611, -3.963, -4.342, -4.710, -5.115, -5.298, -5.809};
	
  double geS[21] = {0.01, 0.005, 0.005, 0.005, 0.005,0.005,0.005,0.005,0.005,0.005,
  	0.005,0.005,0.005,0.005,0.005,0.005,0.005,0.005,0.005,0.005,0.005};

  double *x, *y, *w;
  
/*
Function: int gsl_fit_wlinear (const double * x, const size_t xstride, const double * w, const size_t wstride, const double * y, const size_t ystride, size_t n, double * c0, double * c1, double * cov00, double * cov01, double * cov11, double * chisq)

 This function computes the best-fit linear regression coefficients (c0,c1) of the model Y = c_0 + c_1 X for the weighted datasets (x, y), two vectors of length n with strides xstride and ystride. The vector w, of length n and stride wstride, specifies the weight of each datapoint. The weight is the reciprocal of the variance for each datapoint in y  
	w(i)=1/(sigma(i))^2
	 */
	
/*
 *   Silicio
 */
  	x=siI;
  	y=siV;
	w=(double *) malloc(nsi*sizeof(double));
	
	for (i=0; i<nsi; i++) {
		w[i]=(1/pow(siS[i], 2))*1000;
		y[i]=y[i]*1000;
	}
	scrivi("silicio", x, y, w, nsi);
	free(w);

/*    gsl_fit_wlinear (x, 1, w, 1, y, 1, n, 
                   &c0, &c1, &cov00, &cov01, &cov11, 
                   &chisq);
 */				   
				   
/*
 *   Germanio
 */
  	x=geI;
  	y=geV;
	w=(double *) malloc(nge*sizeof(double));
  
	for (i=0; i<nge; i++) {
		w[i]=(1/pow(geS[i], 2))*1000;
		y[i]=y[i]*1000;
	}
	scrivi("germanio", x, y, w, nge);

/*    gsl_fit_wlinear (x, 1, w, 1, y, 1, n, 
                   &c0, &c1, &cov00, &cov01, &cov11, 
                   &chisq);
 */				   
				   
  return 0;
}


int scrivi (const char* elemento, double *x, double *y, double *w, int n) {
  
  char nomeout[14], nomedat[14], nomelin[14], nomelog[14];
  double c0, c1, cov00, cov01, cov11, chisq, I0, sigma, sigma0, sigma1, yf, yf_err;
  int i;
  
  FILE *risultati, *dati, *lineari, *loga;
  
  sprintf(nomeout, "%s%s", elemento, ".out");
  sprintf(nomedat, "%s%s", elemento, ".dat");
  sprintf(nomelin, "%s%s", elemento, ".lin");
  sprintf(nomelog, "%s%s", elemento, ".log");

  risultati = fopen(nomeout, "w");
  dati = fopen(nomedat, "w");
  lineari = fopen(nomelin, "w");
  loga = fopen(nomelog, "w");

  gsl_fit_wlinear (x, 1, w, 1, y, 1, n, &c0, &c1, &cov00, &cov01, &cov11, &chisq);
  I0 = exp(-c0/c1);
  sigma0=sqrt(cov00);
  sigma1=sqrt(cov11);
  sigma = I0*sqrt(pow(sigma0/c1, 2)+pow(c0*sigma1/pow(c1,2),2));
  fprintf (risultati, "# %s su %d campioni\n", elemento, n);
  fprintf (risultati, "# Retta di best fit y=A+Bx: V = ((%g +/- %g)+ (%g +/- %g) ln(I)) [mV]\n", c0, sigma0, c1, sigma1);
  fprintf (risultati, "# A = (%g +/- %g)[mV]; I_0 = %g +/- %g [mA]\n", c0, sigma0, I0, sigma);
  fprintf (risultati, "# B = etaV_T = (%g +/- %g) [mV];\n", c1, sigma1);
  fprintf (risultati, "# Matrice di covarianza:\n");
  fprintf (risultati, "# [%g, %g\n# %g, %g]\n", cov00, cov01, cov01, cov11);
  fprintf (risultati, "# chisq = %g\n", chisq);
  fprintf (risultati, "# Caratt. per GNUplot: %g*exp(%.5f*x)\n", I0, 1/c1);
  fprintf (risultati, "# Retta per GNUplot: %g + %g*x\n", -c0/c1, 1/c1);
  fflush(risultati);
  fclose(risultati);

//  fprintf (dati, "V [mV]\t\tI [mA]\t\tstd_dev(V)\n\n");

  for (i = 0; i < n; i++) {
    gsl_fit_linear_est (x[i], c0, c1, cov00, cov01, cov11, &yf, &yf_err);
    fprintf (dati, "%g\t\t%g\t\t%g\n", y[i], exp(x[i]), yf_err);  
  }
  fflush(dati);
  fclose(dati);

  for (i = 0; i < n; i++) {
    gsl_fit_linear_est (x[i], c0, c1, cov00, cov01, cov11, &yf, &yf_err);
    fprintf (lineari, "%g\t\t%g\n", y[i], x[i]);  
  }
  fflush(lineari);
  fclose(lineari);

  for (i = 0; i < n; i++) {
    gsl_fit_linear_est (x[i], c0, c1, cov00, cov01, cov11, &yf, &yf_err);
    fprintf (loga, "%g\t\t%g\n", y[i], log10(exp(x[i])));  
  }
  fflush(loga);
  fclose(loga);
  
  return 0;

}