#include <stdio.h>
FILE *fp;
char  h[144][80],head[72][80],f1[30],g[9];
short b[4096][4072];
float a[4096][4072],c[4032][4096],d1[20],d2[20],xd1,xd2;
int   n1,n2;
main(ac,av)
int ac; char *av[];
{
  int i,j,k,id;
  float x;
  if(ac<2){
     printf("\n\t*** addfit (overscan,statsky,x-y,) ******\n");
     printf("\n\tUsage: e4 in_fits\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); }
  n1=4072;   n2=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,n1*n2,2,fp);  fclose(fp);
  swap2(b,n1*n2*2);
  for(i=0;i<n2;i++)for(j=0;j<n1;j++)a[i][j]=b[i][j]+32768.;
  for(i=0;i<n2;i++){
    for(j=0;j<20;j++){ d1[j]=a[i][4032+j]; d2[j]=a[i][4052+j]; }
    xd1=xmedian(d1,20);   xd2=xmedian(d2,20);
    for(j=0;j<2016;j++)c[4031-j][i]=a[i][j]-xd1;
    for(j=2016;j<4032;j++)c[4031-j][i]=a[i][j]-xd2;
  }
  n1=4032;
  k=indexpos(head,"SKYADU  ",72);
  stat4sky(c,&xd1); 
  sprintf(&head[k][9]," %20.2f ",xd1); head[k][31]='/';
  swap4(c,n1*n2*4);
  av[1][0]='e';
  fp=fopen(av[1],"wb");
  fwrite(head,72,80,fp);
  fwrite(c,n1*n2,4,fp);
  fclose(fp);
  printf("%s  produced !\n",av[1]);
}

stat4sky(c,sky)
float c[][4096],*sky;
{
  int i,j,k;
  float sky1,sky2,sky3,sky4,s;
  short d[100*100];
  k=0; for(i=1816;i<1916;i++)for(j=1848;j<1948;j++)d[k++]=c[i][j]/2.;  //1
  white_black(d,1,k,&sky1,&s);  sky1*=2;
  k=0; for(i=1816;i<1916;i++)for(j=2148;j<2248;j++)d[k++]=c[i][j]/2.;  //2
  white_black(d,1,k,&sky2,&s);  sky2*=2;
  k=0; for(i=2116;i<2216;i++)for(j=1848;j<1948;j++)d[k++]=c[i][j]/2.;  //3
  white_black(d,1,k,&sky3,&s);  sky3*=2;
  k=0; for(i=2116;i<2216;i++)for(j=2148;j<2248;j++)d[k++]=c[i][j]/2.;  //4
  white_black(d,1,k,&sky4,&s);  sky4*=2;
  printf("4 block's skyadu: %9.2f %9.2f %9.2f %9.2f\n",sky1,sky2,sky3,sky4);
  *sky=(sky1+sky2+sky3+sky4)/4.;
 if(sky1==0.)*sky=sky2;                 // ???
}

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

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];
}