#include <stdio.h>
#include <math.h>
#define n 30000              // usl star number
#define bok 48.              // standard SCHMIDT system,type=1  
#define r05  1.0             // 
float  buf[13000*13000];
char   head[72][80],f1[30],b[30],c[80],d[150], mark[n],c0=34,c1=37,cc;
FILE   *fp,*fp1,*fp2;
double a8[8],b8[8],xx[n],yy[n],x,y,x3,x4,r,cx,cy,pi;  // position usl
float  usl[n],ber[n],sky,see,gain,x1,x2,ava,sig, z1,z2;   // mag
int    n1,n2, i0,i1,ic, m,pg=0,cg=0; // int ccd;
float  scale,r0,rr;
main(int ac,char **av)
{
  int i,j,k;
  if(ac<2){
    printf("\n\t******* use USL to calibrating & check matched***********\n");
    printf("\n\tUsage: c9 *.fit [!,#,%,$]\n");
    printf("\n\t ! plot fit, # gif fit, %c plot match, $ gif match\n\n",c1);
    exit(0);
  }
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  system("cp /vega2/rhbin/uefault.param default.param");
  system("cp /vega2/rhbin/uefault.conv default.conv");
  system("cp /vega2/rhbin/uefault.nnw default.nnw");
  cc=av[ac-1][0]; if(cc=='!' || cc=='#' || cc=='%' || cc=='$'){ ac--; pg=1; }
                
  for(i=1;i<ac;i++){
    strcpy(f1,av[i]);
    printf("%3d: %s ******\n",i,f1);
    fp=fopen(f1,"rb+");
    if(fp==0){ printf("not found!\n");  continue; }
    fread(head,72,80,fp);
    k=indexpos(head,"A81     ",72);
    for(j=0;j<8;j++)sscanf(&head[k+j][10],"%lf",&a8[j]);
    xytoad(a8,b8);
    sscanf(&head[3][20],"%d",&n1);
    sscanf(&head[4][20],"%d",&n2);
    k=indexpos(head,"SCALE   ",72);
    if(k==72){ printf("\n\t Scale not found!\n"); exit(0); }
    else sscanf(&head[k][20],"%f",&scale);
/*
    k=indexpos(head,"CCD_NO: ",72); sscanf(&head[k][15],"%d",&ccd);
    k=indexpos(head,"SEEING  ",72); sscanf(&head[k][15],"%f",&see);
    k=indexpos(head,"GAIN    ",72); 
    sscanf(&head[k][10],"%f %f %lf %lf",&x1,&x2,&x3,&x4); gain=(x1+x2+x3+x4)/4.
    
    strncpy(b,&head[6][53],25);
    for(j=0;j<9;j++)if(b[j]==32)break;
    if(j!=9){ printf("not coordinated!\n"); continue; }
*/
    see=2.0; gain=1.8;
    for(j=0;j<13;j++){ b[j]=head[11][17+j]; b[j+13]=head[12][17+j]; } b[26]=0;
    k=n1; if(n1<n2)k=n2;
    z1=k*scale; z1/=60.; z1+=0.5; k=z1;
    r0=r05/scale;
    sprintf(c,"usl %s %d c9.usl 22.0",b,k);
    system(c);                    
    m=0;
    fp1=fopen("c9.usl","r");
l10:
    fgets(d,90,fp1); if(feof(fp1))goto l20;
    if(d[0]=='#')goto l10;
    sscanf(&d[64],"%f",&usl[m]);  mark[m]=0;
    chs(&d[6],&x3);  chs(&d[19],&x4);
    rad_xy(a8[6],a8[7],x3,x4,&x,&y,1,0);
    if(x<10. || x>n1-10.)goto l10;
    if(y<10. || y>n2-10.)goto l10;
//    if(ccd==4)if(x<n1/2-10. || y>n2/2-10.)goto l10;
    xx[m]=x; yy[m]=y;
    if(m<n)m++;
    goto l10;
l20:
    fclose(fp1); 
    if(m<2){ printf("less U standard star in this area\n"); continue; }
  
    fp2=fopen("/vega2/rhbin/uefault.sex","r");
    fp1=fopen("default.sex","w");
l30:
    fgets(c,80,fp2); if(feof(fp2))goto l40;
    if(c[0]=='G'&&c[1]=='A'&&c[2]=='I'){sprintf(&c[6],"%4.1f",gain);c[10]=9;}
    if(c[0]=='S'&&c[1]=='E'&&c[2]=='E'){sprintf(&c[12],"%5.2f",see);c[17]=9;}
    fputs(c,fp1); goto l30;
l40:
    fclose(fp1);
    sprintf(c,"sex %s",f1);
    system(c);                           // get test.cat
    fp1=fopen("test.cat","r");
    rr=0.;
l50:
    fgets(d,150,fp1);  if(feof(fp1))goto l60;
    if(d[0]=='#')goto l50;
    sscanf(&d[48],"%f %f %lf %lf",&x1,&x2,&x3,&x4);
    if(x1>25.)goto l50;
    r=r0;                      // err 1 pixel, 0.454"
    for(j=0;j<m;j++){
      x=xx[j]-x3; x*=x; 
      y=yy[j]-x4; y*=y; 
      x=sqrt(x+y);
      if(x<r){r=x; k=j; }
    } 
    if(r!=r0){
      ber[k]=usl[k]-x1; mark[k]=2;
      if(usl[k]>21. || usl[k]<16.)mark[k]=1;
//      printf("%5.2f",r*scale);
      rr+=r*scale;
    }
    goto l50;
l60:
    fclose(fp1);
    for(j=k=0;j<m;j++)if(mark[j]) k++; 
    printf("\nusl_star: %d matched star: %d  rms: %5.2f%c",m,k,rr/k,c0);
    printf("   (matched aperture 1.00%c)\n",c0);
    printf("by using usl mag: 16.~21m\n\tmagnitude match: ");
l70:
    x1=x2=0.;
    for(j=k=0;j<m;j++)if(mark[j]==2){ k++; x1+=ber[j]; }
    printf(" %4d ",k);
    ava=x1/k;  
    for(j=0;j<m;j++)if(mark[j]==2) x2+=(ber[j]-ava)*(ber[j]-ava);
    sig=sqrt(x2/k); 
    printf(" av: %6.3f %6.3f\n",ava,sig);
    sig*=3.;
    for(j=ic=0;j<m;j++)if(mark[j]==2){
      x2=ber[j]-ava; if(x2<0.)x2=-x2;
      if(x2>sig){ mark[j]=1; ic++; }
    }
    if(ic){ printf("\t  3_sigma_match: "); goto l70;  }
    sig/=3.;  r=pow(10,(ava-4.)*0.4);
    sprintf(&head[29][0],"CALIBRAT=%21.3f / rms:%5.3f  usl:%4d  matched:%4d  W:%7.4lf",
     ava,sig,m,k,r);    for(j=77;j<80;j++)head[29][j]=32;
    printf("CALIBRAT=%8.3f  Weight: %7.4lf\n",ava,r);
    if(k<6){ for(j=0;j<80;j++)head[29][j]=32;
              printf("******* less matched (<6), not write down to file\n");
              continue;
           }
// mark==2, 16-21 mastched star, use to cal seeing
    fread(buf,n1*n2,4,fp); swap4(buf,n1*n2*4);
    j=indexpos(head,"SKYADU  ",72); sscanf(&head[j][20],"%f",&sky);
    seeing(buf,n1,n2,xx,yy,mark,m,sky,&see,scale);
    printf("SKY: %7.2f  Seeing:  %5.2f\n",sky,see);
    sprintf(&head[j+1][0],"SEEING  =%21.3f /",see); head[j+1][32]=32;

    fseek(fp,0,0); fwrite(head,72,80,fp);
    fclose(fp);
    if(pg){
      i0=0; i1=1;
      if(cc=='!'|| cc=='%')pgbegin_(&i0,"/xw",&i1,&i1,3);
      if(cc=='#'|| cc=='$')pgbegin_(&i0,"/jg",&i1,&i1,3);
      if(cc=='!' || cc=='#'){
        x1=800., x2=600.;        pgpap_(&x1,&x2);
        x1=14.; x2=23.; z1=x1-ava; z2=x2-ava;
        pgsci_(&i1); pgenv_(&x1,&x2,&z1,&z2,&i0,&i0);
        sprintf(c,"file: %s   green points: %d",f1,k);
        pglabel_("USL mag","Bertin mag",c,7,10,strlen(c));
        sprintf(c,"calibration:%9.3f    rms:%9.3f",ava,sig);
        x1=14.5; z1=z2-1.; pgtext_(&x1,&z1,c,strlen(c));
      sprintf(c,"Weight(FLUX):%8.3f",r);
      z1=z2-1.7; pgtext_(&x1,&z1,c,strlen(c));
      sprintf(c,"SkyAdu:%8.2f",sky);
      z1=z2-2.4; pgtext_(&x1,&z1,c,strlen(c));
      sprintf(c,"Seeing:%8.2f",see);
      z1=z2-3.1; pgtext_(&x1,&z1,c,strlen(c));
        for(j=0;j<m;j++)if(mark[j]){
          ic=mark[j]*2-1; pgsci_(&ic);
          x1=usl[j]; x2=usl[j]-ber[j];
          pgpoint_(&i1,&x1,&x2,&i1);
        }
      }
      if(cc=='%' || cc=='$'){
        x1=600.; x2=600.;        pgpap_(&x1,&x2);
        x1=z1=1.; x2=n1; z2=n2;
        pgsci_(&i1); pgenv_(&x1,&x2,&z1,&z2,&i0,&i0);
        sprintf(d,"%s   0.5%c matched objects %d",f1,c0,k);
        pglabel_("PIXEL(R.A.)","PIXEL(Dec.)",d,11,11,strlen(d));
        for(j=0;j<m;j++)if(mark[j]){
//          ic=mark[j]*2-1; pgsci_(&ic);
          z1=xx[j]; z2=yy[j]; pgpoint_(&i1,&z1,&z2,&i1);
        }
      }
      pgend_();
    }    
  }
}

seeing(map,n1,n2,xx,yy,mark,m,sky,see,scale)
char  *mark;
int   n1,n2,m;
float map[n2][n1],sky,*see,scale;
double *xx,*yy;
{
  short d[2000];
  float z,s1,s2;
  int i,j,k,ix,iy,b[113],l,ip;

  for(l=ip=0;ip<m;ip++)if(mark[ip]==2){
    ix=xx[ip]+0.5; iy=yy[ip]+0.5;
    for(j=0,i=ix-14;i<=ix+14;i++,j+=4){
    z=0.; for(k=-2;k<=2;k++)z+=map[iy+k][i]-sky; 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;
    for(j=0;j<112;j++)if(b[j]<0)b[j]=0;
    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+=map[i][ix+k]-sky; 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;
    for(j=0;j<112;j++)if(b[j]<0)b[j]=0;
    histat(b,&z,&s2,112);
    k=(s1+s2)/8.*2.355*100.+0.5;          // expand 100
    if(k >100 && k<999 )d[l++]=k;
    if(l>=2000)break;
  }
  white_black(d,1,l,&z,&s1); z*=0.01;        // shrink 100
  *see=z*scale;
}

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

xytoad(x,a)
double *x,*a;
{
  double z;
  z=x[0]*x[3]-x[2]*x[1];
  a[0]= x[3]/z;
  a[1]=-x[1]/z;
  a[2]=-x[2]/z;
  a[3]= x[0]/z;
  a[4]= (x[2]*x[5]-x[4]*x[3])/z;
  a[5]= (x[4]*x[1]-x[0]*x[5])/z;
}

rad_xy(ac,dc,ra,de,x,y,type,k)
double ac,dc,ra,de,*x,*y;  int type,k;
{
  double xi,xn, tmp,rar,der;
  rar=ra*cx; der=de*cy;
/*                         // same discribe formula as following 3 lines
  double sd,cd,td,co;      // china bai_ke_quan_shu (astronmy) p.552
  sd=sin(dc);  cd=cos(dc);  td=tan(der);
  co=cos(rar-ac);  tmp=sd*td+cd*co;
  xi=sin(rar-ac)/tmp;
  xn=(cd*td-sd*co)/tmp;
*/
  tmp=atan(tan(der)/cos(rar-ac));
  xi=cos(tmp)*tan(rar-ac)/cos(tmp-dc);
  xn=tan(tmp-dc);
  if(type==1){     // schmidt TELESCOPE
    tmp=sqrt(xi*xi+xn*xn);
    if(tmp!=0.){
      tmp=atan(tmp)/tmp;
      xi*=tmp; xn*=tmp;
    }
  }
  if(type==2){    //  BOK telescope
    tmp=1.+bok*(xi*xi+xn*xn);
    xi*=tmp; xn*=tmp;
  }
  if(k==0){
    *x=b8[0]*xi+b8[2]*xn+b8[4];
    *y=b8[1]*xi+b8[3]*xn+b8[5];
  }
  if(k==1){ *x=xi; *y=xn;}
}

xy_rad(ac,dc,x,y,ra,de,type,k)
double ac,dc,x,y,*ra,*de;  int type,k;
{
  double xi,xn,tmp,tmp1;
  double x1,y1,r1,d1;
  int i;
  xi=a8[0]*x+a8[2]*y+a8[4];
  xn=a8[1]*x+a8[3]*y+a8[5];
  if(type==1){
    tmp=sqrt(xi*xi+xn*xn);
    if(tmp!=0.){
      tmp=tan(tmp)/tmp;
      xi*=tmp;  xn*=tmp;
    }
  }
  if(type==2){
    tmp=1.+bok*(xi*xi+xn*xn);
    for(i=0;i<4;i++){
      x1=xi/tmp; y1=xn/tmp;
      tmp=tan(dc);  tmp1=1.-y1*tmp;
      r1=atan(x1/cos(dc)/tmp1);
      d1=atan((y1+tmp)*cos(r1)/tmp1);
      r1+=ac;
      tmp=atan(tan(d1)/cos(r1-ac));
      x1=cos(tmp)*tan(r1-ac)/cos(tmp-dc);
      y1=tan(tmp-dc);
      tmp=1.+bok*(x1*x1+y1*y1);
//                                 printf("%lf\n",tmp);
    }
    xi/=tmp; xn/=tmp;
  }
  tmp=tan(dc);  tmp1=1.-xn*tmp;
  *ra=atan(xi/cos(dc)/tmp1);
  *de=atan((xn+tmp)*cos(*ra)/tmp1);
  *ra+=ac;
  tmp=(*de)-dc; if(tmp<0.)tmp=-tmp;
  if(tmp>2.){  *de=-(*de);  *ra+=pi; }
  if(*ra<0.)*ra+=(pi+pi);
  if(k==0){ *ra/=cx; *de/=cy; }
}