/*********************************************************************************************************


 Program: dafat -Data Analyzer For the Atom Test- version 1.0
 
 Author : Matteo Melani
 
 Date   : 22 April 1998.




*********************************************************************************************************/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dafat_10.h"




void main(void)
{
  int m,i,k,l,g,d,s,s0,s1,s2,s3,s4,s5,s6,st,sc,sn,sp,kll,cn;
   
  for (i=0; i<64; i++)
	{ 
	  x_n[i]=i;
	  sig_y[i]=0.0003;
	}
  system("clear");
  write_header_on_screen();
  printf("Enter file name to process\a :");   /* get file to process */
  fgets(line,sizeof(line),stdin);
  sscanf (line,"%s", &input_fn);  
  sscanf(input_fn,"waf%d",&wafer_id);
  printf("\tProcessing data in file %s\n",input_fn);

  
  strcat(res_fn, input_fn);                 /* make the names of the output file */
  strcat(ntu_fn, input_fn);
  strcat(sum_fn, input_fn);      

  
  input_file = fopen(input_fn ,"r");        /* open files to read and write*/
  res_file   = fopen(res_fn,"w");
  sum_file   = fopen(sum_fn,"w");
  ntu_file   = fopen(ntu_fn,"w");
  db_file    = fopen("db.txt","a");

  
  fprintf(res_file ,"\n\n\n\n******************");/* write heder on results output file */
  fprintf(res_file ,"**************************");
  fprintf(res_file ,"************************\n");
  fprintf(res_file ,"\n\twafer #%d\n\n",wafer_id);
  fprintf(res_file ,"**************************");
  fprintf(res_file ,"**************************");
  fprintf(res_file ,"****************\n\n\n\n\n");
  
  fprintf(sum_file,"LOT ID=%s\tWAFER ID=%d\n",lot_id, wafer_id); /* write heder on summary results output file */
  fprintf(sum_file,"CHIP ID\tCONTACT\tPOWER\tANALOG\t");
  fprintf(sum_file,"DACS\tDIGITAL\tCAL_INJ\tCH.MASK\t");
  fprintf(sum_file,"NOISE\tHF\tTMIN\tCLASS\tLT FLAG\n");           
  

 for(cn=0; cn<86; cn++)
    {
                               
      for (i=0; i<4; i++)             /* throw away header */ 
	  fgets(temp,sizeof(temp),input_file);
      for (i=0; i<681; i++)           /* read data from input file */     
	{                     
	  fgets(data_line, sizeof(data_line),input_file);  
	  sscanf(data_line,"%s %d %E %E %E", &data[i].test_name, &data[i].chip_id,
		 &data[i].raw_data, &data[i].spec_lo,&data[i].spec_hi);
	} 
      for (i=0; i<3; i++)             /* throw away empty lines */ 
	fgets(temp,sizeof(temp),input_file);   

      
                     
      /* contacts test */
      contacts_flag=0; 
      contacts_test=1;
      s0=0;
      for (i=0; i<25 ; i++) 
	{
	  if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	    {
	      contacts_flag=1;
	      contacts_test=0;
	      strcpy(error_contacts[s0].test_name, data[i].test_name);   
	      error_contacts[s0].raw_data = data[i].raw_data;
	      s0++;
	    } 
	}


      /* power test*/
      power_flag=0;
      power_test=1;
      s1=0;
      for (i=25; i<27 ; i++)
	{
	  if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	    {
	      power_flag=1;
	      power_test=0;
	      strcpy(error_power[s1].test_name, data[i].test_name);
	      error_power[s1].raw_data = data[i].raw_data;
	      s1++;

	    }
                     
	}
      
     /* correct some raw data */
      data[35].spec_hi=2.4;
      data[41].spec_hi=2.4;
      data[47].spec_hi=2.4;
      data[53].spec_hi=2.4;
      data[60].spec_hi=0.1;

     /* analog test */ 
      analog_flag=0;
      analog_test=1;
      s2=0;
      for(i=27;i<61;i++)
	{
	  if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	    {
	      analog_flag=1; 
	      analog_test=0;
	      strcpy(error_analog[s2].test_name, data[i].test_name);
	      error_analog[s2].raw_data = data[i].raw_data;
	      s2++;
	    }  
	} 


      /* get cal_vref and cal_base */
      cal_base=data[62].raw_data;
      cal_vref=data[63].raw_data;

      /* get cal_in, VTP and VTN for n- and p- side */                       
     
      sc=0;
      for (i=64; i<448; i=i+6)
	{    
	  cal_in_n[sc]=data[i].raw_data;
	  VTP_n[sc]=data[i+1].raw_data;
	  VTN_n[sc]=data[i+2].raw_data;
	  sc++;
	}             
      sc=0;
      for (i=448; i<577; i=i+6)
	{
	  cal_in_p[sc]=data[i].raw_data;
	  VTP_p[sc]=data[i+1].raw_data;
	  VTN_p[sc]=data[i+2].raw_data;
	  sc++;
	}
      x_p[0]=0;     
     m_VTP_n = average(VTP_n,64.0);
     m_VTP_p = average(VTP_p,22.0);
      
     /* DACs test */
     
      fit (x_n,VTN_n,63,sig_y,&a_VTN_n,&b_VTN_n,&chi2_VTN_n);
      fit (x_p,VTN_p,22,sig_y,&a_VTN_p,&b_VTN_p,&chi2_VTN_p); 
      fit (x_n,cal_in_n,63,sig_y,&a_cal_n,&b_cal_n,&chi2_cal_n); 
      fit (x_p,cal_in_p,22,sig_y,&a_cal_p,&b_cal_p,&chi2_cal_p); 
    
   

      if (     chi2_cal_p>0.1 && chi2_cal_p<30  
            && chi2_cal_n>0.1 && chi2_cal_n<30
            && chi2_VTN_p>0.1 && chi2_VTN_p<30 
            && chi2_VTN_n>0.1 && chi2_VTN_n<30  ) 
        chi2_test=100;
      else
        chi2_test=0;
       if ( chi2_test==100 
            && chi2_cal_p<10 && chi2_cal_n<10 
            && chi2_VTN_p<10 && chi2_VTN_n<10  )
        chi2_test=1;

       if (   b_VTN_n<0.12  && b_VTN_n>0.008
           && b_VTN_p<0.008 && b_VTN_p>0.004
           && b_cal_p<1.0   && b_cal_p>0.006
           && b_cal_n<1.0   && b_cal_n>0.006            
           && chi2_test>0 )
          dacs_test=chi2_test;
       else
          dacs_test=0;   
      
               
                 
     num_dig_failure=0;

      /* digial functinality tests (ATOM1to6-->ATOM42) */       
      digital_func_flag=0;
      lt_flag_digital_1a=0;
      lt_flag_digital_1b=0;
      s3=0;                
      for(i=577; i<630; i++)
	{
	  if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	    {
	      lt_flag_digital_1b++; 
	      digital_func_flag=1;
	      num_dig_failure++;
	      for (l=0; l<33; l++)  /* check if the error is due to the left token */                                 
		if (strcmp(data[i].test_name, left_token_vectors[l].lt_test_name)==0 )
		  for (g=0; g<7; g++)
		    if (data[i].raw_data==left_token_vectors[l].error_location[g])
		      {
			lt_flag_digital_1a++; 
			lt_flag_digital_1b--;
		      }     
                           
	      strcpy(error_digital_func[s3].test_name, data[i].test_name);
	      error_digital_func[s3].raw_data = data[i].raw_data;
	      s3++;
	    }
	  if (lt_flag_digital_1b)
	    digital_func_test=0;
	  else if (lt_flag_digital_1a)
	    digital_func_test=2;
	  else 
	    digital_func_test=1;
	}



      /* CAL_INJ test (ATOM44_C-->ATOM48C_40) */
      cal_inj_flag=0;
      cal_inj_test=1;
      s4=0;
      s5=0;
      counter_1o5=0;
      for(i=630; i<635; i++)
	{
	  if (data[i].raw_data > data[i].spec_lo && data[i].raw_data < data[i].spec_hi)
	    {
	      strcpy(good_cal_inj[s5].test_name, data[i].test_name); 
	      counter_1o5++;
	      s5++;
	    }                    
	  else 
	    {
	      strcpy(error_cal_inj[s4].test_name, data[i].test_name);
	      error_cal_inj[s4].raw_data = data[i].raw_data;
	      s4++;

	    }
	  if (counter_1o5 < 1)
	    { 
	      cal_inj_flag=1;
	      cal_inj_test=0;
	    }

	} 
                      


      /* channel mask test (ATOM49C_40)*/                
      chnl_mask_flag=0;
      chnl_mask_test=1;
      i=635;
      if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	{
	  chnl_mask_flag=1;
	  num_dig_failure++;
	  chnl_mask_test=0;
	}    
 


                  
      /* 45D_20 test (ATOM45D_20) */
      _45D_20_flag=0;
      _45D_20_test=1;
      i=636;
      if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	{
	  _45D_20_flag=1;
	  num_dig_failure++;
	  _45D_20_test=0;
	}
 



      /* noise test (ATOM50_40) */
      noise_flag=0;
      noise_test=1;
      i=637;
      if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	{
	  noise_flag=1;
	  if (data[i].raw_data==1411)
	    noise_test=2;
	  else
	    noise_test=0;
	}
               
                                            
      /* HF test (A7TO9HF-->ADAC37HF)*/       
      hf_flag=0;
      hf_test=1;
      s6=0;
      for(i=642; i<645; i++)
	{
	  if (data[i].raw_data < data[i].spec_lo || data[i].raw_data > data[i].spec_hi)
	    {
	      hf_flag=1;
	      num_dig_failure++;
	      hf_test=0;
	      strcpy(error_hf[s6].test_name, data[i].test_name);
	      error_hf[s6].raw_data = data[i].raw_data;
	      s6++;
	    }
	}


      /* tmin test */
      tmin_flag=0;
      tmin_test=1;
      st=0;
      for(i=646; i<681; i=i+4)
	{
	  for (l=0;l<2; l++)
	    {
	      if (data[i+l].raw_data < data[i+l].spec_lo || data[i+l].raw_data > TMIN_HI_SPEC)
		{
		  tmin_flag=1;
		  tmin_test=0;
		  strcpy(error_tmin[st].test_name, data[i+l].test_name);
		  error_tmin[st].raw_data = data[i+l].raw_data;
		  st++;
		}
	    }
	}


      /* decide chip class on the base of the tests results */   
      all_test_but_45D_20 = digital_func_test*cal_inj_test
                           *chnl_mask_test*noise_test*hf_test;
 
      if (( contacts_test*power_test*analog_test
            *digital_test_all_but4*dacs_test) 
           && _45D_20_test==0)
	   _45D_20_test=100;
 
      all_tests =  contacts_test*power_test*analog_test*dacs_test
                  *digital_func_test*cal_inj_test*chnl_mask_test
                  *_45D_20_test*noise_test*hf_test*tmin_test;
                   
 
      if (all_tests==1)
	{
	  chip_class='A';
          lt_flag=0;
	  num_A_chips++;
	  wafer_data[d].good_chip_id=data[0].chip_id;
	  d++;
	}
      else if (all_tests==200 || all_tests==400 || all_tests==20000 || all_tests==40000)
	{
	  chip_class='B';
	  lt_flag=1;
	  num_Blt_chips++;
	  wafer_data[d].pglt_chip_id=data[0].chip_id;
	  d++;
	}
      else if (all_tests==2 || all_tests==4)
	{
	  chip_class='A';
	  lt_flag=1;
	  num_Alt_chips++;
	  wafer_data[d].lt_chip_id=data[0].chip_id;
	  d++;
	}
      else if (all_tests==100 || all_tests==10000) 
	{
	  chip_class='B';
	  lt_flag=0;
	  num_B_chips++;
	  wafer_data[d].pg_chip_id=data[0].chip_id;
	  d++;
	} 
      else
	{
	  chip_class='C';
          lt_flag=0;
	  num_C_chips++;
	  wafer_data[d].bad_chip_id=data[0].chip_id;        
	  d++; 
	}




      /* print results on file */
      fprintf(res_file,"\n\n\n\n\n---------- CHIP ID:%d ---------",data[1].chip_id);
      fprintf(res_file,"--------- CLASS = %c && %d----------\n\n",chip_class, lt_flag);
 
      /* print results for contacts test */
      fprintf(res_file,"\tCONTACTS TEST = %d\n", contacts_test);
      if (contacts_flag)
	{
	  fprintf(res_file,"\t\t error detected in test:\t\tLocation:\n");
	  for (k=0; k<s0 ; k++)
	    fprintf(res_file,"\t\t%s\t\t%.4E\n",error_contacts[k].test_name,
                                          error_contacts[k].raw_data);
          fprintf(res_file,"\n");
	}


      /* print results for power test */
      fprintf(res_file,"\tPOWER TEST    = %d\n", power_test);
      if (power_flag)
	{
	  fprintf(res_file,"\t\terror detected in test:\t\tLocation:\n"); 
	  for (k=0; k<s1 ; k++)
	    fprintf(res_file,"\t\t%s\t\t%.4E\n",error_power[k].test_name,
                                          error_power[k].raw_data);
          fprintf(res_file,"\n");
	}

      /* print results for analog test */
      fprintf(res_file,"\tANALOG TEST   = %d\n", analog_test);
      if (analog_flag)
	{
	  fprintf(res_file,"\t\terror detected in test:\t\tValue:\n");
	  for (k=0; k<s2 ; k++)
	    fprintf(res_file,"\t\t%s\t\t  %.4E\n",error_analog[k].test_name,
                                            error_analog[k].raw_data);
	} 
       
      /* print results for DACSs test and fits */
       fprintf(res_file,"\tDACs TEST = %d\n", dacs_test); 
       if (f_flag)
	fprintf(res_file,"\n\tSome or all data are not in range! Fit procedure aborted!\n");
       else
	{
	 
	  fprintf(res_file,"\t    CAL_IN N-side : a = %.4f  b(slope) = %.4f\tchi2 = %.4f\n", a_cal_n, b_cal_n, chi2_cal_n);
          fprintf(res_file,"\t    CAL_IN P-side : a = %.4f  b(slope) = %.4f\tchi2 = %.4f\n", a_cal_p, b_cal_p, chi2_cal_p);	 
          fprintf(res_file,"\t    VTN N-side    : a = %.4f  b(slope) = %.4f\tchi2 = %.4f\n", a_VTN_n, b_VTN_n, chi2_VTN_n);
	  fprintf(res_file,"\t    VTN P-side    : a = %.4f  b(slope) = %.4f\tchi2 = %.4f\n", a_VTN_p, b_VTN_p, chi2_VTN_p);
	  fprintf(res_file,"\t    VTP N-side = %.4f\n",m_VTP_n);
          fprintf(res_file,"\t    VTP p-side = %.4f\n",m_VTP_p);
          fprintf(res_file,"\t    CHI2 TEST = %d\n",chi2_test);
	  fprintf(res_file,"\t    CAL_BASE = %.4f V", cal_base);
	  fprintf(res_file,"\t    CAL_VREF = %.4f V\n\n", cal_vref);
	}

      /* print results for digital functinality test */            
      fprintf(res_file,"\tDIGITAL FUNC. TEST = %d\n",digital_func_test);
      if (digital_func_flag)
	{
	  fprintf(res_file,"\t\terror detected in test:\t\tLocation:\n");
	  for (k=0; k<s3 ; k++)
	    fprintf(res_file,"\t\t%s\t\t%.0f\n",error_digital_func[k].test_name, 
                                          error_digital_func[k].raw_data);
           fprintf(res_file,"\n");
	}

      /* printf results for cal_inj test */
      fprintf(res_file,"\tCAL_INJ TEST = %d\n",cal_inj_test);
      if (cal_inj_flag)
	{
	  fprintf(res_file,"\t\terror detected in test:\t\tLocation:\n");
	  for (k=0; k<s4 ; k++)
	    fprintf(res_file,"\t\t%s\t\t%.0f\n",error_cal_inj[k].test_name,
                                          error_cal_inj[k].raw_data);
           fprintf(res_file,"\n");
	}
      else
	{
	  fprintf(res_file,"\t\ttest passed (out of 5): \n");
	  for (k=0; k<s5; k++)
	    fprintf(res_file,"\t\t%s\n",good_cal_inj[k].test_name);
           fprintf(res_file,"\n");
	}
 /* print results for 45D_20 test */
      fprintf(res_file,"\t45D_20 TEST = %d \n",_45D_20_test);
       if (_45D_20_flag)
	{
	  fprintf(res_file,"\t\terror detected in test:\t\tLocation:\n");
	  fprintf(res_file,"\t\t%s\t\t%.0f\n",data[636].test_name ,
                                                        data[636].raw_data);
        fprintf(res_file,"\n");
	}

      /* print results for HF test */
       fprintf(res_file,"\tHF TEST    = %d \n",hf_test);
      if (hf_flag)
	{
	  fprintf(res_file,"\t\terror detected in test:\t\tLocation:\n");
	  for (k=0; k<s6 ; k++)
	    fprintf(res_file,"\t\t%s\t\t%.0f\n",error_hf[k].test_name,
                                          error_hf[k].raw_data);
          fprintf(res_file,"\n");
	}
     

      /* print results for channel mask test */
      fprintf(res_file,"\tCH.MASK TEST = %d\n",chnl_mask_test);
      if (chnl_mask_flag)
        {
 	  fprintf(res_file,"\t\terror detected at location:%.0f\n",data[635].raw_data);
          fprintf(res_file,"\n");
	}
      /* print results for noise test */
      fprintf(res_file,"\tNOISE TEST   = %d\n",noise_test);
      if (noise_flag)
	{
	  fprintf(res_file,"\t\terror detected at location:%.0f\n",data[637].raw_data);
          fprintf(res_file,"\n");
	}

      /* print results for tmin test */
      fprintf(res_file,"\tTMIN TEST    = %d\n",tmin_test);
      if (tmin_flag)
	{
	  fprintf(res_file,"\t\terror detected in test:\t\tLocation:\n");
	  for (k=0; k<st ; k++)
	    fprintf(res_file,"\t\t%s\t\t%.4E\n",error_tmin[k].test_name,
                                          error_tmin[k].raw_data);
           fprintf(res_file,"\n");     
    	} 
    
 
 /* print summary results file*/
    fprintf(sum_file,"%d\t%d\t%d",data[0].chip_id,contacts_test, power_test);
    fprintf(sum_file,"\t%d\t%d\t%d",analog_test,dacs_test,digital_func_test);
    fprintf(sum_file,"\t%d\t%d",cal_inj_test*_45D_20_test,chnl_mask_test);
    fprintf(sum_file,"\t%d\t%d\t%d",noise_test,hf_test,tmin_test);
    fprintf(sum_file,"\t%c\t%d\n",chip_class,lt_flag);
     
 
      /* print results on data base file */
	  fprintf(db_file,"20886-01\t%d\t\t%d"  , wafer_id, data[0].chip_id);
	  fprintf(db_file,"\t\t%.4f\t\t%.4f\t\t%.4f" , a_VTN_n, b_VTN_n, chi2_VTN_n);
	  fprintf(db_file,"\t\t%.4f\t\t%.4f\t\t%.4f" , a_cal_n, b_cal_n, chi2_cal_n);
	  fprintf(db_file,"\t\t%.4f\t\t%.4f\t\t%.4f" , a_VTN_p, b_VTN_p, chi2_VTN_p);
	  fprintf(db_file,"\t\t%.4f\t\t%.4f\t\t%.4f" , a_cal_p, b_cal_p, chi2_cal_p);
	  fprintf(db_file,"\t\t%.4f\t\t%.4f\t\t%.4f\t\t%.4f" , m_VTP_n,m_VTP_p, cal_base, cal_vref);
	  fprintf(db_file,"\t\t%.4f\t\t%.4f\t\t%.4f" , data[38].raw_data, data[27].raw_data, data[28].raw_data);
	  fprintf(db_file,"\t\t%.4f\t\t%.4f\t\t%.4f" , data[30].raw_data, data[32].raw_data, data[34].raw_data); 
	  fprintf(db_file,"\t\t%c\t\t%d\n",chip_class,lt_flag);
	

      /* print data on ntuple file */
      fprintf(ntu_file, "\t%10.4f \t %10.4f \t %10.4f\t %10.4f",b_VTN_n,b_VTN_p ,b_cal_n, b_cal_p );
      fprintf(ntu_file, "\t%10.4f \t %10.4f \t %10.4f\t %10.4f",a_VTN_n,a_VTN_p ,a_cal_n, a_cal_p );
      fprintf(ntu_file, "\t%10.4f \t %10.4f \t %10.4f",cal_base,cal_vref,data[65].raw_data);
      fprintf(ntu_file, "\t%10.4f \t %10.4f \t %10.4f\t%c\n",data[25].raw_data,data[26].raw_data,data[34].raw_data,chip_class);
              

    }                                                                 
      printf("\tProcessing terminated!\n\a");
         

    
  /* print the wafer data on the  results on res file */
   
  fprintf(res_file,"\n\n\n\n\n\n+++++++++ SUMMARY OF RESULTS ++++++++++++++++++++++++");
  fprintf(res_file,"(%s) ++++++++++++++++\n\n\n\n\n", input_file);
  fprintf(res_file,"\tA\tA&LT\tB\tB&LT\tC\n");
  fprintf(res_file,"\n---------------------------------------------------------------\n");
  for (i=0; i<86; i++)
    fprintf(res_file,"\t%d\t%d\t%d\t%d\t%d\n", wafer_data[i].good_chip_id,wafer_data[i].lt_chip_id,
	    wafer_data[i].pg_chip_id,wafer_data[i].pglt_chip_id,wafer_data[i].bad_chip_id);
  fprintf(res_file,"\n----------------------------------------------------------------\n");
  fprintf(res_file,"\t%d \t%d \t%d \t%d \t%d\n\n\n\n",num_A_chips, num_C_chips, num_Alt_chips, num_B_chips, num_Blt_chips);
  fprintf(res_file,"NUMBER OF CHIPS TESTED = 86\n\n");
  fprintf(res_file,"NUMBER OF A    CHIPS  = %d\n\n",num_A_chips);
  fprintf(res_file,"NUMBER OF A&LT CHIPS  = %d\n\n",num_Alt_chips); 
  fprintf(res_file,"NUMBER OF B    CHIPS  = %d\n\n",num_B_chips);
  fprintf(res_file,"NUMBER OF B&LT CHIPS  = %d\n\n",num_Blt_chips);
  fprintf(res_file,"NUMBER OF C    CHIPS  = %d\n\n",num_C_chips);
            
  yield=(((float) num_A_chips)/86)*100;  
  yield_lt=(((float) num_A_chips+num_Alt_chips)/86)*100;
  fprintf(res_file,"\nYIELD = %.2f %% ",yield); 
  fprintf(res_file,"\t\t(YIELD WITH LT CHIPS= %.2f %%) \n\n\n",yield_lt);


  /* close all files before exit */

  fclose(input_file);
  fclose(res_file);
  fclose(sum_file);
  fclose(db_file );
  fclose(sum_file);
  fclose(ntu_file);

  return ;


}




/************************************************************************************** 

     Function: write_header_on_screen
  

        Usage: the function write a header to present the program on the sdtout
                           
***************************************************************************************/
void write_header_on_screen(void)
{
  printf("**********************************************************\n");
  printf("*                                                        *\n");
  printf("*        DAFAT  - Data Analyzer For the AToM Test        *\n");
  printf("*                                                        *\n");
  printf("*                                                        *\n");
  printf("*  Version: 1.0                                          *\n");
  printf("*  Author : Matteo Melani                                *\n");
  printf("*           Lawrence Berkeley National Laboratory        *\n");
  printf("*           e-mail: mmelani@lbl.gov                      *\n");
  printf("*                                                        *\n");
  printf("**********************************************************\n");
  return;
}

void fit (float x[], float y[], int ndata, float sig[], float *a, float *b, float *chi2) 
{ 
  int i;  
  float wt,t; 
  float sxoss=0.0; 
  float sx=0.0;
  float sy=0.0;
  float st2=0.0; 
  float ss;
  float sigdat; 
  float siga,sigb; 

  *b=0.0;  
  ss=0.0; 

  for (i=0; i<ndata; i++) 
    {  
      if ( fabs(y[i]) < 0.000001)
 	  y[i]= 1.0; 
 
      wt=1.0/(sig[i]*sig[i]); 
      ss += wt; 
      sx += x[i]*wt;
      sy += y[i]*wt; 
    } 
  sxoss=sx/ss;
  for (i=0; i<ndata;i++) 
    { 
      t=(x[i]-sxoss)/sig[i]; 
      st2 += t*t; 
      *b += t*y[i]/sig[i]; 
    }

  *b /= st2;                            
  *a=(sy-sx*(*b))/ss;
  siga=sqrt((1.0+sx*sx/(ss*st2))/ss); 
  sigb=sqrt(1.0/st2); 
  *chi2=0.0; 
  for (i=0;i<ndata;i++) 
    *chi2 += (((y[i]-(*a)-(*b)*x[i])/sig[i])*((y[i]-(*a)-(*b)*x[i])/sig[i])) /(ndata-2); 
  if (*chi2>98)             /* give me a reasonable chi2! */ 
    *chi2=99.0; 
  else if (*chi2<0.1)
    *chi2=0.09;   

  return ;  

} 

 

 

/****************************************************************************

     Function: average


        Usage: the function calculates the average of an array of float numbers
              on number. Given the array and the number of elements in the array
              it returns the average.

*****************************************************************************/
float average (float dati[],float num_dati)
{
  int i;
  float sum_dati=0;
  float average=0;
  
  for (i=0; i< num_dati; i++)
    sum_dati += dati[i];
  if (num_dati !=0.0)
    { 
      average=sum_dati/num_dati;
      return average;
    }
  else
    {
      printf("\nError! Dividing by 0 in function average!\n");
      return 0.0;
    }

}