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

#define PI 3.14159265358979323846264338328
#define NBcheckpts 31341

double ak[200];
double f = 0.9;
double OAD = 0.085;


double xsample[31341];
double ysample[NBcheckpts];
double f1sample[NBcheckpts];
double OAP1sample[NBcheckpts];
double OAT1sample[NBcheckpts];
double M1sample[NBcheckpts];

double CX = 6.4339e-07;
double CXX = 1.15783e-06;
double CYY = 1.10203e-06;
double C1 = 5.20341e-07;
double C2 = 3.42046;
double C3 = 0.852328;

int main()
{
  FILE *fp;
  long i,k;
  double f1,OAP1,OAT1;

  double f1_err,OAP1_err,OAT1_err;

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

  // load ak
  fp = fopen("../PIAAM1_cosfitcoeff_200.dat","r");
  for(i=0;i<200;i++)
    {
      fscanf(fp,"%ld %lf\n",&tmpl,&tmp1);
      ak[tmpl] = tmp1;
    }
  fclose(fp);
  
  // load check points
  if((fp = fopen("../PIAAM1_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;
      f1sample[i] = tmp4;
      OAP1sample[i] = tmp5;
      OAT1sample[i] = tmp6;
      M1sample[i] = tmp3;
    }
  fclose(fp);

  fp = fopen("fitresidualM1.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);

      f1 = r*r/4.0/f;
      for(k=0;k<200;k++)
	f1 += ak[k]*cos(k*(r/0.045)*PI);
      
      OAP1 = ((OAD-x)*(OAD-x)+y*y-OAD*OAD)/(2.0*(f+sqrt(f*f+OAD*OAD)))-(x*x+y*y)/(4.0*f);

      OAT1 = CX*X+CXX*X*X+CYY*Y*Y+C1*sin(C2*R)*(1.0+C3*R*R)*x/r;
      
      fprintf(fp,"%.9f %.9f %g %g %g %g\n",x,y,f1sample[i]-f1,OAP1sample[i]-OAP1,OAT1sample[i]-OAT1,M1sample[i]-(f1+OAP1+OAT1));
    } 
  fclose(fp);

  return(0);
}