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

/*
  convert from pixel to wavelength using fit-function 
       for 2-prism spectroscopy of high-speed spectrograph

  ver 1.0                             developped by M.Isogai on 2008 07/09
  ver 1.1  change start and end condition           modified on 2008 07/10
  ver 1.2  change end condition (add parameter)              on 2008 07/17
  ver 1.3  change calculation range (i)                      on 2008 07/18

  execute as
   %  cpwv13 shift[pix] [output-file-name]

*/

double fitfunc(double a, double b, double c, double d, double e, double x)
{
  double out;
  out = a/pow( fabs(b+x) ,d) + c + e*x;
  return out;
}

int main(int argc, char *argv[])
{
  FILE *outfile;
  int i,imax,is,ie;
  double a,b,c,d,e;
  double diff,y1,y2,y3,x;

  if (argc != 3)
    {
      printf("Help: execute as: cpwv13 shift[pix] [output-file-name]\n");
      exit(1);
    }

  diff = atof(argv[1]);
  outfile = fopen(argv[2] ,"w");

  // best-fit parameter with ZEMAX-data-file
  a = 2227.07;
  b = -279.986;
  c = 29.3246;
  d = 0.366173;
  e = -0.272483;

  // start and end of pixel position
  is=121;
  ie=290;

  for(i=is;i<ie+1;i++)
    {
      x = i+diff;
      y1 = fitfunc(a,b,c,d,e,x);
      y2 = fitfunc(a,b,c,d,e,x+0.5);
      y3 = fitfunc(a,b,c,d,e,x-0.5);

      if( y1 > 349. && y1 < 1100. && y2-y3 > 0.)
	{
	  fprintf(outfile,"%3d  %8.3f  %8.4f  %8.2f %8.3f\n",i,y1,y2-y3, y1*10., (y2-y3)*10.);
	}
    }
  fclose(outfile);
}