#include <stdio.h>
#define n1 4096
#define n2 4096
short d[1024*1024];
float a[n2*n1],b[n2*n1],c[n2*n1];
char head[72][80],head0[72][80];
FILE *fp,*fp1;
main(int ac, char **av)
{
  int i,j,k,j1,j2,aftsig;
  char f1[60],cc;
  float x,y,sky,sig,y1,y2,y3,oldsig;
  j=n1*n2;
  if(ac<2){ printf("\n\t *** q*.fit  = p*.fit- (SF_?.fit- 1.) *const *****\n");
            printf("\n\t SF_?.fit = super_?.fit/F*?.fit (use div_fit)\n");
            printf("\n\t psc1 p*.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,4,j,fp);
  if(k!=j){ printf("%s not a 4k*4k_short file!\n",av[1]); exit(0); } 
  swap4(a,4*k); fclose(fp);
  strcpy(f1,av[1]); i=strlen(f1); cc=f1[i-8];
  if(ac>2) printf("color: %c\n",cc);
  strcpy(f1,"S0_?.log");  f1[3]=cc;
  fp1=fopen(f1,"a+");
  strcpy(f1,"/vega2/rhbin/s0_?.fit"); f1[16]=cc;
  fp=fopen(f1,"rb"); if(fp==0){ printf("%s not found!\n",f1); 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",f1); exit(0); } 
  swap4(b,4*k); fclose(fp);
  y2=100;  
  for(x=-10;x<=900;x+=10){
    for(i=0;i<k;i++)c[i]=a[i]-(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);
    if(ac>2)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);
  fprintf(fp1,"%s ",av[1]);
  fprintf(fp1,"con=%6.1f sigma=%6.2f (%6.2f%6.2f)   min:const=%6.1f\n",y,y2,y1,y3,x);
  if(x>0.)for(i=0;i<k;i++)c[i]=a[i]-(b[i]-1.)*x;
  else {x=0.;  printf("\t\t\t no change at all\n");
          fprintf(fp1,"\t\t\t no change at all\n"); }
  strcpy(f1,av[1]); f1[0]='q';
  sprintf(&head[0][33],"adjust const.=%7.3f",x); head[0][54]=32;
  fp=fopen(f1,"wb"); fwrite(head,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;
}