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

#define PI 3.14159265358979323846264338328
#define NBcheckpts 31341

#define NBcscoeff 1225

double bk[200];

long CScoeff_i[NBcscoeff];
long CScoeff_j[NBcscoeff];
double CScoeff_cos[NBcscoeff];
double CScoeff_sin[NBcscoeff];

double f = 0.9;
double OAD = 0.085;


double xsample[31341];
double ysample[NBcheckpts];
double f2sample[NBcheckpts];
double OAP2sample[NBcheckpts];
double OAT2sample[NBcheckpts];
double M2sample[NBcheckpts];


int main()
{
  FILE *fp;
  long i,k;
  double f2,OAP2,OAT2;

  double x,y,r,X,Y,R;
  double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
  long tmpl,tmpl1,tmpl2;

  // load bk
  fp = fopen("../PIAAM2_cosfitcoeff_200.dat","r");
  for(i=0;i<200;i++)
    {
      fscanf(fp,"%ld %lf\n",&tmpl,&tmp1);
      bk[tmpl] = tmp1;
    }
  fclose(fp);
  
  // load CScoeff
  if((fp = fopen("../PIAAM2_OAT25.dat","r"))==NULL)
    {
      printf("ERROR: could not read file\n");
      exit(0);
    }

  for(i=0;i<NBcscoeff;i++)
    {
      fscanf(fp,"%ld %ld %lf %lf\n",&tmpl1,&tmpl2,&tmp1,&tmp2);
      CScoeff_i[i] = tmpl1;
      CScoeff_j[i] = tmpl2;
      CScoeff_cos[i] = tmp1;
      CScoeff_sin[i] = tmp2;
      //printf("%ld %ld %g %g\n",CScoeff_i[i],CScoeff_j[i],CScoeff_cos[i],CScoeff_sin[i]);
      //      exit(0);
    }
  fclose(fp);
  
  // load check points
  if((fp = fopen("../PIAAM2_check_fine.dat","r"))==NULL)
    {
      printf("ERROR: could not read file\n");
      exit(0);
    }

  for(i=0;i<NBcheckpts;i++)
    {
      fscanf(fp,"%lf %lf %lf %lf %lf %lf\n",&tmp1,&tmp2,&tmp3,&tmp4,&tmp5,&tmp6);
      xsample[i] = tmp1;
      ysample[i] = tmp2;
      f2sample[i] = tmp4;
      OAP2sample[i] = tmp5;
      OAT2sample[i] = tmp6;
      M2sample[i] = tmp3;
    }
  fclose(fp);

  fp = fopen("fitresidualM2.dat","w");
  for(i=0;i<NBcheckpts;i++)
    {
      x = xsample[i];
      y = ysample[i];
      r = sqrt(x*x+y*y);
      X = x/0.045;
      Y = y/0.045;
      R = sqrt(X*X+Y*Y);

      f2 = -r*r/4.0/f;
      for(k=0;k<200;k++)
	f2 += bk[k]*cos(k*(r/0.045)*PI);
      
      OAP2 = -((OAD+x)*(OAD+x)+y*y-OAD*OAD)/(2.0*(f+sqrt(f*f+OAD*OAD)))+(x*x+y*y)/(4.0*f);
      
      OAT2 = 0.0;
      for(k=0;k<NBcscoeff;k++)
	OAT2 -= 0.5*(CScoeff_cos[k]*cos(X*CScoeff_i[k]+Y*CScoeff_j[k])+CScoeff_sin[k]*sin(X*CScoeff_i[k]+Y*CScoeff_j[k]));
      
      fprintf(fp,"%.9f %.9f %g %g %g %g\n",x,y,f2sample[i]-f2,OAP2sample[i]-OAP2,OAT2sample[i]-OAT2,M2sample[i]-(f2+OAP2+OAT2));
    } 
  fclose(fp);

  return(0);
}