/* last reading s, for new file use it's head	*/
/*				 1994,5,5	*/
/* 1994,5,27 for sun */
/*  s(i)=(s(i)+f64)/f(i)           94.9	ignore bias		*/ 
/* 94,10,4 add *.fit */
/* 94,10,5 modify BZERO --> PZERO, BSCALE --> PSCALE */
/* 94,11,17 fitshead */
/* 95,3,28 add
    ipos=indexpos(head,"END     ");
    if(ipos==72)fread(flat,80,36,fp3);
*/
#include <stdio.h>
#include <sys/types.h>
#include <sys/stat.h>

#define size  4096
#define nbyte size*size*4

char head[72][80],dummy[36*80];
char a[nbyte];
char f1[30];
FILE *fp,*fp1,*fp3;
//float flat[size][size];
short ccd[size+64];

main(ac,av)
int ac; char *av[];                        
{
  float f64,fmval;
  int   i,ipos;			/* here, int is 4 bytes length */
  int   n1,n2,c1,c2,n11,n22,c11,c22;
  if(ac<3)stop("\tUsage:  flaten datafile flatfile   (for any size,60222)\0");

  strcpy(f1,av[2]);				/* for flat */
  if(inqh(f1,&n1,&n2,&c1,&c2,&f64,&fmval)==0)stop("FLAT not found\0");
  if(f64!=0)stop("do SUPERFLAT first!\0");
  if(n1!=size || n2!=size)stop("FLAT not 4k*4k\0");
  printf("  %s\t\tmean:%9.2f\n",f1,fmval);
  fp3=fopen(f1,"rb"); fread(head,80,72,fp3);
    ipos=indexpos(head,"END     ",72);
    if(ipos==72)fread(dummy,80,36,fp3);
  strcpy(f1,av[1]);				/* for data */
  if(inqh(f1,&n11,&n22,&c11,&c22,&f64,&fmval)==0)stop("DATA not found\0");
  printf("size: %d %d %d %d     %d %d %d %d\n",n1,n2,c1,c2,n11,n22,c11,c22);
  if(c11<c1 || c22<c2 || c11+n11>c1+n1 || c22+n22>c2+n2){
    printf("Error: data frame out of FLAT frame !!!\n");
    exit(0);
  }    
  fp1=fopen(f1,"rb"); fread(head,80,72,fp1);
    ipos=indexpos(head,"END     ",72);
    if(ipos==72){
      fread(dummy,80,36,fp1);
      for(i=0;i<80;i++)head[71][i]=32;
      strncpy(&head[71][0],"END",3);
    }
  div_flat(a,f64,n1,n2,c1,c2,n11,n22,c11,c22); printf("\n");
  strncpy(&head[1][27],"-32",3);	/* real*4 */
  i=index(f1,'.');
  strcpy(&f1[i],".fit");
  f1[0]='p';
  ipos=indexpos(head,"STATUS  ",72);
  strncpy(&head[ipos][11],"SUPER STAR",10);
  ipos=indexpos(head,"BZERO   ",72);
  if(ipos==72)ipos=indexpos(head,"PZERO   ",72);
  strncpy(&head[ipos][20],"  0.0000",8);
  unlink(f1);
  printf("written data:  %s\n",f1);
  fp=fopen(f1,"wb");
  fwrite(head,80,72,fp);
  swap4(a,4*n1*n2);
  fwrite(a,4,n1*n2,fp);
  fclose(fp);
  exit(0);
}

div_flat(a,f64,n1,n2,c1,c2,n11,n22,c11,c22)
float a[],f64;
int n1,n2,c1,c2,n11,n22,c11,c22;
{
  float *p,*flat;
  int i,j,ki,kj;
  flat=a;
  fread(flat,4,n1*n2,fp3);
  swap4(flat,4*n1*n2);
  fclose(fp3);
  p = a; 
  kj=c11-c1;  
  for(i=0;i<n22;i++){
    if(i%100==0){ printf("."); fflush(stdout); }
    fread(ccd,2,n11,fp1); swap2(ccd,2*n11);
    ki=(i+c22-c2)*size+kj;
//   for(j=0;j<n11;j++) *p++ =(ccd[j]+f64)/flat[ki][j+kj];
     for(j=0;j<n11;j++) *p++ =(ccd[j]+f64)/flat[ki+j];
  }
  fclose(fp1); 
}
	
inqh(f1,n1,n2,c1,c2,f64,fm)
char f1[];
int *n1,*n2,*c1,*c2;
float *f64,*fm;
{
  int k,ipos;
  k=index(f1,'.');
  if(k==-1) strcat(f1,".fit");
  if((fp=fopen(f1,"rb"))==NULL)return(0);
  fread(head,80,72,fp); fclose(fp);
  sscanf(&head[3][21],"%d",n1); 	/* &*n1 == n1 */
  sscanf(&head[4][21],"%d",n2); 
  ipos=indexpos(head,"CRVAL1   ",72);
  sscanf(&head[ipos][21],"%d",c1); 
  sscanf(&head[ipos+1][21],"%d",c2); 
  *fm= -99.;
  ipos=indexpos(head,"VOLT7   ",72);
  if(ipos==72)ipos=indexpos(head,"MEANVAL ",72);
  if(nindex(&head[ipos][0],"MEANVAL")==0)sscanf(&head[ipos][20],"%f",fm);
  ipos=indexpos(head,"BZERO   ",72);
  if(ipos==72)ipos=indexpos(head,"PZERO   ",72);
  sscanf(&head[ipos][20],"%f",f64); 
  return(1);
}