#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "comphys.h"
#include "comphys.c"
void sort(unsigned long n, float arr[]);
float *scany,*scanx,*factorx,*factory;
float *t;
float fx(float x[]);
float fy(float x[]);

/* Note: Codes 0 = blank image, 1 = star found, but not on slit,
   2 = star found and it is on (but may be decentered) the slit.
*/

int main(int argc, char *argv[])
{
  float Backx[382] = 
    {235461,235160,234701,235266,234043,234154,234886,233602,233176,234855,
     233675,232935,232735,232930,232889,232849,232969,232821,232438,232695,
     232539,232605,232184,232331,232406,232690,231730,232186,231565,231736,
     231971,232244,231330,232770,231183,232494,231715,231783,231614,231281,
     231530,231747,230893,231272,230294,231107,231758,230633,231038,231749,
     231188,230847,231380,230410,230759,231026,231740,230935,230674,230543,
     230711,230639,230444,230328,230461,229911,230484,230191,229703,229998,
     230086,229912,229958,230210,229946,229823,230065,229563,229671,229605,
     230011,230004,229175,229406,229836,229482,229575,229126,229531,229237,
     229361,228963,229784,228960,230381,228477,229911,229054,229318,229012,
     228763,228737,228667,228822,229061,229122,228375,228532,229040,228659,
     228396,228529,228558,228421,228274,228366,229155,229424,228704,227996,
     228294,228304,228276,228330,228590,228392,228415,228208,227462,227758,
     228319,228361,227939,228026,228216,228268,227842,227671,228374,227717,
     228061,227682,227703,227893,227883,227623,227735,227351,227628,228085,
     227801,227788,227652,226932,227980,227536,228448,227670,227259,228028,
     226832,227946,227345,227146,227625,227147,227726,227137,227255,227463,
     227326,227663,230613,227171,226864,226904,227435,227227,227239,227455,
     227688,227883,227073,226866,227292,227251,227056,226909,227081,226585,
     226939,227142,227110,228345,227361,226810,228572,226526,227229,227032,
     226862,227351,226874,227089,227013,227095,226922,226851,227329,226890,
     226773,227940,226981,226944,227088,226351,227358,226625,227058,226798,
     226714,227061,226563,226261,226393,226408,226705,226962,227482,226253,
     226643,226411,226293,226724,226576,226724,226459,226683,226429,226740,
     226478,226405,226880,226510,226712,226590,226427,226704,226159,226864,
     226497,226841,226586,226108,226277,226276,226870,226198,226784,226292,
     226912,226296,226358,226525,226183,225863,226502,226469,226788,226359,
     225910,226362,227672,226697,226371,226690,226255,226059,226511,226086,
     227535,225927,226681,226271,226427,226164,225964,226272,226340,225807,
     226102,225994,226240,227177,226716,225996,226668,226199,225779,225877,
     226066,226696,226133,225673,226346,226177,225893,227194,225574,226304,
     225587,226596,226711,226299,226122,227826,226007,226355,226372,226030,
     225952,227796,226684,226356,226743,226058,225957,225745,226143,225615,
     226286,226418,225718,226331,225742,225881,228213,225928,226093,226096,
     226763,226352,226206,230222,226511,226029,225794,226657,225864,225903,
     226274,226055,226178,225706,225670,227656,226283,226242,226515,225627,
     225889,226686,225528,226012,226501,225885,225918,225774,226296,226387,
     225786,225830,225824,225276,225596,225785,225720,225899,225884,226357,
     226484,225657};
  float Backy[255] = 
    {340257,340915,340292,340312,340391,340295,340513,340388,340289,340755,
     340253,340653,341498,340715,341152,340827,340804,340736,340281,340558,
     340634,340831,340897,340784,341013,340734,341030,340832,340953,340732,
     341322,341457,341203,340420,341117,340598,340644,341503,341457,340923,
     341064,340554,341403,340943,341108,341254,340915,341330,340461,341285,
     341133,341860,341030,340900,341019,340844,341537,341303,342181,341179,
     340570,341474,341485,341020,341445,341466,340987,342347,341196,340830,
     341941,341499,340840,340603,342877,341465,341603,341410,340919,340893,
     341478,340863,341005,341467,341349,340685,342714,341966,341478,341165,
     341280,341787,341827,341080,341395,341285,341019,341196,341656,341233,
     341282,341311,341154,341614,343214,340930,341805,341566,341288,342110,
     342405,341746,342335,341698,341459,341452,341431,341626,341010,341502,
     341478,341508,341571,341837,341625,342412,342559,341981,341395,342384,
     341634,341408,341220,342258,341547,341716,341657,342493,341811,341243,
     341333,341996,342335,341243,342608,342299,342382,341790,342269,341837,
     341327,342084,341644,342013,341669,341337,341574,341620,342180,341822,
     341122,342190,344282,342619,341649,341484,342208,341887,342376,341812,
     341640,342438,341769,341826,341723,341598,342115,342494,341801,342321,
     342579,341946,343508,342010,341676,341814,341514,341932,342451,342200,
     342163,342217,342210,342150,342144,341892,341597,342345,342184,342135,
     341750,341547,341828,342168,342249,342212,341902,341791,341822,341889,
     341911,341668,342245,341958,342060,342149,342249,342071,342173,342594,
     343053,342105,342678,341832,342105,342810,342893,342020,343485,342066,
     342316,342891,341927,341986,342331,342019,342458,341738,341696,341662,
     342702,342197,342113,342073,341970,342397,342290,342193,342317,342042,
     342074,344292,342411,342325,342325};
  float **array;
  float *stat,*tempy,*backx,*backy;
  float median,max,maxm1,sumx,sumy,sumM,cx,cy,smax,smax2,xmax2,xmin,xav;
  int xmax;
  float medianx,mediany;
  float smay,ymax,ymax2;
  float *p,**xi,fret;
  float xcoord,ycoord;
  int iter,k;
  float A,B,C,D,ypos;
  float threshold = 1500.0;
  int i,j,l,I,J,fwhm,err,a,b,c,m,n;
  int bx,by;
  int minfound = 0;
  char buffer[5000],buf[5000],*tmp,infile[80];
  FILE *in,*out,*outy,*bad;

  backx = Backx - 1;
  backy = Backy - 1;
  array = matrix(1,382,1,255);
  scany = vector(1,255);
  tempy = vector(1,255);
  scanx = vector(1,382);
  factorx = vector(1,382);
  factory = vector(1,255);
  t = vector(1,382);
  stat = vector(1,97410);
  p = vector(1,4);
  xi = matrix(1,4,1,4);

  out = fopen("scanx.out","w");
  outy = fopen("scany.out","w");

  strcpy(infile,argv[1]);

  sprintf(buf,"imlist %s > temp.dat",infile);
  err = system(buf);

  for(i=1;i<=382;i++) t[i] = (float)i;

  in = fopen("temp.dat","r");
  if(in == NULL) {
    printf("Cannot open temp.dat\n");
    exit(1);
  }

  fgets(buf,4000,in);
  fgets(buf,4000,in);

  for(j=255;j>=1;j--) {
    fgets(buffer,4000,in);
    tmp = strtok(buffer," ");
    for(i=1;i<=382;i++) {
      tmp = strtok(NULL," ");
      array[i][j] = atof(tmp);
    }
  }
  fclose(in);

  /* Correct for badpixels by reading badpix.dat file.  Skip if it can't
     find the file */

  bad = fopen("badpix.dat","r");
  if(bad != NULL) {
    while(fscanf(bad,"%d,%d",&bx,&by) != EOF) {
      array[bx][by] = (array[bx-1][by-1] + array[bx][by-1] + array[bx+1][by-1]
                    + array[bx-1][by] + array[bx+1][by] + array[bx-1][by+1] + 
		       array[bx][by+1] + array[bx+1][by+1])/8.0;
    }
    fclose(bad);
  }


  /* Get y-direction scan by adding up rows in x-direction.  This scans
    over the entire chip.  We may want to change that if we can routinely
    get the star within, say, 100 pixels of the slit, so that if there
    are other stars in the field there is less chance of confusion */
  for(j=1;j<=255;j++) tempy[j] = 0.0;

  for(j=1;j<=255;j++) {
    tempy[j] = 0.0;
    for(i=1;i<=382;i++) {
	tempy[j] += array[i][j];
    }
  }
  
  //background subtraction -- y scan.

  for(i=1;i<=255;i++) factory[i] = tempy[i]/backy[i];
  sort(255,factory);
  mediany = factory[128];

  /* The median -- mediany -- is a scaling function that scales the
     stored background values to the background in the image.  This is because
     the background is a function of the temperature */

  for(i=1;i<=255;i++) tempy[i] = tempy[i] - mediany*backy[i];


  // Appy a boxcar smooth to profile, to smooth out dust spots
  scany[1] = tempy[1];
  scany[2] = tempy[2];
  for(j=3;j<254;j++) scany[j] = 
      (tempy[j-2]+tempy[j-1]+tempy[j]+tempy[j+1]+tempy[j+2])/5.0;
  scany[254] = tempy[254];
  scany[255] = tempy[255];

  for(i=1;i<=255;i++) fprintf(outy,"%d %f\n",i,scany[i]);
  fclose(outy);

  /* Get x-direction scan by adding up columns.  Again, this is over the
     entire chip */
  for(i=1;i<=382;i++) scanx[i] = 0.0;

  for(i=1;i<=382;i++) {
    scanx[i] = 0.0;
    for(j=1;j<=255;j++) {
      scanx[i] += array[i][j];
    }
  }

  // Background subtraction, x scan

  for(i=1;i<=382;i++) factorx[i] = scanx[i]/backx[i];
  sort(382,factorx);
  medianx = factorx[191];
  for(i=1;i<=382;i++) scanx[i] = scanx[i] - medianx*backx[i];

  for(i=1;i<=382;i++) fprintf(out,"%d %f\n",i,scanx[i]);
     fclose(out);

  // Find the x-coordinate of the star

  smax = 0.0;
  for(j=1;j<=382;j++) {
    if(scanx[j] > smax) {
      smax = scanx[j];
      xmax = j;
    }
  }

  // No star detected if maximum < threshold
  if(smax < threshold) {
    printf("0.00,0.00,0,0.00\n");
    exit(0);
  }

  // If star is in border region, return no star
  if(xmax < 15 || xmax > 367) {
    printf("0.00,0.00,0,0.00\n");
    exit(0);
  }

  /* Check -- did we just discover a cosmic ray?  Check values +/- 2 pixels
     from this max.  If their average is > 0.5*smax, we have a star.  If not,
     we have a cosmic ray.  Fix the cosmic ray and then try again. */

  xav = (scanx[xmax-2] + scanx[xmax+2])/2.0;
  while(xav < 0.5*smax) {
    scanx[xmax] = 0;
    smax = 0.0;
    for(j=1;j<=382;j++) {
      if(scanx[j] > smax) {
	smax = scanx[j];
        xmax = j;
      }
    }
    xav = (scanx[xmax-2] + scanx[xmax+2])/2.0;
  }

  //  printf("Xmax = %f  smax = %e\n",xmax,smax);

  /* The Initial Point */
  p[1] = A = smax;
  p[2] = B = xmax;
  p[3] = C = 10.0;
  p[4] = D = 0.0;

  for(i=1;i<=4;i++) {
    for(j=1;j<=4;j++) {
      if(i==j) xi[i][j] = 1.0;
      else xi[i][j] = 0.0;
    }
  }

  /* Invoking Powell */

  powell(p,xi,4,0.001,&iter,&fret,fx);

  xcoord = p[2];

  /* Note: if the powell fit for the xcoordinate of the star gives a
     negative value, or a large positive value > 382, then the gaussian
     fit is probably fitting a tail of the gaussian to a small gradient
     across the image.  This probably means that the image is blank.  Interpret
     the field as blank if the xcoord is within a margin of 15 pixels or
     outside of the chip.  flag = 0 means blank image. */ 

  if(xcoord < 15.0 || xcoord > 367.0) {
    printf("0.00,0.00,0,0.00\n");
    exit(0);
  }

  // printf("xcoord = %f\n",p[2]);

  /* Now the y-direction.  This involves determining whether or not the
     star is on the slit */


  // find maximum point

  ymax = -100.0;
  for(j=1;j<=255;j++) {
    if(scany[j] > ymax) {
      ymax = scany[j];
      l = j;
    }
  }

  // Now, search for a minimum
  // Search first in positive direction
  a = l;
  b = l+1;
  c = l+2;
  minfound = 0;
  /* Note that if the star is more than 12 pixels off the slit, it
     won't find the minimum.  This may need adjustment. */
  for(k=0;k<=12;k++) {
    if(scany[b+k] <= scany[a+k] && scany[b+k] <= scany[c+k]) {
      minfound = 1;
      xmin = m = b+k;
      break;
    }
  }
  // If not found, search in negative direction
  if(minfound == 0) {
    for(k=0;k<=12;k++) {
      if(scany[b-k] <= scany[a-k] && scany[b-k] <= scany[c-k]) {
	minfound = 1;
	xmin = m = b-k;
	break;
      }
    }
  }
  // printf("minfound = %d  %f\n",minfound,xmin);
  if(minfound == 0) {
    // Star not on slit
    // Find y-coordinate using powell's method
    p[1] = A = ymax;
    p[2] = B = l;
    p[3] = C = 10.0;
    p[4] = D = 0.0;
    powell(p,xi,4,0.001,&iter,&fret,fy);
    /* Print out position of slit, flag = 1 to indicate not on slit, and 
       another 0 because we can't determine the position of slit from this
       image.  If the star is near the edge, the powell fit can be faulty,
       so compare the location powell found with B, and adjust back if 
       necessary  */
    if(fabs(p[2]-B) > 5.0) p[2] = B;
    printf("%f,%f,1,0.00\n",xcoord,p[2]);
  } else {  // Look for a maximum -- the other half of the star
    if(m <= l) { // Look in the negative direction
      a = m-2;
      b = m-1;
      c = m;
      for(k=0;k<=10;k++) {
	if(scany[b-k] >= scany[a-k] && scany[b-k] >= scany[c-k]) {
	  n = b-k;
	  ymax2 = scany[n];
	  break;
	}
      } 
    } else { // Look in the positive direction
       a = m;
       b = m+1;
       c = m+2;
       for(k=0;k<=10;k++) {
	 if(scany[b+k] >= scany[a+k] && scany[b+k] >= scany[c+k]) {
	   n = b+k;
	   ymax2 = scany[n];
	   break;
	 }
       }
      }
    // Use a weighted mean to determine y position of star when on slit
    ypos = (n*ymax2 + l*ymax)/(ymax2+ymax);
    /* print out star position, flag = 2 to indicate on slit, and position 
       of slit.  However, if xmin is less than 50 pixels from the edge, then
       you did not find the slit */
    if(xmin < 50 || xmin > 205) {
      printf("%f,%f,1,0.00\n",xcoord,ypos);
    } else printf("%f,%f,2,%f\n",xcoord,ypos,xmin);
  }

  free_matrix(array,1,382,1,255);
  free_vector(scany,1,255);
  free_vector(tempy,1,255);
  free_vector(scanx,1,382);
  free_vector(t,1,382);
  free_vector(stat,1,97410);
  free_vector(p,1,4);
  free_matrix(xi,1,4,1,4);
  return(0);
}

float fx(float x[])
{
  float y1,chi2;
  int i;

  chi2 = 0.0;
  for(i=1;i<=382;i++) {
    y1 = x[1]*exp(-(t[i]-x[2])*(t[i]-x[2])/(2.0*x[3]*x[3])) + x[4];
    chi2 += (y1-scanx[i])*(y1-scanx[i]);
  }
  // printf("\nx2 = %f chi2 = %e",x[2],chi2);
  return(chi2);
}

float fy(float x[])
{
  float y1,chi2;
  int i;

  chi2 = 0.0;
  for(i=1;i<=255;i++) {
    y1 = x[1]*exp(-(t[i]-x[2])*(t[i]-x[2])/(2.0*x[3]*x[3])) + x[4];
    chi2 += (y1-scany[i])*(y1-scany[i]);
  }
  // printf("\nchi2 = %e",chi2);
  return(chi2);
}