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

float xtalkb(float ts);
float gatti(float x,float h,float xtalk,float *qq);

float chifit_gatti(float xa,float ts,float *err,float *d,int *u);

float ht={0.495};
float h;

float chifit_gatti(float xa,float ts,float *err,float *d,int *u){
int i;
 float xtalk;
// float qerr=3.8;
float qerr=2.7;
float syserr=.008;
float qq[3];
float xn,xnh,xnl;
float step;
int ifin,istep;
float n1,n2,n3;
float chi2last,chi2min,chi2,chi2t,chi2l,chi2h,chi20;
float chi2t1,chi2t2,chi2t3;
float x,xl,xh,x0;
float errl,errh;
float erl,erh;
  h=ht/(.970+xa*(1.58-.970)/65.);
  ifin=0;
  istep=0;
  chi2min=1.0e12;
  chi2last=1.0e12;
  x=-0.5;
  step=0.5;
NEXT:
  istep=istep+1;
  if(istep>10000)goto FIN;
  chi2=0.0;
  if(u[0]==1){
    xtalk=xtalkb(ts-50);
    gatti(x,h,xtalk,qq);
    xn=(qq[0]*d[0]+qq[1]*d[1]+qq[2]*d[2])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[0]-xn*qq[0])*(d[0]-xn*qq[0])+(d[1]-xn*qq[1])*(d[1]-xn*qq[1])+(d[2]-xn*qq[2])*(d[2]-xn*qq[2]);
    n1=d[0]+d[1]+d[2];
    chi2=chi2+chi2t/(qerr*qerr+(syserr*n1)*(syserr*n1));
  }
  if(u[1]==1){
    xtalk=xtalkb(ts);
    gatti(x,h,xtalk,qq);
    xn=(qq[0]*d[3]+qq[1]*d[4]+qq[2]*d[5])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[3]-xn*qq[0])*(d[3]-xn*qq[0])+(d[4]-xn*qq[1])*(d[4]-xn*qq[1])+(d[5]-xn*qq[2])*(d[5]-xn*qq[2]);
    n2=d[3]+d[4]+d[5];
    chi2=chi2+chi2t/(qerr*qerr+(syserr*n2)*(syserr*n2));
  }
  if(u[2]==1){
    xtalk=xtalkb(ts+50);
    gatti(x,h,xtalk,qq);
    xn=(qq[0]*d[6]+qq[1]*d[7]+qq[2]*d[8])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[6]-xn*qq[0])*(d[6]-xn*qq[0])+(d[7]-xn*qq[1])*(d[7]-xn*qq[1])+(d[8]-xn*qq[2])*(d[8]-xn*qq[2]);
    n3=d[6]+d[7]+d[8];
    chi2=chi2+chi2t/(qerr*qerr+(syserr*n3)*(syserr*n3));
  }
  // printf(" x chi2 step  %f %f %f \n",x,chi2,step);
  if(chi2<chi2last){
    chi2last=chi2;
    if(chi2min>chi2)chi2min=chi2;
    x=x+step;
    goto NEXT;
  }
  if(step<0.0001){ifin=1;goto FIN;}
  x=x-2*step;
  step=step/2.;
  chi2last=1.0e12;
  goto NEXT;
FIN:
  x=-x;
  /* now calculate error */
  x0=-x;
  errl=.001;
  errh=.001;
  for(i=0;i<3;i++){
  xl=-x-errl;
  xh=-x+errh;
  chi2l=0.0;
  if(u[0]==1){
    xtalk=xtalkb(ts-50);
    gatti(xl,h,xtalk,qq);
    xn=(qq[0]*d[0]+qq[1]*d[1]+qq[2]*d[2])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[0]-xn*qq[0])*(d[0]-xn*qq[0])+(d[1]-xn*qq[1])*(d[1]-xn*qq[1])+(d[2]-xn*qq[2])*(d[2]-xn*qq[2]);
    n1=d[0]+d[1]+d[2];
    chi2l=chi2l+chi2t/(qerr*qerr+(syserr*n1)*(syserr*n1));
  }
  if(u[1]==1){
    xtalk=xtalkb(ts);
    gatti(xl,h,xtalk,qq);
    xn=(qq[0]*d[3]+qq[1]*d[4]+qq[2]*d[5])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[3]-xn*qq[0])*(d[3]-xn*qq[0])+(d[4]-xn*qq[1])*(d[4]-xn*qq[1])+(d[5]-xn*qq[2])*(d[5]-xn*qq[2]);
    n2=d[3]+d[4]+d[5];
    chi2l=chi2l+chi2t/(qerr*qerr+(syserr*n2)*(syserr*n2));
  }
  if(u[2]==1){
    xtalk=xtalkb(ts+50);
    gatti(xl,h,xtalk,qq);
    xn=(qq[0]*d[6]+qq[1]*d[7]+qq[2]*d[8])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[6]-xn*qq[0])*(d[6]-xn*qq[0])+(d[7]-xn*qq[1])*(d[7]-xn*qq[1])+(d[8]-xn*qq[2])*(d[8]-xn*qq[2]);
    n3=d[6]+d[7]+d[8];
    chi2l=chi2l+chi2t/(qerr*qerr+(syserr*n3)*(syserr*n3));
  }
  chi2l=chi2l-chi2;
  chi2h=0.0;
    if(u[0]==1){
    xtalk=xtalkb(ts-50);
    gatti(xh,h,xtalk,qq);
    xn=(qq[0]*d[0]+qq[1]*d[1]+qq[2]*d[2])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[0]-xn*qq[0])*(d[0]-xn*qq[0])+(d[1]-xn*qq[1])*(d[1]-xn*qq[1])+(d[2]-xn*qq[2])*(d[2]-xn*qq[2]);
    n1=d[0]+d[1]+d[2];
    chi2h=chi2h+chi2t/(qerr*qerr+(syserr*n1)*(syserr*n1));
  }
  if(u[1]==1){
    xtalk=xtalkb(ts);
    gatti(xh,h,xtalk,qq);
    xn=(qq[0]*d[3]+qq[1]*d[4]+qq[2]*d[5])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[3]-xn*qq[0])*(d[3]-xn*qq[0])+(d[4]-xn*qq[1])*(d[4]-xn*qq[1])+(d[5]-xn*qq[2])*(d[5]-xn*qq[2]);
    n2=d[3]+d[4]+d[5];
    chi2h=chi2h+chi2t/(qerr*qerr+(syserr*n2)*(syserr*n2));
  }
  if(u[2]==1){
    xtalk=xtalkb(ts+50);
    gatti(xh,h,xtalk,qq);
    xn=(qq[0]*d[6]+qq[1]*d[7]+qq[2]*d[8])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[6]-xn*qq[0])*(d[6]-xn*qq[0])+(d[7]-xn*qq[1])*(d[7]-xn*qq[1])+(d[8]-xn*qq[2])*(d[8]-xn*qq[2]);
    n3=d[6]+d[7]+d[8];
    chi2h=chi2h+chi2t/(qerr*qerr+(syserr*n3)*(syserr*n3));
  }
  chi2h=chi2h-chi2;
  errl=errl*sqrt(1.0/chi2l);
  errh=errh*sqrt(1.0/chi2h);
  // printf(" i errl %d %f \n",i,errl);
  }
  err[0]=errl;
  err[1]=errh;
  err[2]=chi2;  
  /*   x0=x0-errl;
    chi20=0.0;
    if(u[0]==1){
    xtalk=xtalkb(ts-50);
    gatti(x0,h,xtalk,qq);
    xn=(qq[0]*d[0]+qq[1]*d[1]+qq[2]*d[2])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[0]-xn*qq[0])*(d[0]-xn*qq[0])+(d[1]-xn*qq[1])*(d[1]-xn*qq[1])+(d[2]-xn*qq[2])*(d[2]-xn*qq[2]);
    chi2t1=chi2t;
    n1=d[0]+d[1]+d[2];
    chi20=chi20+(chi2t/(qerr*qerr+(syserr*n1)*(syserr*n1)));
  }
  if(u[1]==1){
    xtalk=xtalkb(ts);
    gatti(x0,h,xtalk,qq);
    xn=(qq[0]*d[3]+qq[1]*d[4]+qq[2]*d[5])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[3]-xn*qq[0])*(d[3]-xn*qq[0])+(d[4]-xn*qq[1])*(d[4]-xn*qq[1])+(d[5]-xn*qq[2])*(d[5]-xn*qq[2]);
    chi2t2=chi2t;
    n2=d[3]+d[4]+d[5];
    chi20=chi20+(chi2t/(qerr*qerr+(syserr*n2)*(syserr*n2)));
  }
  if(u[2]==1){
    xtalk=xtalkb(ts+50);
    gatti(x0,h,xtalk,qq);
    xn=(qq[0]*d[6]+qq[1]*d[7]+qq[2]*d[8])/(qq[0]*qq[0]+qq[1]*qq[1]+qq[2]*qq[2]);
    chi2t=(d[6]-xn*qq[0])*(d[6]-xn*qq[0])+(d[7]-xn*qq[1])*(d[7]-xn*qq[1])+(d[8]-xn*qq[2])*(d[8]-xn*qq[2]);
    chi2t3=chi2t;
    n3=d[6]+d[7]+d[8];
    chi20=chi20+(chi2t/(qerr*qerr+(syserr*n3)*(syserr*n3)));
  }
  printf(" chi2 chi20 %f %f \n",chi2,chi20);
  fprintf(fp4," %f %f \n",chi2,chi20);
  // fprintf(fp4," %f %f %f %d %d %d \n",chi2t1,chi2t2,chi2t3,d[0]+d[1]+d[2],d[3]+d[4]+d[5],d[6]+d[7]+d[8]); */ 
  return x;
}