#include <stdio.h>
#define n1 4096
#define n2 4096
short a[n2*n1],d[1024*1024];
float b[n2*n1],c[n2*n1];
char head[72][80],head0[72][80];
FILE *fp;
main(int ac, char **av)
{
  int i,j,k,j1,j2,aftsig;
  char f1[60];
  float x,y,bzero,sky,sig,y1,y2,y3,oldsig;
  j=n1*n2;
  if(ac<3){ printf("\n\t *** q*.fit  = d*.fit-  (super_?.fit- 1.) *const *****\n");
            printf("\n\t dsc d*.fit super_?.fit\n\n");
            exit(0); }
  fp=fopen(av[1],"rb"); if(fp==0){ printf("%s not found!\n",av[1]); exit(0);}
  fread(head,72,80,fp); k=fread(a,2,j,fp);
  if(k!=j){ printf("%s not a 4k*4k_short file!\n",av[1]); exit(0); } 
  swap2(a,2*k); fclose(fp);
  k=indexpos(head,"BZERO   ",72);
  sscanf(&head[k][20],"%f",&bzero); 
  fp=fopen(av[2],"rb"); if(fp==0){ printf("%s not found!\n",av[2]); exit(0);}
  fread(head0,72,80,fp); k=fread(b,4,j,fp);
  if(k!=j){ printf("%s not a 4k*4k_float file!\n",av[2]); exit(0); } 
  swap4(b,4*k); fclose(fp);
  y2=100;  
  for(x=-10;x<=500;x+=10){
    for(i=0;i<k;i++)c[i]=(a[i]+bzero)-(b[i]-1.)*x;
    j=0; for(j1=2048-500;j1<2048+500;j1++)  
         for(j2=2048-500;j2<2048+500;j2++)  
         d[j++]=c[j1*4096+j2];
    white_black(d,1,j,&sky,&sig);
//    printf("x=%6.1f  sky=%6.1f  sigma=%6.2f\n",x,sky,sig);
    if(aftsig){ y3=sig; aftsig=0; }
    if(sig<y2){ y2=sig; y=x; y1=oldsig; aftsig=1; }
    oldsig=sig;
    if(sig>1.5*y2)break;
  }
  paowu(y1,y2,y3,y,10.,&x);
  printf("%s ",av[1]);
  printf("con=%6.1f sigma=%6.2f (%6.2f%6.2f)   min:const=%6.1f\n",y,y2,y1,y3,x);
  for(i=0;i<k;i++)c[i]=(a[i]+bzero)-(b[i]-1.)*x;
  strcpy(f1,av[1]); f1[0]='q';
  sprintf(&head0[0][33],"adjust const.=%7.3f",x); head0[0][54]=32;
  fp=fopen(f1,"wb"); fwrite(head0,72,80,fp); swap4(c,4*k); 
  fwrite(c,4,k,fp); fclose(fp);
}

paowu(y1,y2,y3,x2,step,x)
float y1,y2,y3,x2,step,*x;
{
  float y;
  y=(y1-y3)/((y1+y3)*.5-y2)*0.25;
  *x=x2+y*step;
}