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

#define NBprofpts 80

float prof[NBprofpts];
float INVRANDMAX = 1.0/RAND_MAX;
float a1,b1,a2,b2;

float ran1()
{
  float value;

  value = INVRANDMAX*rand();

  return(value);
}

float fitvalue()
{
  long i;
  float value = 0.0;
  float x;
  float tmp;

  for(i=0;i<NBprofpts;i++)
    {
      x = 1.0*i/NBprofpts;
      tmp = prof[i] - exp( a1*pow(x,b1) + a2*pow(x,b2));
      value += tmp*tmp;
    }

  return(value);
}

int main()
{
  FILE *fp;
  long i;
  float a1b,b1b,a2b,b2b;
  float tmpf;
  long tmpl;
  float value,best_value;
  float range;

  srand(time(NULL));

  fp = fopen("APLCapo_2.000000.1024.ref.prof","r");
  for(i=0;i<NBprofpts;i++)
    {
      fscanf(fp,"%ld %f\n",&tmpl,&tmpf);
      prof[i] = tmpf;
    }
  fclose(fp);

  a1b = -2.758451;
  b1b = 2.003837;
  a2b = -0.996055;
  b2b = 5.910419;

  a1 = a1b;
  b1 = b1b;
  a2 = a2b;
  b2 = b2b;

  value = fitvalue();
  printf("fit value = %.20f\n",value);
  best_value = value;
  printf("a1b = %f;\n",a1b);
  printf("b1b = %f;\n",b1b);
  printf("a2b = %f;\n",a2b);
  printf("b2b = %f;\n",b2b);
  printf("f(x) = exp(%f*x**%f + %f*x**%f)\n",a1b,b1b,a2b,b2b);

  while(1)
    {
      range = 0.0001;
      a1 = a1b + range*(0.5-ran1());
      b1 = b1b + range*(0.5-ran1());
      a2 = a2b + range*(0.5-ran1());
      b2 = b2b + range*(0.5-ran1());
      value = fitvalue();

      if(value<best_value)
	{
	  fprintf(stderr,"NEW BEST POINT %.20f - > %.20f\n",best_value,value);
	  best_value = value;
	  a1b = a1;
	  b1b = b1;
	  a2b = a2;
	  b2b = b2; 
	  printf("a1b = %f;\n",a1b);
	  printf("b1b = %f;\n",b1b);
	  printf("a2b = %f;\n",a2b);
	  printf("b2b = %f;\n",b2b);
	  printf("f(x) = exp(%f*x**%f + %f*x**%f)\n",a1b,b1b,a2b,b2b);
	}
    }
  
  return(0);
}