#include <stdio.h>
#define n1  4096
#define n2  4032
#define f48 348         // max join_do_flat file is 348
char  head[72][80], f1[30];
float a[n2][n1],ccd[n1*f48],fm[f48], gg;
FILE  *fp,*f[f48];		
int   imode,ip=0;

main(ac,av)
int ac; char *av[];                        
{
  int i,j,k;			 
  if(ac<3){
    printf("\n\t**** get a median(mean) bias after e4(-overscan)\n");
    printf("\n\tUsage: b4 batch outfit [mode]\n\n");
    exit(0);
  }
  strcpy(f1,av[1]);		/* batch file */
  fp=fopen(f1,"r");
  if(fp==0){ printf("\n\tbatch_file not found: %s\n",f1); exit(0); }
  for(i=0;i<f48;i++){
    fscanf(fp,"%s\n",f1); printf("\n%4d:   %s  ",++ip,f1);
    if(f1[0]!='e'){ printf("\n\tfile must be e*_?.fit\n\n"); exit(0); }
    f[i]=fopen(f1,"rb"); if(f[i]==0){ printf("not found!\n"); exit(0); }
    fread(head,80,72,f[i]);
    sscanf(&head[3][24],"%d",&j);
    sscanf(&head[4][24],"%d",&k);
    if(j!=n1 || k!=n2){ printf("not a 4096*4032 file\n"); exit(0); }
    k=indexpos(head,"SKYADU  ",72);
    sscanf(&head[k][20],"%f",&fm[i]);  printf("%9.2f",fm[i]);
    if(feof(fp))break;
  } printf("\n");
  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      ");
 
  work(a);
  strncpy(&head[1][27],"-32",3);	/* real*4 */
  strcpy(f1,av[2]); i=nindex(f1,".");
  if(i>0)strcpy(&f1[i],".fit"); else strcat(f1,".fit");
  k=indexpos(head,"SKYADU  ",72);
  gg=0.; for(i=0;i<ip;i++)gg+=fm[i];  
  sprintf(&head[k][9]," %20.2f ",gg); head[k++][31]='/';
  sprintf(&head[k][0],"SUP_BIAS=%3d biasfiles",ip); 
  if(imode == 1)strncpy(&head[k][21],"_median",7);
  if(imode == 0)strncpy(&head[k][21],"_mean  ",7);
  unlink(f1);
  printf("produce: %s\n",f1);
  fp=fopen(f1,"wb");
  fwrite(head,80,72,fp);
  swap4(a,n1*n2*4);
  fwrite(a,4,n1*n2,fp);
  fclose(fp);
}

work(a)
float *a;
{
  float bb[f48],x;
  int i,j,k,l;
  for(i=0;i<58;i++)printf("-"); printf("\n");
  for(i=0;i<n2;i++){
    if(i%99==0){ printf("."); fflush(stdout); }
    for(j=0;j<ip;j++){ k=j*n1; fread(&ccd[k],4,n1,f[j]); }
    swap4(ccd,n1*ip*4);		
    for(j=0;j<n1;j++){
      l=j; for(k=0;k<ip;k++){  bb[k]=ccd[l]; l+=n1; }
      if(imode == 0)sup_mean(bb,ip,&x);
      if(imode == 1)sup_median(bb,ip,&x);
      *a++ =x;
    }
  }
  for(j=0;j<ip;j++)fclose(f[j]); 
}
	
sup_mean(a,nn,mean)
float *a,*mean; int nn;
{
  int i; float x;
  x=0.; for(i=0;i<nn;i++) x+=a[i];
  *mean= x/nn;
}
	
sup_median(a,nn,medi)
float *a,*medi; int nn;
{
  float x; int i,j;
  for(i=0;i<nn-1;i++)for(j=i+1;j<nn;j++)if( *(a+i) > *(a+j)){
    x= *(a+i); *(a+i) = *(a+j); *(a+j)=x;
  }
  i=nn/2;
  *medi = (i+i==nn) ? (*(a+i) + *(a+i-1))*0.5 : *(a+i);
}