/* ************************************************************
c no weight, assume exp time is equal, otherwise as it a weight
c fortran superdark.for is differ this, and this is the same of superbias.
c format of "super.par"
c ----------------------------- at top_line and the first char is "g"
c g = 3.2
c r = 6.0
c nsigma = 3.0
c ## sigma= sqrt((r/g)**2+<I>/g)
c ------------------------------ "super.par" text contain 4 lines.
c 1994,10,4 add *.fit filename
c 94,10,5 modify BZERO --> PZERO, BSCALE --> PSCALE 
c 94,11,17 fitshead
c 95,12,3   g=4.1, r =12., ns=4
*************************************************************** */
#include <stdio.h>

#define size  2048
#define nbyte size*size*4

char head[72][80];
char a[nbyte];
char f1[30];
short ccd[size*18];
float f64[18];
float gg,rr,ss;
FILE *fp,*f[18];			/*  f[18] for flat */
int imode,ip=0;

main(ac,av)
int ac; char *av[];                        
{
  float fmval;
  int i,n1,n2,ipos;			/* here, int is 4 bytes length */
  if(ac<3)stop("\tUsage:  superdark batch_file outputfile [mode]\0");

  strcpy(f1,av[1]);		/* batch file */
  fp=fopen(f1,"r");
  for(i=0;i<18;i++){
    ip++;
    fscanf(fp,"%s\n",f1); printf("\t%d:\t%s\n",ip,f1);
    if(inqh(f1,&n1,&n2,&f64[i])==0)stop(" not found!\0");
    if(n1!=size || n2!=size)stop("SIZE not 2k*2k\0");
    f[i]=fopen(f1,"rb"); fread(head,80,72,f[i]);
    if(feof(fp) != 0)break;
  }
  fclose(fp);
  if(ip < 4)stop("too less files\0");
  imode=0;
  if(ac>3){
    strcpy(f1,av[3]);
    if(nindex(f1,"med")==0 || nindex(f1,"MED")==0)imode=1;
  } else if(ip > 9)imode=1;
  if(imode == 1)printf("Median    ");
  if(imode == 0)printf("Mean      ");
  getpar(&gg,&rr,&ss);
  bmean18(a,ccd);  printf("\n");
  strncpy(&head[1][27],"-32",3);	/* real*4 */
  strcpy(f1,av[2]);
  i=index(f1,'.');
  if(i>0)strcpy(&f1[i],".fit");
  else strcat(f1,".fit");
  xmean(a,&fmval);
  printf("%s    Mean value: %8.2f\n",f1,fmval);
  ipos=indexpos(head,"VOLT7   ");
  if(ipos==72)ipos=indexpos(head,"MEANVAL ");
  strncpy(&head[ipos][0],"MEANVAL =",9);
  strncpy(&head[ipos][33],"MEAN VALUE OF DARK  FIELD",25);
  sprintf(&head[ipos][22],"%8.2f",fmval); head[ipos][30]=32;
  ipos=indexpos(head,"STATUS  ");
  strncpy(&head[ipos][11],"SUPER DARK",10);
  if(imode == 1)strncpy(&head[ipos][21],"_median",7);
  if(imode == 0)strncpy(&head[ipos][21],"_mean  ",7);
  ipos=indexpos(head,"PZERO   ");
  if(ipos==72)ipos=indexpos(head,"BZERO   ");
  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);
  fwrite(a,4,n1*n2,fp);
  fclose(fp);
  exit(0);
}

bmean18(a,ccd)
float *a;
short ccd[];
{
  float bb[18];
  float rgs,x;
  int i,j,k,l;
  printf("g= %6.2f   r= %6.2f  nsigma= %6.2f\n",gg,rr,ss);
  for(i=0;i<69;i++)printf("-"); printf("\n");
  rgs=ss*rr/gg;
  for(i=0;i<size;i++){
    if(i%30==0){ printf("."); fflush(stdout); }
    for(j=0;j<ip;j++){
      k=j*size;
      fread(&ccd[k],2,size,f[j]);
    }
    for(j=0;j<size;j++){
      l=j;
      for(k=0;k<ip;k++,l+=size)bb[k]=ccd[l]+f64[k];
      if(imode == 0)sup_mean(bb,ip,rgs,&x);
      if(imode == 1)sup_median(bb,ip,&x);
      *a++ =x;
    }
  }
  for(j=0;i<ip;j++)fclose(f[j]);
}
	
sup_mean(a,n,rgs,mean)
float a[],rgs,*mean; 
int n;
{
  int i,j,m;
  float x;
  for(i=0;i<n-1;i++)for(j=i+1;j<n;j++)if(a[i]>a[j]){
    x=a[i]; a[i]=a[j]; a[j]=x;
  }
  j=n;
l10:
  m=j; i=j/2;
  x= (i+i==j) ? (a[i]+a[i-1])*.5 : a[i];
  j=0;
  for(i=0;i<m;i++){
    a[j++]=a[i];
    if((a[i]-x) > rgs)j--;
  }
  if(m!=j && j>2)goto l10;
/* get mean value */
  x=0.;
  for(i=0;i<m;i++)x+=a[i]; 
  *mean= x/m;
}

xmean(a,x)
/* take 2:2,2:3,3:2,3:3  1/4 of whole ccd frame do statistic */
float a[2048][2048], *x;
{ 
  int i,j;
  float y=0;
  for(i=512;i<1536;i++)for(j=512;j<1536;j++)y+=a[i][j];
  *x = y/(1024.*1024.);
}
	
inqh(f1,n1,n2,f64)
char f1[];
int *n1,*n2;
float *f64;
{
  int k,ipos;
  FILE *fp;
  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,"PZERO   ");
  if(ipos==72)ipos=indexpos(head,"BZERO   ");
  sscanf(&head[ipos][20],"%f",f64);
  return(1);
}

stop(s)
char s[];
{ printf("%s\n",s); exit(0); }
  
index(a,c)
char a[],c;
{
  int i,j;
  j=strlen(a);
  for(i=0;i<j;i++)if(a[i]==c)break;
  return (i<j)?i:-1; 
}

nindex(a,b)
char a[],b[];
{
  int i,j,k;
  char c;
  j=strlen(a);
  k=strlen(b);
  c=b[0];
  for(i=0;i<j;i++)if(a[i]==c && strncmp(&a[i],b,k)==0)break;
  return (i<j)?i:-1; 
}

sup_median(a,n,medi)
float *a,*medi; int n;
{
  float x; int i,j;
  for(i=0;i<n-1;i++)for(j=i+1;j<n;j++)if( *(a+i) > *(a+j)){
    x= *(a+i); *(a+i) = *(a+j); *(a+j)=x;
  }
  i=n/2;
  *medi = (i+i==n) ? (*(a+i) + *(a+i-1))*0.5 : *(a+i);
}

getpar(g,r,nsigma)
float *g,*r,*nsigma;
{
  FILE *fp;
  char c; int i;
  if((fp=fopen("super.par","r"))==NULL){
    printf("\t I cannot find file 'super.par'\n");
    *g=4.1; *r=12.0; *nsigma=4.0;
    printf("\t assume: g=4.1; r=12.0; nsigma=4.0\n");
  } else {
    fscanf(fp,"%c%c%c%f\n",&c,&c,&c,g);
    fscanf(fp,"%c%c%c%f\n",&c,&c,&c,r);
    for(i=0;i<9;i++)fscanf(fp,"%c",&c);
    fscanf(fp,"%f",nsigma);
    fclose(fp);
  }
}

indexpos(head,f1)
char head[72][80],f1[8];
{
  int i,j;
  for(i=0;i<72;i++){
    for(j=0;j<8;j++) if(head[i][j]!=f1[j])goto l10;
    return i;
l10:
    continue;
  } 	
  return i;
}