//Solutions to ex9.6
#include <stdio.h>
#include <math.h>
#include "comphys.h"
#include "comphys.c"
#define NDATA 2000
#define MA 4
void decay(float x, float a[], float *y, float dyda[], int na);

main()
{
  float *x,*y,*sig,**covar,**alpha,*a,chisq,ochisq,alambda;
  int *ia,ma;
//  char name[80];
  FILE *in,*init,*out;
  int i,k,ndata;

  x = vector(1,NDATA);
  y = vector(1,NDATA);
  sig = vector(1,NDATA);
  a = vector(1,MA);
  covar = matrix(1,MA,1,MA);
  alpha = matrix(1,MA,1,MA);
  ia = ivector(1,MA);

  ma = MA;
  /* Beginning estimates */

  a[1] = 11.0;
  a[2] = 0.20;
  a[3] = 3.1;
  a[4] = 0.5;

  if((in = fopen("decay.out","r")) == NULL) nrerror("Cannot find input file");

  k = 1;
  while(fscanf(in,"%f %f",&x[k],&y[k]) != EOF) {
    sig[k] = 1.0;
    k++;
  }
  ndata = k-1;
  fclose(in);
  
  if((init = fopen("init.dat","r")) == NULL) nrerror("Cannot find input file");
  fscanf(init,"%f,%f,%f,%f",&a[1],&a[2],&a[3],&a[4]);
  fclose(init);
  
  for(i=1;i<=MA;i++) ia[i] = 1;

  alambda = -1.0;
  chisq = 0.0;
  mrqmin(x,y,sig,ndata,a,ia,ma,covar,alpha,&chisq,decay,&alambda);

  if((out = fopen("output.dat","w")) == NULL) nrerror("Cannot find input file");
  
  fprintf(out,"\nchisq = %f\n",chisq);
  ochisq = chisq + 2.0;

  while(fabs((ochisq - chisq)/ochisq) > 0.001 ) {
    ochisq = chisq;
    mrqmin(x,y,sig,ndata,a,ia,ma,covar,alpha,&chisq,decay,&alambda);
  }

  alambda = 0.0;
  mrqmin(x,y,sig,ndata,a,ia,ma,covar,alpha,&chisq,decay,&alambda);
  fprintf(out,"%f %f\n%f %f\n%f %f\n%f %f\n",
               a[1],sqrt(covar[1][1])/(ndata-ma),
               a[2],sqrt(covar[2][2])/(ndata-ma),
               a[3],sqrt(covar[3][3])/(ndata-ma),
               a[4],sqrt(covar[4][4])/(ndata-ma));
  fclose(out);
}

void decay(float x, float a[], float *y, float dyda[], int na)
{
  *y = a[1]*exp(-a[2]*x)*sin(a[3]*x + a[4]);
  dyda[1] = exp(-a[2]*x)*sin(a[3]*x + a[4]);
  dyda[2] = -x*a[1]*exp(-a[2]*x)*sin(a[3]*x + a[4]);
  dyda[3] = x*a[1]*exp(-a[2]*x)*cos(a[3]*x + a[4]);
  dyda[4] = a[1]*exp(-a[2]*x)*cos(a[3]*x + a[4]);
  return;
}