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

#define size  1024*1024

char head[72][80];
char f1[30];
float bias[size],flat[size];		/* at last, bias as out put */
short ccd[size];
FILE *fp;

main(ac,av)
int ac; char *av[];                        
{
  float f64,fmval;
  int i,k,n1,n2,n,m1,m2,ic,ir,ipos;	/* here, int is 4 bytes length */

  if(ac<3)stop("\tUsage:  flaten_300 datafile flatfile\0");

  strcpy(f1,av[1]);		/* data */
  if(inqh(f1,&n1,&n2,&f64,&fmval)==0)stop("DATA not found!\0");
  ipos=indexpos(head,"CRVAL1  ",72);
  sscanf(&head[ipos][23],"%d",&ic); 
  sscanf(&head[ipos+1][23],"%d",&ir); 
  m2 = (n1*n2>size) ? size/n1 :n2;
  if(n2 != m2)printf(" I change row number, because of your big file\n");
  m1=n1;
  n=m1*m2;

  strcpy(f1,av[2]);				/* for flat */
  k=inqh(f1,&n1,&n2,&f64,&fmval);
  if(k==0)stop("FLAT not found\0");
  if(f64!=0)stop("do SUPERFLAT first!\0");
  if(n1!=2048 || n2!=2048)stop("FLAT not 2k*2k\0");
  printf("  %s\t\tmean:%9.2f\n",f1,fmval);
  get300(f1,flat,m1,m2,ic,ir,k);

  strcpy(f1,av[1]);				/* for data */
  i=inqh(f1,&n1,&n2,&f64,&fmval);
  fp=fopen(f1,"rb"); fread(head,80,72,fp);
    ipos=indexpos(head,"END     ",72);
    if(ipos==72){
      fread(ccd,80,36,fp);
      for(i=0;i<80;i++)head[71][i]=32;
      strncpy(&head[71][0],"END",3);
    }
  fread(ccd,2,n,fp);
  swap2(ccd,2*n);
  fclose(fp);

  for(i=0;i<n;i++) bias[i] =(ccd[i]+f64)/flat[i];
  strncpy(&head[1][27],"-32",3);	/* real*4 */
  sprintf(&head[4][22],"%8d",m2); head[4][30]=32;
  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,"PZERO   ",72);
  if(ipos==72)ipos=indexpos(head,"BZERO   ",72);
  strncpy(&head[ipos][20],"  0.",4);
  head[ipos][0]='P'; head[ipos+1][0]='P'; 
  unlink(f1);
  printf("written data:  %s\n",f1);
  fp=fopen(f1,"wb");
  fwrite(head,80,72,fp);
  swap4(bias,4*n);
  fwrite(bias,4,n,fp);
  fclose(fp);
  exit(0);
}

get300(f1,a,n1,n2,ic,ir,l)
char f1[];
float a[];
int n1,n2,ic,ir,l;
{
  long i;
  int j,k=0;
  fp=fopen(f1,"rb"); 
  for(j=0;j<n2;j++,k+=n1){
    i=l*2880+((ir+j)*2048+ic)*4;
    fseek(fp,i,0);
    fread(&a[k],4,n1,fp);
    swap4(&a[k],4*n1);
  }
  fclose(fp);
}

inqh(f1,n1,n2,f64,fm)
char f1[];
int *n1,*n2;
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); 
  *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); 
  ipos=indexpos(head,"END     ",72);
  if(ipos==72)return(3);
  return(2);
}