#include <stdio.h>
//  1  3
//  2  4
FILE *fp;
char  h[144][80],head[72][80],f1[30],g[9],f3[30];
short b[4096][4072];
float a[4096][4072],c[4032][4096],d1[20],d2[20],d[4032][4096];
float xd1[4096],xd2[4096],yd1[4096],yd2[4096],xx[4096],zd1[4096],zd2[4096];
int   n1,n2,n11,n22;
char  osmark[4],over;
int   idate;
float avsky,aslope;
main(ac,av)
int ac; char *av[];
{
  int i,j,k,id;
  int i1,i2;
  float x,x1,x2,x3,x4,y1,y2;
  if(ac<3){
     printf("\n\t*** addfit (overscan,statsky,x-y; last -bias) ******\n");
     printf("\n\tUsage: e4b in_fits bias_fit [!]\n");
     printf("\n\tuse 4 kind of methods to calaulate overscan\n");
     printf("\n\t! plot curve\n");
     exit(0);
  }
  fp=fopen(av[1],"rb");
  if(fp==0){ printf("\n\tinput file not found ! %s\n\n",av[1]); exit(0); }
  sscanf(&av[1][1],"%d",&k);
  idate=0;
  if(k>5700)idate=k;            // 201109   5822-5828
                                // 201110   5859-5866
  n2=4072;   n1=4096;
  fread(h,108,80,fp);
  sscanf(&h[1][28],"%d",&k); 
  if(k!=16){ printf("\n\tinput file must_be *_?.fit(16)\n\n"); exit(0); }
  sscanf(&h[2][28],"%d",&k); 
  if(k!= 2){ printf("\n\tinput file must_be *_?.fit(16)\n\n"); exit(0); }
  k=indexpos(h,"END     ",108);
  if(k==108)fread(&h[k][0],36,80,fp); 
  for(i=0;i<72;i++)for(j=0;j<80;j++)head[i][j]=32;
  strcpy(&head[ 0][0], "SIMPLE  =                    T ");
  strcpy(&head[ 1][0], "BITPIX  =                  -32 ");
  strcpy(&head[ 2][0], "NAXIS   =                    2 ");
  strcpy(&head[ 3][0], "NAXIS1  =                 4096 ");
  strcpy(&head[ 4][0], "NAXIS2  =                 4032 ");
  for(i=0;i<5;i++)head[i][31]='/'; 
  puthead(head, 7,h,"CCD_NO: ");    
  puthead(head, 8,h,"DATE-OBS");
  puthead(head, 9,h,"TIME-OBS");
  puthead(head,10,h,"EXPTIME ");
  puthead(head,11,h,"RA      ");
  puthead(head,12,h,"DEC     ");
  puthead(head,13,h,"HA      ");
  puthead(head,14,h,"EPOCH   ");
  puthead(head,15,h,"FILTER  ");
  puthead(head,16,h,"FOCUSVAL");
  puthead(head,17,h,"OBJECT  ");
  puthead(head,18,h,"INSTRUME");
  puthead(head,19,h,"RDNOISE ");
  puthead(head,20,h,"GAIN    ");
  puthead(head,21,h,"FILENAME");
  puthead(head,22,h,"OBSERVER");
  puthead(head,23,h,"SITEELEV");
  puthead(head,24,h,"CAMTEMP ");
  puthead(head,25,h,"READ_SPD");
  puthead(head,26,h,"TEMPS1  ");
  puthead(head,27,h,"WEATHER0");
  puthead(head,28,h,"WEATHER3");
  puthead(head,30,h,"RDNOCAL ");
  puthead(head,31,h,"GAINCAL ");
  puthead(head,32,h,"SKYADU  ");
  strcpy(&head[71][0], "END"); head[71][3]=32;
  k=indexpos(head,"RA      ",72);
  if(head[5][0]==32){
  for(i=0;i<2;i++)for(j=0;j<80;j++)head[5+i][j]=head[k+i][j];
  head[5][2]=head[6][3]='-';
  head[5][3]=head[6][4]='O';
  head[5][4]=head[6][5]='B';
  head[5][5]=head[6][6]='S';
  }

  sscanf(&head[7][28],"%d",&id);

  strcpy(g,"RDNOIS? ");
  if(id==1 || id==2){
    for(i=0;i<4;i++){
      g[6]=id*4+45+i; k=indexpos(h,g,108); 
      for(j=0;j<3;j++)head[19][12+j+i*5]=h[k][27+j];
    }
  }
  if(id==3){
    k=indexpos(h,"RDNOIS9 ",108); for(j=0;j<3;j++)head[19][12+j]=h[k][27+j];
    k=indexpos(h,"RDNOIS10",108); for(j=0;j<3;j++)head[19][17+j]=h[k][27+j];
    k=indexpos(h,"RDNOIS11",108); for(j=0;j<3;j++)head[19][22+j]=h[k][27+j];
    k=indexpos(h,"RDNOIS12",108); for(j=0;j<3;j++)head[19][27+j]=h[k][27+j];
  }
  if(id==4){
    k=indexpos(h,"RDNOIS13",108); for(j=0;j<3;j++)head[19][12+j]=h[k][27+j];
    k=indexpos(h,"RDNOIS14",108); for(j=0;j<3;j++)head[19][17+j]=h[k][27+j];
    k=indexpos(h,"RDNOIS15",108); for(j=0;j<3;j++)head[19][22+j]=h[k][27+j];
    k=indexpos(h,"RDNOIS16",108); for(j=0;j<3;j++)head[19][27+j]=h[k][27+j];
  }
  k=(id-1)*4+1;
  sprintf(&head[19][31],"/ Read noise (e-) for amp %2d--%2d",k,k+3);
  head[19][63]=32;
//
  strcpy(g,"GAIN?   ");
  if(id==1 || id==2){
    for(i=0;i<4;i++){
      g[4]=id*4+45+i; k=indexpos(h,g,108); 
      for(j=0;j<3;j++)head[20][12+j+i*5]=h[k][27+j];
    }
  }
 if(idate==0){
  if(id==1){
    sprintf(&head[30][11]," 7.1  6.9  7.3 13.7 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"1.79 1.83 1.76 1.82 / IDL results"); head[31][44]=32;
  }
  if(id==2){
    sprintf(&head[30][11]," 9.2  9.6  9.0  9.2 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"2.26 2.17 2.24 2.20 / IDL results"); head[31][44]=32;
  }
  if(id==3){
    k=indexpos(h,"GAIN9   ",108); for(j=0;j<3;j++)head[20][12+j]=h[k][27+j];
    k=indexpos(h,"GAIN10  ",108); for(j=0;j<3;j++)head[20][17+j]=h[k][27+j];
    k=indexpos(h,"GAIN11  ",108); for(j=0;j<3;j++)head[20][22+j]=h[k][27+j];
    k=indexpos(h,"GAIN12  ",108); for(j=0;j<3;j++)head[20][27+j]=h[k][27+j];
    sprintf(&head[30][11],"11.3  9.0 17.4 11.5 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"2.29 2.29 2.26 2.30 / IDL results"); head[31][44]=32;
  }
  if(id==4){
    k=indexpos(h,"GAIN13  ",108); for(j=0;j<3;j++)head[20][12+j]=h[k][27+j];
    k=indexpos(h,"GAIN14  ",108); for(j=0;j<3;j++)head[20][17+j]=h[k][27+j];
    k=indexpos(h,"GAIN15  ",108); for(j=0;j<3;j++)head[20][22+j]=h[k][27+j];
    k=indexpos(h,"GAIN16  ",108); for(j=0;j<3;j++)head[20][27+j]=h[k][27+j];
    sprintf(&head[30][11]," 8.0  8.0  8.0  8.0 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"2.10 2.10 2.10 2.10 / IDL results"); head[31][44]=32;
  }
 }
 if(idate){
  if(id==1){
    sprintf(&head[30][11]," 6.5  7.0  7.0  7.5 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"1.78 1.74 1.55 1.76 / IDL results"); head[31][44]=32;
  }
  if(id==2){
    sprintf(&head[30][11]," 8.1  8.1  8.6  8.6 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"1.84 1.80 1.80 1.81 / IDL results"); head[31][44]=32;
  }
  if(id==3){
    k=indexpos(h,"GAIN9   ",108); for(j=0;j<3;j++)head[20][12+j]=h[k][27+j];
    k=indexpos(h,"GAIN10  ",108); for(j=0;j<3;j++)head[20][17+j]=h[k][27+j];
    k=indexpos(h,"GAIN11  ",108); for(j=0;j<3;j++)head[20][22+j]=h[k][27+j];
    k=indexpos(h,"GAIN12  ",108); for(j=0;j<3;j++)head[20][27+j]=h[k][27+j];
    sprintf(&head[30][11],"10.8  9.3 10.2  9.3 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"1.70 1.70 1.73 1.85 / IDL results"); head[31][44]=32;
  }
  if(id==4){
    k=indexpos(h,"GAIN13  ",108); for(j=0;j<3;j++)head[20][12+j]=h[k][27+j];
    k=indexpos(h,"GAIN14  ",108); for(j=0;j<3;j++)head[20][17+j]=h[k][27+j];
    k=indexpos(h,"GAIN15  ",108); for(j=0;j<3;j++)head[20][22+j]=h[k][27+j];
    k=indexpos(h,"GAIN16  ",108); for(j=0;j<3;j++)head[20][27+j]=h[k][27+j];
    sprintf(&head[30][11],"13.3  7.7  7.1  6.6 / IDL results"); head[30][44]=32;
    sprintf(&head[31][11],"1.65 1.72 1.85 1.89 / IDL results"); head[31][44]=32;
  }
 }
  k=(id-1)*4+1;
  sprintf(&head[20][31],"/ Conversion gain (e-/DN) for amp %2d--%2d",k,k+3);
  head[20][71]=32;
//
  for(j=30;j<48;j++)if(head[22][j]=='/')break; j+=2;
  k=indexpos(h,"OPERATOR",108);
  for(i=0;i<30;i++)head[22][i+j]=h[k][i];
//
  k=indexpos(h,"SITELAT ",108);
  for(i=0;i<23;i++)head[23][i+33]=h[k][i];
  k=indexpos(h,"SITELONG",108);
  for(i=0;i<23;i++)head[23][i+56]=h[k][i];
//
  k=indexpos(h,"WEATHER1",108);
  for(i=0;i<20;i++)head[27][i+29]=h[k][i+10];
  k=indexpos(h,"WEATHER2",108);
  for(i=0;i<21;i++)head[27][i+52]=h[k][i+10];
//
  k=indexpos(h,"WEATHER4",108);
  for(i=0;i<24;i++)head[28][i+29]=h[k][i+10];
  k=indexpos(h,"WEATHER5",108);
  for(i=0;i<20;i++)head[28][i+56]=h[k][i+10];
//
//
  fread(b,n2*n1,2,fp);  fclose(fp);
  swap2(b,n2*n1*2);
  for(i=0;i<n1;i++)for(j=0;j<n2;j++)a[i][j]=b[i][j]+32768.;
  for(i=0;i<n1;i++){
    for(j=0;j<20;j++){ d1[j]=a[i][4032+j]; d2[j]=a[i][4052+j]; }
    xd1[i]=xmedian(d1,20);   xd2[i]=xmedian(d2,20);
  }
// CCD_4 xd2[n11-n1] !=0;
// overscan take median again 
  n11=n1/2; j=99; k=499;  
  for(i=0;i<n1;i++){
    if(i<n11){
      yd1[i]=xmed2(xd1,n11,i,j); yd2[i]=xmed2(xd2,n11,i,j);
      zd1[i]=xmed2(xd1,n11,i,k); zd2[i]=xmed2(xd2,n11,i,k);
    } else {
      yd1[i]=xmed2(&xd1[n11],n11,i-n11,j); yd2[i]=xmed2(&xd2[n11],n11,i-n11,j);
      zd1[i]=xmed2(&xd1[n11],n11,i-n11,k); zd2[i]=xmed2(&xd2[n11],n11,i-n11,k);
    }
  }
  n2=4032; n22=n2/2;                  //  1  3
  k=0;                                //  2  4
// xd~ori; yd~win9; zd~win99;  input
// yd is out, 1: win9 2: win99 3: win999 4 white_black;
  j=0; 
  avsky=0.;         // while mark=4, force balence & slop adjust
// j=0: normal; j=1,4 force mark==4; j=5: force !=4; j=6: force use a value;
//  idate==0
  if(strncmp(av[1],"d5473_0026_1.fit",16)==0)j=3;
  if(strncmp(av[1],"d5508_0087_3.fit",16)==0)j=2;   // force mark=4
  if(strncmp(av[1],"d5508_0089_1.fit",16)==0)j=4;   // force mark=4
  if(strncmp(av[1],"d5534_0098_3.fit",16)==0)insert(yd2,zd2,3200,3800);
  if(strncmp(av[1],"d5535_0132_1.fit",16)==0)j=4;   // force mark=4
// idate==1
  if(strncmp(av[1],"d5824_0015_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5824_0039_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5824_0039_3.fit",16)==0){ avsky=0.49; aslope=0.016; }
  if(strncmp(av[1],"d5824_0057_3.fit",16)==0)j=6;
  if(strncmp(av[1],"d5825_0047_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5825_0051_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5825_0057_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5825_0061_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5825_0067_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5826_0102_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5827_0018_3.fit",16)==0){ avsky=0.80; aslope=0.0; }
  if(strncmp(av[1],"d5827_0092_4.fit",16)==0)j=5;
  if(strncmp(av[1],"d5827_0094_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5827_0096_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5827_0151_3.fit",16)==0)j=6;
  if(strncmp(av[1],"d5828_0054_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5828_0070_3.fit",16)==0)j=6;
// 111024-31    bgright star + black strip
  if(strncmp(av[1],"d5859_0062_3.fit",16)==0)j=6;
  if(strncmp(av[1],"d5859_0141_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5860_0013_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5860_0046_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5860_0052_3.fit",16)==0)j=6;
  if(strncmp(av[1],"d5861_0077_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5861_0109_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5862_0026_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5862_0042_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5862_0090_3.fit",16)==0)j=6;
  if(strncmp(av[1],"d5862_0120_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5862_0121_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5862_0128_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5862_0174_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5864_0058_2.fit",16)==0)j=6;
  if(strncmp(av[1],"d5864_0059_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5864_0121_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5864_0133_3.fit",16)==0)j=6;
  if(strncmp(av[1],"d5864_0168_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5864_0172_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5865_0149_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5865_0165_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5866_0106_4.fit",16)==0)j=6;
  if(strncmp(av[1],"d5866_0173_1.fit",16)==0)j=6;
  if(strncmp(av[1],"d5866_0189_3.fit",16)==0)j=6;
  if(strncmp(av[1],"d5866_0204_2.fit",16)==0)j=6;
   
  osmark[k++]=judge(yd1,zd1,n11,0,0,ac,j,1);
  osmark[k++]=judge(yd2,zd2,n11,0,n22,ac,j,2);
  osmark[k++]=judge(&yd1[n11],&zd1[n11],n11,n11,0,ac,j,3);
  osmark[k++]=judge(&yd2[n11],&zd2[n11],n11,n11,n22,ac,j,4);
// always use win499
  if(idate==0){
    if(j!=1)adjustslope(yd1,n11,1);          //1
    if(j!=2)adjustslope(yd2,n11,1);          //2
    if(j!=3)adjustslope(&yd1[n11],n11,0);    //3
    if(j!=4)adjustslope(&yd2[n11],n11,0);    //4
  }
  for(i=0;i<4;i++){ head[18][26+i]=osmark[i]; printf("%c",osmark[i]); }
  if(ac>3){
    for(i=0;i<n1;i++)xx[i]=i;
    for(i=0;i<n1;i++){ zd1[i]=yd1[i]-50; zd2[i]=yd2[i]-50; }
    i=0; i1=1; i2=2; pgbegin_(&i,"/xw",&i1,&i2,3);
    x=1.5; pgsch_(&x); x1=-60.; x2=4096+60.; 
    y1=(xmedian(zd1,2048)+xmedian(&zd1[2048],2048))/2.-200;  y2=y1+600.;
    pgenv_(&x1,&x2,&y1,&y2,&i,&i);
      pglabel_("","",av[1],0,0,strlen(av[1]));
    i2=3; pgsci_(&i2); pgpoint_(&n1,xx,xd1,&i1);
    i2=1; pgsci_(&i2); pgpoint_(&n1,xx,zd1,&i1);
    x=y2-100; y1=100; y2=y1+n11; 
    pgtext_(&y1,&x,"gate 1",6); pgtext_(&y2,&x,"gate 3",6);
    y1=xmedian(zd2,2048); y2=xmedian(&zd2[2048],2048);  
      if(id!=4)y1=(y1+y2)/2.-200.; else y1=y2-200.;          y2=y1+600.;
    pgenv_(&x1,&x2,&y1,&y2,&i,&i);
    i2=3; pgsci_(&i2); pgpoint_(&n1,xx,xd2,&i1);
    i2=1; pgsci_(&i2); pgpoint_(&n1,xx,zd2,&i1);
    x=y2-100; y1=100; y2=y1+n11; 
    pgtext_(&y1,&x,"gate 2",6); pgtext_(&y2,&x,"gate 4",6);
    pgend_();
  }
// some 111024~111031 *_4.fit, last 4 gate, need not modify overscan
  if(idate>=5859 && idate<=5866){   // 111024~~111031
    j=0;                            // no bright star, but strip
    if(strncmp(av[1],"d5859_0062_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5859_0154_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5859_0160_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5860_0067_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0068_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0073_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0075_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0081_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0085_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0100_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0101_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0104_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0105_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0106_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5860_0115_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5862_0170_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5863_0130_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5863_0139_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5863_0153_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5864_0110_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5864_0179_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5864_0181_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0049_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0051_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0061_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0066_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0073_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0084_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0085_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0089_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0092_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0098_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0114_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0136_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0138_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0149_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0155_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0168_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0169_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5865_0180_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0038_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0077_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0100_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0101_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0104_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0108_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0111_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0158_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0165_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5866_0166_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5866_0167_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5866_0168_4.fit",16)==0)j=3;
    if(strncmp(av[1],"d5866_0173_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0174_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0175_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0214_4.fit",16)==0)j=4;
    if(strncmp(av[1],"d5866_0218_4.fit",16)==0)j=4;
    if(j==3)for(i=2048;i<n1;i++)yd1[i]=xd1[i];
    if(j==4)for(i=2048;i<n1;i++)yd2[i]=xd2[i];
  }
// 111024~111031  some 1/4 ccd need add a value 
  if(idate>=5859 && idate<=5866){   // 111024~~111031
    j=k=0;
    if(strncmp(av[1],"d5862_0104_1.fit",16)==0){ j=4; x4=2.; } // force bright
    if(strncmp(av[1],"d5863_0065_2.fit",16)==0){ j=4; x4=6.; } // force bright
    if(strncmp(av[1],"d5863_0131_3.fit",16)==0){ j=2; x2=6.; } // force bright
    if(strncmp(av[1],"d5864_0022_1.fit",16)==0){ j=4; x4=3.; } // force bright
    if(strncmp(av[1],"d5864_0022_1.fit",16)==0){ k=1; x1=2.; } // force bright
    if(strncmp(av[1],"d5864_0048_3.fit",16)==0){ k=1; x1=4.; } // force bright
    if(strncmp(av[1],"d5864_0132_4.fit",16)==0){ k=1; x1=6.; } // force bright
    if(strncmp(av[1],"d5865_0023_2.fit",16)==0){ k=1; x1=-2.; } 
    if(strncmp(av[1],"d5865_0023_2.fit",16)==0){ j=4; x4=-5.; } 
    if(strncmp(av[1],"d5866_0037_2.fit",16)==0){ j=4; x4=-4.; } 
    if(strncmp(av[1],"d5866_0051_2.fit",16)==0){ j=4; x4=-4.; } 
    if(strncmp(av[1],"d5866_0055_2.fit",16)==0){ j=4; x4=-3.; } 
    if(strncmp(av[1],"d5866_0056_2.fit",16)==0){ j=4; x4=-3.; } 
    if(k==1)for(i=0; i<n1/2;i++)yd1[i]-=x1;
    if(k==3)for(i=2048;i<n1;i++)yd1[i]-=x3;
    if(j==2)for(i=0; i<n1/2;i++)yd2[i]-=x2;
    if(j==4)for(i=2048;i<n1;i++)yd2[i]-=x4;
  }
// overscan correct
  for(i=0;i<n1;i++){
    for(j=0;j<2016;j++)c[4031-j][i]=a[i][j]-yd1[i];
    for(j=2016;j<4032;j++)c[4031-j][i]=a[i][j]-yd2[i];
  }
// -bias
  fp=fopen(av[2],"rb");
  if(fp==0){ printf("\n\t%s not found!\n",av[2]); exit(0); }
  fseek(fp,72*80,0);
  fread(d,n2*n1,4,fp);
  fclose(fp);
  swap4(d,n2*n1*4);
  for(i=0;i<n2;i++)for(j=0;j<n1;j++)c[i][j]-=d[i][j];
// -bias end
  stat4sky(c,&x,avsky); 
  k=indexpos(head,"SKYADU  ",72);
  sprintf(&head[k][9]," %20.2f / -bias only",x); head[k][43]=32;
  swap4(c,n2*n1*4);
  av[1][0]='e';
  fp=fopen(av[1],"wb");
  fwrite(head,72,80,fp);
  fwrite(c,n2*n1,4,fp);
  fclose(fp);
  printf("  %s  produced !\n",av[1]);
}

adjustslope(a,n,m)                  // a[n], if(m)reverise
float *a; int n,m;
/*
#include <stdio.h>
#define n 24
main()
{
  int i,j,k,a[n]={3,2,2,3,4,4,5,6,7,8,9,8,7,6,6,7,7,8,9,11,14,13,14,15};
  for(i=0;i<n;i++)printf("%d ",a[i]); printf("\n");
  k=a[0]; for(i=1;i<n;i++)if(k>a[i])k=a[i];  a[0]=k; 
  for(i=1;i<n;i++){ if(a[i]>=a[i-1]) continue; 
    k=a[i]; for(j=1;j<i;j++)if(a[j]>k)a[j]=k;
  }
  for(i=0;i<n;i++)printf("%d ",a[i]); printf("\n");
}
*/
{
  int i,j;
  float x;
  if(m)for(i=0;i<n/2;i++){ x=a[i]; a[i]=a[n-1-i]; a[n-1-i]=x; }
  x=a[0]; for(i=1;i<n;i++)if(x>a[i])x=a[i];  a[0]=x; 
  for(i=1;i<n;i++){ if(a[i]>=a[i-1]) continue; 
    x=a[i]; for(j=1;j<i;j++)if(a[j]>x)a[j]=x;
  }
  if(m)for(i=0;i<n/2;i++){ x=a[i]; a[i]=a[n-1-i]; a[n-1-i]=x; }
}

insert(a,b,k1,k2)
float *a,*b; int k1,k2;
{
  float y1,y2; int i;
  y1=a[k1]; y2=a[k2];
  for(i=k1+1;i<k2;i++)a[i]=b[i]=y2-(y2-y1)*(k2-i)/(k2-k1);
}

judge(yd,zd,n,n11,n22,ac,f4,gate)           // input yd-win99; zd~win499
float *yd,*zd; int n,n11,n22,ac,f4,gate;    // output yd
{
  int i,i1,i2,j,j1,j2,k;
  float x,s[2048],e[2048];
  short f[2016];
  if(f4==gate)goto lf4;
  j=999;
  for(i=0;i<n;i++)e[i]=xmed2(yd,n,i,j);  
  x=(e[512]+e[1024]+e[1536])/3.;
  if(f4==6)for(i=0;i<n;i++)zd[i]=x;  // force a value
  j1=j2=0;
  for(i=99;i<n-99;i++){
    i1=yd[i]-zd[i]; if(i1<0)i1=-i1; if(i1>2)j1+=(i1-2)*(i1-2);
    i2=zd[i]-e[i];  if(i2<0)i2=-i2; if(i2>1)j2+=(i2-1)*(i2-1);
  }
  if(ac>3)printf(" %d_%d ",j1,j2);
  for(i=0;i<n;i++)yd[i]=zd[i];      // always output, use 499
  if(j1==0) return('1');
  if(j2<=1) return('2'); 
  j=0; for(i=0;i<n-1;i++)j+=(e[i]-e[i+1])*(e[i]-e[i+1]);
  k=0; for(i=0;i<n-1;i++)if(e[i]!=e[i+1])k++;
  if(j==k)  return('3'); 
     
  if(f4==5 || f4==6)  return('3');    //   force return;
lf4:
  for(i=0;i<n;i++){
    for(j=0;j<2016;j++)e[j]=a[i+n11][j+n22];
    k=0; for(j=0;j<2016;j++)if(e[j]<32000.)f[k++]=e[j];
    white_black(f,1,k,&s[i],&x);
  }                                         // header find balance in 20 
  for(i=0;i<n;i++)zd[i]=xmed2(s,n,i,9); 
  for(i=0;i<n;i++){ e[i]=yd[i]; s[i]=zd[i]; }
  sortx1(e,n); sortx1(s,n);
  x=xmedian(e,20)-xmedian(s,20);
  for(i=0;i<n;i++)yd[i]=zd[i]+x;
          return('4');
}

stat4sky(c,sky,avsky)
float c[][4096],*sky, avsky;
{
  int i,j,k,l=0,n11,n22;
  float sk[4],s;
  short d[200*200];
  k=0; for(i=2016;i<2416;i+=2)for(j=1648;j<2048;j+=2)d[k++]=c[i][j]/2.;  //1
  white_black(d,1,k,&sk[l],&s);  sk[l++]*=2;
  k=0; for(i=1616;i<2016;i+=2)for(j=1648;j<2048;j+=2)d[k++]=c[i][j]/2.;  //2
  white_black(d,1,k,&sk[l],&s);  sk[l++]*=2;
  k=0; for(i=2016;i<2416;i+=2)for(j=2048;j<2448;j+=2)d[k++]=c[i][j]/2.;  //3
  white_black(d,1,k,&sk[l],&s);  sk[l++]*=2;
  k=0; for(i=1616;i<2016;i+=2)for(j=2048;j<2448;j+=2)d[k++]=c[i][j]/2.;  //4
  white_black(d,1,k,&sk[l],&s);  sk[l++]*=2;
  printf(" block's skyadu: %9.2f %9.2f %9.2f %9.2f\n",sk[0],sk[1],sk[2],sk[3]);
  s=0.; j=0;
  for(i=0;i<4;i++)if(osmark[i]!='4'){ s+=sk[i]; j++; }
  *sky=s/j;
  if(sk[1]==0.){ *sky=sk[3]; return; }                // for idate==0 CCD_4
  for(l=0;l<4;l++)if(osmark[l]=='4'){
    s=*sky-sk[l]; n11=n22=0;         // 0~2, 1~3 background as same
    if(l==0){ n22=2016;          if(osmark[1]!='4')s=sk[1]-sk[0]; }
    if(l==1)                     if(osmark[0]!='4')s=sk[0]-sk[1];
    if(l==2){ n22=2016; n11=2048;if(osmark[3]!='4')s=sk[3]-sk[2]; } 
    if(l==3){ n11=2048;          if(osmark[2]!='4')s=sk[2]-sk[3]; }
    if(avsky!=0.){
      s*=avsky;
      printf("marked 4, adjust a little*****%f\n",s);
      for(j=0;j<2016;j++)for(i=0;i<2048;i++)c[j+n22][i+n11]+=s-i*aslope;
    }
  }
}
/*
#include <stdio.h>
main(int ac, char **av)
{
  int  i,j,j1,j2,j3,k,c,n;
  if(ac<4){ printf("\n\t***** a dimension window minw********\n"); exit(0); }
  sscanf(av[1],"%d",&n); sscanf(av[2],"%d",&c); sscanf(av[3],"%d",&j3);
  for(i=0;i<n;i++){
    j1=j2=c;
    k=2*i+1;         if(k<c)j1=k; 
    k=(n-i)*2-1;     if(k<c)j2=k;
    j=(j1>j2)?j2:j1; if(j<j3)j=j3;
    printf("%2d %2d     %2d\n",i,j,c);
  }
}
*/
xmed2(b,n,i1,c)
float *b; int n,i1,c;
{
  int i,j,k,i2;
  float a[999],x;                     // alway in odd
// adjust window
  int j1,j2;
  j1=j2=c;
  k=2*i1+1;     if(k<c)j1=k;
  k=(n-i1)*2-1; if(k<c)j2=k;
  j=(j1>j2)?j2:j1; if(j<9)j=9;        // assume min_window is 9
// adjust window end
  k=j/2;
  i2=i1-k; if(i2<0)i2=0;
  if(i1+k>=n-1)i2=n-j-1;
  for(i=0;i<j;i++)a[i]=b[i2+i]; 
  sortx1(a,j);
  return(a[k]);
}

xmedian(b,n)
float *b; int n;
{
  int i,j;
  float x,a[4096];
  for(i=0;i<n;i++)a[i]=b[i];
  if(n<=20){
    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;
    }
  } else sortx1(a,n);
  i=n%2; j=n/2; 
  if(i)return(a[j]);
  return((a[j]+a[j-1])/2);
}

sortx1(a,n)
float *a; int n;
{
  int i,j,k=2,ii,ifin;
  float x;
  while(k<n)k*=2;
  k=(3*k)/4-1; if(n<k)k=n;
l20:
  k/=2;        ifin=n-k;
  for(ii=0;ii<ifin;ii++){
    i=ii;      j=i+k;
    if(a[i]<=a[j])continue;
    x = a[j]; 
l40:
    a[j]= a[i]; 
    j=i;       i-=k;
    if(i<0)goto l60;
    if(a[i]>x)goto l40;
l60:
    a[j]=x;  
  }
  if(k>1)goto l20;
}

puthead(head,k,h,f)
char head[][80],h[][80],*f;
int  k;
{
  int i,j;
  j=indexpos(h,f,140);
  if(j==140){ sprintf(&head[k][0],"%s=",f); head[k][9]=32; }
  else for(i=0;i<80;i++)head[k][i]=h[j][i];
}