// divided 4k_pixel to 4*4, add 5*5 to a pixel of 2k_fit
// in order to NO lost any adu
#include <stdio.h>
#include <math.h>
#define n 4096
#define m 2048
FILE *fp;
float a[n*4][n*4],b[n][n],c[m][m];
char head[72][80];
main(int ac, char **av)
{
  int i,j,k,i1,j1,ii,jj, nxc,nyc,nx,ny;
  float x;
  double w,z,pi,cx,cy;
  char d[30];
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  if(ac<3){
    printf("\n\t ****** z4kto2k, keep ADU value ********");
    printf("\n\t                               2011,2,15\n");
    printf("\n\t Usage: infile(4k) outfile(2k) [2048 2048]\n");
    printf("\n\t        1280 <  IX,IY  < 2816\n\n");
    exit(0);
  }
  nxc=nyc=2048;
  if(ac>3)sscanf(av[3],"%d",&nxc);
  if(ac>4)sscanf(av[4],"%d",&nyc);
  if(nxc<1280 || nyc<1280 || nxc>2816 || nyc>2816){
    printf("\nerror: IX,IY outof ccd frame\n\n"); exit(0); 
  }
  fp=fopen(av[1],"rb"); if(fp==0){ printf("\ninfile not found!\n\n"); exit(0); }
  fread(head,72,80,fp); fread(b,n*n,4,fp); fclose(fp); swap4(b,n*n*4);
  for(j=0;j<n;j++)for(i=0;i<n;i++){
    x=b[j][i]/16.;  i1=i*4; j1=j*4;
    for(jj=j1;jj<j1+4;jj++)for(ii=i1;ii<i1+4;ii++)a[jj][ii]=x;
  }
  nx=(nxc-1280)*4; ny=(nyc-1280)*4;
  for(j=0;j<m;j++)for(i=0;i<m;i++){
    i1=i*5+nx; j1=j*5+ny; x=0.;
    for(jj=j1;jj<j1+5;jj++)for(ii=i1;ii<i1+5;ii++)x+=a[jj][ii];
    c[j][i]=x;
  } k=m;
  sprintf(&head[0][64],"from 4k*4k"); head[0][74]=32;
  sprintf(&head[3][24],"%6d",k);      head[3][30]=32;
  sprintf(&head[4][24],"%6d",k);      head[4][30]=32;
  k=indexpos(head,"VOLT2   ",72);
  if(k!=72){
    sscanf(&head[k][20],"%f",&x); x/=0.64;
    sprintf(&head[k][20],"%10.2f",x);   head[k][30]=32;
  }
  k=indexpos(head,"SKYADU  ",72);
  if(k!=72){
    sscanf(&head[k][20],"%f",&x); x/=0.64;
    sprintf(&head[k][20],"%10.2f",x);   head[k][30]=32;
  }
  k=indexpos(head,"INSTRUME",72);
  if(k!=72){
    for(i=10;i<30;i++)if(head[k][i]<='9' && head[k][i]>='0')break;
    head[k][i]='3';                   // assume sort,  4k --> 2k
  }
  k=indexpos(head,"A81     ",72);
  if(k!=72){
    for(i=0;i<6;i++){
      sscanf(&head[k][10],"%le",&z); if(i<4)z/=0.8; else z/=1.6;
      sprintf(&head[k][10],"%20.12le",z);   head[k++][30]=32;
    }
    sscanf(&head[k][10],"%lf",&w); 
    sscanf(&head[k+1][10],"%lf",&z); 
//     printf("%10.7lf %10.7lf\n",w,z);
    w+=(2048-nxc)*1.32/3600.*cy/cos(z)-0.000004;
    z+=(2048-nyc)*1.44/3600.*cy       -0.000004;
//     printf("%10.7lf %10.7lf\n",w,z);
    sprintf(&head[k][10],"%20.7lf",w);   head[k++][30]=32;
    sprintf(&head[k][10],"%20.7lf",z);   head[k++][30]=32;
  }
  swap4(c,m*m*4);
  fp=fopen(av[2],"wb"); fwrite(head,72,80,fp); fwrite(c,m*m,4,fp); fclose(fp);
  printf("\nproduced:  %s  OK!\n\n",av[2]);
  k=indexpos(head,"A81     ",72);
  if(k!=72){ sprintf(d,"coord9 %s !",av[2]); system(d); } 
}