#include <stdio.h>
#include <math.h>
// SL array
#define ns  20000           // SL all data in 31'
double sx[ns],sy[ns];
int    si1[ns],si2[ns];
float  sf1[ns],sf2[ns],sf3[ns],sf4[ns];
// our
#define os   5000
double ox[os],oy[os];
float  o[os][13];
short  om[os],kk[os];        // kk[] keep SL_num
short  bp[os];
// common
#define sig2   1.0           // dx*dx+dy*dy
double cx,cy,pi,dc;
char c1[14],c2[14];
FILE *fp,*fp1,*fp2;
#define n1  4096
#define n2  4032
char bad1[n2][n1],bad2[n2][n1],bad3[n2][n1],bad4[n2][n1];
main(int ac,char **av)
{
  int    i,j,k,l,m,n,ip;
  double xx,yy,zz,sig;
  float  dx,dy,zx,zy,t0,t1,f0;
  char   a[80],b[128],f1[20],f2[20],c,d[180];    //d_164
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  if(ac<2){ printf("\n\t*** add SL_u(psf,model) to our cat ****\n");
            printf("\n\t                   jiang 2011,4,13\n");
            printf("\n\tusl_chart: field 31'  mag_limit 24m\n");
            printf("\n\tp*.cat ---> q*.cat (add bad_pixel mark)\n");
            printf("\n\tUsage: um ! [24.]\n\n");
            exit(0);
  }
// read bad_pixel table
  readbad(bad1,1); readbad(bad2,2); readbad(bad3,3); readbad(bad4,4);
// NOTICE  x-1,y-1 == bad[y][x];
/*
                 *
               * * *
             * * o * *     o is star center
               * * *       if meet *, regard bad
                 *
    i=bad3[1][52]; j=bad3[1][53]; k=bad3[1][54];
    printf("%d %d %d\n",i,j,k);
*/
  f0=24.; if(ac>2)sscanf(av[2],"%f",&f0);
  system("ls p*.cat >um.tmp");
  fp1=fopen("um.tmp","r");
l00:
  fgets(f1,20,fp1); if(feof(fp1)) exit(0);
  k=strlen(f1); f1[k-1]=0;
// our
  n=0;
  fp=fopen(f1,"r");
  strcpy(f2,f1); f2[0]='q';  fp2=fopen(f2,"w");
  fgets(b,128,fp);   fprintf(fp2,"%s",b);
  fgets(b,128,fp);   fprintf(fp2,"%s",b);  
  chs(&b[67],&dc); dc=cos(dc*cy)*54000.; 
  b[79]=0; sprintf(a,"usl %s 31 um.cat %4.1f",&b[54],f0);
  system(a);
  for(i=0;i<9;i++){ fgets(b,128,fp);   fprintf(fp2,"%s",b); }
  fgets(b,128,fp);  fprintf(fp2,"# 1 bad_pixel     c1\n"); 
  for(i=2;i<=17;i++){ fgets(b,128,fp);   fprintf(fp2,"%s",b); }
  fprintf(fp2,"#18 SL_number     c20\n");
  fprintf(fp2,"#19 SL_u_psf      f7.3\n");
  fprintf(fp2,"#20 SL_u_psf_err  f6.3\n");
  fprintf(fp2,"#21 SL_u_model    f7.3\n");
  fprintf(fp2,"#22 SL_u_model_er f6.3\n");

l30:
  fgets(b,128,fp); if(feof(fp))goto l40;
  if(b[0]=='#')fprintf(fp2,"%s",b);
  if(b[0]=='#' || b[82]=='9')goto l30;
  chs(&b[5],&ox[n]); chs(&b[18],&oy[n]);
  sscanf(&b[31],"%f %f",&o[n][0],&o[n][1]);              // x.y
  if(b[51]!=32)sscanf(&b[47],"%f %f %f",&o[n][2],&o[n][3],&o[n][4]);
  else o[n][2]=o[n][3]=o[n][4]=0.;
  sscanf(&b[68],"%f %f %f %f %f %f %f %f %d", &o[n][5],&o[n][6],&o[n][7],
                 &o[n][8],&o[n][9],&o[n][10],&o[n][11],&o[n][12],&om[n]);
  n++;
  goto l30;
l40:
  fclose(fp);
  sortd2(ox,oy,o,om,n);
// put bad_pixel mark
  k=strlen(f1); k=f1[k-5]-48; 
  if(k==1)putmark(bad1,n);
  if(k==2)putmark(bad2,n);
  if(k==3)putmark(bad3,n);
  if(k==4)putmark(bad4,n);
// sl
  m=0;
  fp=fopen("um.cat","r");
l10:
  fgets(a,80,fp); if(feof(fp))goto l20;
  if(a[0]=='#')goto l10;
  sscanf(&a[32],"%10d%9d %f %f %f %f",
    &si1[m],&si2[m],&sf1[m],&sf2[m],&sf3[m],&sf4[m]);
  chs(&a[6],&sx[m]); chs(&a[19],&sy[m]); m++;
  goto l10;
l20:
  fclose(fp);
// match
  ip=0;
l50:
  t0=t1=0.; k=0;
  for(j=l=0;j<n;j++){
    kk[j]=-1; sig=sig2; 
    for(i=0;i<m;i++){
      yy=sy[i]-oy[j]; yy*=3600.; if(yy>1.)break; 
      zy=yy; yy*=yy; if(yy>sig)continue;
      xx=sx[i]-ox[j]; if(xx> 12.)xx-=24.; if(xx<-12.)xx+=24.;
      xx=xx*dc; zx=xx; xx*=xx;
      zz=xx+yy; if(zz<sig){ sig=zz; kk[j]=i; dx=zx; dy=zy; }	    
    }
    if(kk[j]!=-1){ l++; t0+=dx; t1+=dy; }
    if(bp[j]==1)k++;
  }
  t0/=n; t1/=n;
  printf("%s matched: %d/%d/%d shift: %6.3f %6.3f\n",f1,l,k,n,t0,t1);
  if(ip++<1){
    for(i=0;i<n;i++){ ox[i]+=t0/dc; oy[i]+=t1/3600.; }
    goto l50;
  }
// output q*.cat
  for(i=0;i<n;i++){
    k=kk[i];
    if(k!=-1){ ox[i]=sx[k]; oy[i]=sy[k]; }
    toms2(ox[i],c1,1); toms2(oy[i],c2,0); c=32; if(bp[i]==0)c='*';
    for(j=0;j<180;j++)d[j]=32;
    sprintf(d,"%c%s %s %7.2f %7.2f",c,c1,c2,o[i][0],o[i][1]); d[43]=32;
    if(o[i][2]!=0.)sprintf(&d[43],"%8.2f%7.3f%6.3f",o[i][2],o[i][3],o[i][4]);
    sprintf(&d[64],"%6.2f%8.2f%7.3f%6.3f%5.2f%7.2f%5.2f%6.1f%3d",
    o[i][5],o[i][6],o[i][7],o[i][8],o[i][9],o[i][10],o[i][11],o[i][12],om[i]);
    if(k!=-1){  // add SL data
      sprintf(&d[117]," %10d%09d %6.3f%6.3f%7.3f%6.3f",
             si1[k],si2[k],sf1[k],sf2[k],sf3[k],sf4[k]);
    }
    fprintf(fp2,"%s\n",d);
  }
  fclose(fp2);
  goto l00;
}  

putmark(bad,n)
char bad[n2][n1]; int n;
{
  int i,j,k,ix,iy;
  for(i=0;i<n;i++){ k=1;
    ix=o[i][0]-0.5; iy=o[i][1]-0.5;  // should +0.5 round, bad1[][] is 1 diff coord. so -0.5
    for(j=ix-2;j<=ix+2;j++)if(bad[iy][j]==0)k=0;
    for(j=iy-2;j<=iy+2;j++)if(bad[j][ix]==0)k=0;
    if(bad[iy-1][ix-1]==0)k=0;
    if(bad[iy-1][ix+1]==0)k=0;
    if(bad[iy+1][ix+1]==0)k=0;
    if(bad[iy+1][ix-1]==0)k=0;
//    if(k==0)printf("%7.2f %7.2f\n",o[i][0],o[i][1]);
    bp[i]=k;
  }
}

chs(a,y)
char *a; double *y;
{
  double x;
  char b[128];
  int i,j,k;
  strcpy(b,a); k=strlen(b);
  j=1; for(i=0;i<k;i++){ if(b[i]==32)continue; if(b[i]=='-')j=-1; break; }
  for(i=0;i<k;i++)if(b[i]<'.' || b[i]> '9')b[i]=32;
  i=k=0; x=0.; sscanf(b,"%d %d %lf",&i,&k,&x);
  x=k/60.+x/3600.+i;
  *y=x*j;
}

toms2(aa,c,k)
double aa;           // input
char c[14];          // output
int k;               // input
{
  int  i,j;
  double a,b,x;
  a=aa;
  if(k==1){ if(a>24.)a-=24.; if(a<0.)a+=24.; }
  c[0]=32; if(a<0.){ a=-a; c[0]='-'; }
  i=a; b=(a-i)*60.;
  j=b; x=(b-j)*60.;
  if(j==60){ j=0; i++; }
  if(k==1)sprintf(&c[1],"%2.2d:%2.2d:%6.3f",i,j,x);
  if(k==0)sprintf(&c[1],"%2.2d:%2.2d:%5.2f",i,j,x);
  if(c[7]==32)c[7]='0';   if(c[8]==32)c[8]='0';
  if(c[7]=='6'){
    c[7]='0'; j++; sprintf(&c[4],"%2.2d",j);
    if(c[4]=='6'){
      c[4]='0'; i++; sprintf(&c[1],"%2.2d",i);
    }
    c[3]=':'; c[6]=':';
  }
}

sortd2(x,y,o,om,n)
double *x,*y; float o[][13]; short *om; int n;
{
  int i,j,k=2,ii,ifin;
  double xx,yy;
  float  oo[13];
  short  mm;
  int l;
  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(y[i]<=y[j])continue;
    xx = x[j];   yy= y[j];   mm=om[j];
    for(l=0;l<13;l++)oo[l]=o[j][l];
l40:
    x[j]= x[i];  y[j]= y[i]; om[j]=om[i];
    for(l=0;l<13;l++)o[j][l]=o[i][l];
    j=i;       i-=k;
    if(i<0)goto l60;
    if(y[i]>yy)goto l40;
l60:
    x[j]=xx;     y[j]=yy;    om[j]=mm;
    for(l=0;l<13;l++)o[j][l]=oo[l];
  }
  if(k>1)goto l20;
}

readbad(a,m)
char a[n2][n1]; int m;
{
  int i,j,k;
  char f1[60];
  sprintf(f1,"/line3/uband-data/cat/ubadpixel/ubad_%d.fit",m);
  fp=fopen(f1,"rb"); 
  if(fp==0){ printf("%s not found!\n\n",f1);
    sprintf(f1,"ubad_%d.fit",m);
    fp=fopen(f1,"rb");
    if(fp==0)exit(0);
  }
  fseek(fp,36*80,0); fread(a,n1,n2,fp); fclose(fp);
}