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


int order(int ij, int ik);

void main(argc,argv)
     char **argv;
     int argc;

{ int i=0,j=0,k=0,ordpos,ord1=3,ord2=6,maxbin;
  float binflux[12],maxflux,fluxsp[10][256],fluxsky[10][256];
  short int frame1[65536],frame2[65536];
  FILE *infile1,*infile2,*outfile,*diag;

  printf("%d \n",argc);
  if (argc<2) 
      {printf("\7\n Usage: real3 cm000 cm001\n");
      printf("          ^data  ^background\n");
      printf("    Background can be omitted.\n\n");

               exit(1);
	     }

  infile1 = fopen(argv[1],"r");
  if (argc==3){ infile2 = fopen(argv[2],"r");
                fseek(infile2, 2880, 0);  }

  strcat(argv[1],".sp");
  outfile = fopen(argv[1],"w");

  fseek(infile1, 2880, 0);

  diag = fopen("diag","w");

  /*get read1 from data frame*/

  fread(frame1, 2, 65536, infile1);
  swab(frame1, frame1, 2*65536);

  /*get read2 from data frame*/

  fread(frame2, 2, 65536, infile1);
  swab(frame2, frame2, 2*65536);

  /* create read2 - read1 */

  for (i=0;i<65536;i++){frame2[i]=frame2[i]-frame1[i];}

  if (argc==3){

  /* get read1 from background frame */

  fread(frame1, 2, 65536, infile2);
  swab(frame1, frame1, 2*65536);

  /* subtract background (-read1) from data frame (i.e. add it) */

  for (i=0;i<65536;i++){frame2[i]=frame2[i]+frame1[i];}

  /* get read2 from background frame */

  fread(frame1, 2, 65536, infile2);
  swab(frame1, frame1, 2*65536);

  /* subtract background read2 from data frame */

  for (i=0;i<65536;i++){frame2[i]=frame2[i]-frame1[i];}
  }

  /* frame difference is now ready for analysis */
 
  /* find 4th order spectrum position along slit */

  for (i=0;i<12;i++) {binflux[i]=0.;}

  for (i=100;i<200;i++){
    for (j=4;j<=4;j++){
       ordpos = order(j,i);
       printf("  %4d %4d\n",i,ordpos);
       for (k=0;k<12;k++) {printf("%5d ",frame2[256*(ordpos-6+k)+i]);
                           binflux[k] += frame2[256*(ordpos-6+k)+i];}
       }
      printf("\n");
    }
    for (i=0;i<12;i++){printf("%3d %9.0f \n",i,binflux[i]);}
    maxflux=binflux[0];
    maxbin=0;
    for (i=0;i<12;i++){if (binflux[i]>maxflux) {maxbin=i;maxflux=binflux[i];}}
    printf("maxbin=%3d\n",maxbin);   

 /*extract 3rd through  6th order spectra*/

    for (j=0;j<10;j++){
      for (i=0;i<256;i++){fluxsp[j][i]=0.;}}
    for (i=0;i<256;i++){
      for (j=ord1;j<=ord2;j++){
         ordpos = order(j,i);
         printf("  %4d %4d  ",i,ordpos);
         for (k=-1;k<4;k++) 
               {fluxsp[j][i] += frame2[256*(ordpos-7+maxbin+k)+i];
                printf("%6d ",frame2[256*(ordpos-7+maxbin+k)+i]);
		frame2[256*(ordpos-7+maxbin+k)+i]+=j*1000;
               }
         }
     printf("\n"); 
    }

  /*extract sky*/

  for (j=0;j<10;j++){
      for (i=0;i<256;i++){fluxsky[j][i]=0.;}}
  for (i=0;i<256;i++){
    frame2[i]=i;
    for (j=ord1;j<=ord2;j++){
       ordpos = order(j,i);
       for (k=0;k<10;k++) {
          if (k<maxbin-5 || k>maxbin-1){ 
               fluxsky[j][i] += frame2[256*(ordpos-3+k)+i];
               printf("%6d ",frame2[256*(ordpos-3+k)+i]);
               frame2[256*(ordpos-3+k)+i]+= 9999;}
          }
       printf("\n");
       }
     }

    for(i=0;i<65536;i++){frame1[i]=0;}
    fwrite(frame1,2,65536,diag);
    fwrite(frame2,2,65536,diag);
    fclose(diag);

    /* print results */

    for(j=ord1;j<=ord2;j++){
        for (i=0;i<256;i++) {
                  	  fprintf(outfile,"%4d %7.0f %7.0f %7.0f \n",j*256+i,fluxsp[j][i],fluxsky[j][i],fluxsp[j][i]-fluxsky[j][i]);}
          }


} 

int order(int i, int j)
  {
    float rj,x,y;
    rj=j-1;
    if (i==3) return(61.11 + 0.409639*rj - 0.000139167*rj*rj - .5);
    if (i==4) return(35.97 + 0.438707*rj - 0.000258542*rj*rj - .5);
    if (i==5) return(10.51 + 0.496111*rj - 0.000373005*rj*rj - .5);
    if (i==6) { x= -18.05 + 0.573144*rj - 0.00048111*rj*rj - .5;
                if (x<5.) x=5.;
                return(x);}
  }