#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;
main(ac,av)
int ac; char *av[];
{
  int i,j,k,id;
  int i1,i2;
  float x,x1,x2,y1,y2;
  if(ac<3){
     printf("\n\t*** addfit (overscan,statsky,x-y; last -bias) ******\n");
     printf("\n\tUsage: e4b in_fits bias_fit [ win ] [n]\n");
     printf("\n\tuse 4 kind of methods to calaulate overscan\n");
     printf("\n\tdefault median_win=9, n plot curve\n");
     printf("\n\tif n=1~4, force smooth overscan_curve 1 grade\n\n");
     exit(0);
  }
  fp=fopen(av[1],"rb");
  if(fp==0){ printf("\n\tinput file not found ! %s\n\n",av[1]); exit(0); }
  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(&head[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(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;
  }
  k=(id-1)*4+1;
  sprintf(&head[20][31],"/ Conversion gain (e-/DN) for amp %2d--%2d",k,k+3);
  head[20][71]=32;
//
  k=indexpos(h,"OPERATOR",108);
  for(i=0;i<30;i++)head[22][i+33]=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 
  j=9;   // median_window size
  if(ac>3)sscanf(av[3],"%d",&j);
  printf(" win=%d ",j);
  n11=n1/2;
  k=99;
  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;
  if(ac>4)over=av[4][0]-48; 
  j=0; if(over==1)j=1; osmark[k++]=judge(yd1,zd1,n11,0,0,j);
  j=0; if(over==2)j=1; osmark[k++]=judge(yd2,zd2,n11,0,n22,j);
  j=0; if(over==3)j=1; osmark[k++]=judge(&yd1[n11],&zd1[n11],n11,n11,0,j);
  j=0; if(over==4)j=1; osmark[k++]=judge(&yd2[n11],&zd2[n11],n11,n11,n22,j);
  for(i=0;i<4;i++)head[18][26+i]=osmark[i];
  if(ac>4){
    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_();
  }
// 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); 
  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]);
}

judge(yd,zd,n,n11,n22,os)           // input yd-win9; zd~win99
float *yd,*zd; int n,n11,n22,os;   // output yd
{
  int i,i1,i2,j,k,o;
  float x,y,z[2048],e[2048];
  short f[2016];
  o=os;
  for(i=0;i<n;i++)z[i]=yd[i]-zd[i]; 
  j=0; for(i=0;i<n;i++)if(z[i]>9.)j++;
  if(j==0){ if(o==0){ printf("1"); return('1'); }
            else o=0; }            // normal mark=1
  j=0; for(i=0;i<n-1;i++)j+=(zd[i]-zd[i+1])*(zd[i]-zd[i+1]);
  k=0;  for(i=0;i<n-1;i++)if(zd[i]!=zd[i+1])k++;
  if(j==k){ if(o==0){ for(i=0;i<n;i++)yd[i]=zd[i]; printf("2"); return(50); }
            else o=0; }            //media;  2
  j=999;
  for(i=0;i<n;i++)e[i]=xmed2(zd,n,i,j);  
  for(i=0;i<n;i++)zd[i]=e[i];
  j=0; for(i=0;i<n-1;i++)j+=(zd[i]-zd[i+1])*(zd[i]-zd[i+1]);
  k=0;  for(i=0;i<n-1;i++)if(zd[i]!=zd[i+1])k++;
  if(j==k){ if(o==0){ for(i=0;i<n;i++)yd[i]=zd[i]; printf("3"); return(51); }
            else o=0; }            // media;  3
  for(i=0;i<n;i++){
    for(j=0;j<2016;j++)e[j]=a[i+n11][j+n22];
//    sortx1(e,2016); z[i]=e[n22/4];           // avoid bright star
    k=0; for(j=0;j<2016;j++)if(e[j]<32000.)f[k++]=e[j];
    white_black(f,1,k,&z[i],&x);
  }
  j=9;      printf("4");                  //  background   mark=4
  for(i=0;i<n;i++)zd[i]=xmed2(z,n,i,j); 
//  x,y,e[],z[]
  for(i=0;i<n;i++){ e[i]=yd[i]; z[i]=zd[i]; }
  sortx1(e,n); sortx1(z,n);
  x=xmedian(e,20)-xmedian(z,20);
  for(i=0;i<n;i++)yd[i]=zd[i]+x;
  return(52);
}

stat4sky(c,sky)
float c[][4096],*sky;
{
  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 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]; }
    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;
  }
}
/*
#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,108);
  if(j==108){ sprintf(&head[k][0],"%s=",f); head[k][9]=32; }
  else for(i=0;i<80;i++)head[k][i]=h[j][i];
}