#include <stdio.h>
#include <math.h>
#define gap 30.                // telescape move pixel, 1st_time move 2*gap
float c[4096][4072];
short b[4096][4072];
char  head[72][80];
FILE  *fp2,*fp3;
char a[130];
int  n1=4072,n2=4096;
main(ac,av)
int ac; char *av[];
{
  int i,j,k;
  int ix,iy;
  float x,y,sky,seeing,r;
  char f1[60],f2[60],f3[60];
  if(ac<2){ printf("\n\t ********get center_ADU sky seeing ring_R *******\n");
            printf("\n\t before run: d4 focus.0067.fits  (divide to 4 file)");
            printf("\n\t before run: d4dat focus.0067_?  (produce *_?.dat)\n");
            printf("\n\t Ex: d4f focus.0067_?\n\n");
            exit(0);
          }
  strcpy(f1,av[1]);
  strcpy(f2,f1); strcpy(f3,f1); 
  strcat(f2,".fit");  strcat(f3,".dat"); 
  fp2=fopen(f2,"rb"); if(fp2==0){ printf("%s not found!\n",f2); exit(0); }
  fp3=fopen(f3,"r");  if(fp3==0){ printf("%s not found!\n",f3); exit(0); }
  fread(head,80,72,fp2);   
  sscanf(&head[1][24],"%d",&k);
  if(k!=16){ printf("%s not a integer_fits file\n",f2); exit(0); }
  sscanf(&head[3][24],"%d",&i);
  sscanf(&head[4][24],"%d",&j);
  if(i!=n1 || j!=n2){ printf("NAXIS err!\n"); exit(0); }
  k=indexpos(head,"END     ",72);
  if(k==72){
l00:
    fread(head,36,80,fp2);
    k=indexpos(head,"END     ",36);
    if(k==36)goto l00;
  }
  fread(b,n1*n2,2,fp2);
  fclose(fp2);
  swap2(b,n1*n2*2);
  for(j=0;j<n2;j++)for(i=0;i<n1;i++)c[j][i]=b[j][i]+32768.;

l10:
  fgets(a,130,fp3); if(feof(fp3))goto l20;
  printf("#%s",a);
  if(a[0]=='#')goto l10;
  else printf("# point_center(X,Y)  center_ADU  Sky     Seeing   ring_R\n");
  sscanf(a,"%d %d",&ix,&iy);
  getsky(ix,iy,&sky);
  center(--ix,--iy,&x,&y,sky);       // 1st time    
  if(y<15. ||y >4080.)goto l10;
  if(x<15. ||x >4017.)goto l10;
  see(x,y,sky,&seeing);
  ring(x,y,&r);
  printf("%9.3f %9.3f %8.0f %8.1f %8.2f %8.2f\n",
         x+.5,y+.5,c[iy][ix],sky,seeing,r);
  ix=x+0.5; iy=y-gap*2;
            if(y>n2/2)iy=y+gap*2;
  for(i=0;i<6;i++){
    getsky(ix,iy,&sky);
    center(ix,iy,&x,&y,sky);
    if(x<15. || x >4017.)break;
    if(y<15. || y >4080.)break;
    see(x,y,sky,&seeing);
    ring(x,y,&r);
  printf("%9.3f %9.3f %8.0f %8.1f %8.2f %8.2f\n",
         x+.5,y+.5,c[iy][ix],sky,seeing,r);
    ix=x+0.5; iy=y-gap;
            if(y>n2/2)iy=y+gap;
    if(iy<15 || iy >4080)break;
    if(ix<15 || ix >4017)break;
  } 
  goto l10;
l20:
  fclose(fp3);  
}

getsky(ix,iy,sky)
int ix,iy; float *sky;
{
  int i1,i2,j1,j2, i,j,k;
  short d[300*300];
  float s;
  i1=ix-150; i2=ix+150; j1=iy-150; j2=iy+150;
  if(i1<0)i1=0; if(i2>4063)i2=4063;
  if(j1<0)j1=0; if(j2>4063)j2=4063;
  k=0; for(i=i1;i<i2;i++)for(j=j1;j<j2;j++)d[k++]=c[j][i];
  white_black(d,1,k,sky,&s);
}

center(ix,iy,x,y,sky)
int ix,iy; float *x,*y,sky;
{
  float b[20],s,z;
  int i,j,k,i5=9;
  k=0;
  for(i=ix-i5;i<=ix+i5;i++){
    z=0.;
    for(j=iy-i5;j<=iy+i5;j++)z+=c[j][i]-sky;
    b[k++]=z;
  }
  z=s=0; for(i=0;i<k;i++){ z+=b[i]; s+=b[i]*(i-i5); }
  s/=z; *x=ix+s;
  k=0;
  for(j=iy-i5;j<=iy+i5;j++){
    z=0.;
    for(i=ix-i5;i<=ix+i5;i++)z+=c[j][i]-sky;
    b[k++]=z;
  }
  z=s=0; for(i=0;i<k;i++){ z+=b[i]; s+=b[i]*(i-i5); }
  s/=z; *y=iy+s;
}

see(x,y,white,seeing)
float x,y,white,*seeing; 
// (a81+a84)/2 *180*3600/pi = arcsec / pixel
// bok is 0.45"/pixel
{
  int i,j,k,ix,iy, b[113];   // (29-1)*4+1
  float z,s1,s2;
  ix=x+0.5; iy=y+0.5;
  for(j=0,i=ix-14;i<=ix+14;i++,j+=4){
    z=0.;
    for(k=-2;k<=2;k++)z+=c[iy+k][i]-white;
    b[j]=z;
  }
  for(j=0;j<112;j+=4)b[j+2]=(b[j]+b[j+4])/2.+0.5;
  for(j=0;j<112;j+=2)b[j+1]=(b[j]+b[j+2])/2.+0.5;
  histat(b,&z,&s1,112);
  for(j=0,i=iy-14;i<=iy+14;i++,j+=4){
    z=0.;
    for(k=-2;k<=2;k++)z+=c[i][ix+k]-white;
    b[j]=z;
  }
  for(j=0;j<112;j+=4)b[j+2]=(b[j]+b[j+4])/2.+0.5;
  for(j=0;j<112;j+=2)b[j+1]=(b[j]+b[j+2])/2.+0.5;
  histat(b,&z,&s2,112);
  *seeing=(s1+s2)/8.*2.355*0.45;
  return;   
}

ring(x,y,r)
float x,y,*r;
{
  int i,ix,iy,iz1,iz2,iz3,iz4,jz1,jz2,jz3,jz4;
  float z;
  ix=y+0.5;  iy=x+0.5;
  z=c[ix][iy]; iz1=0; for(i=1;i<12;i++)if(c[ix+i][iy]>z){iz1=i; z=c[ix+i][iy];}
  z=c[ix][iy]; iz2=0; for(i=1;i<12;i++)if(c[ix-i][iy]>z){iz2=i; z=c[ix-i][iy];}
  z=c[ix][iy]; iz3=0; for(i=1;i<12;i++)if(c[ix][iy+i]>z){iz3=i; z=c[ix][iy+i];}
  z=c[ix][iy]; iz4=0; for(i=1;i<12;i++)if(c[ix][iy-i]>z){iz4=i; z=c[ix][iy-i];}
  z=c[ix][iy];jz1=0;for(i=1;i<9;i++)if(c[ix+i][iy+i]>z){jz1=i;z=c[ix+i][iy+i];}
  z=c[ix][iy];jz2=0;for(i=1;i<9;i++)if(c[ix+i][iy-i]>z){jz2=i;z=c[ix+i][iy-i];}
  z=c[ix][iy];jz3=0;for(i=1;i<9;i++)if(c[ix-i][iy+i]>z){jz3=i;z=c[ix-i][iy+i];}
  z=c[ix][iy];jz4=0;for(i=1;i<9;i++)if(c[ix-i][iy-i]>z){jz4=i;z=c[ix-i][iy-i];}
//  printf("%d %d %d %d %d %d %d %d\n",iz1,iz2,iz3,iz4,jz1,jz2,jz3,jz4);
  *r=((iz1+iz2+iz3+iz4)+(jz1+jz2+jz3+jz4)*1.4142)/8.;
  return;
}