#include <stdio.h>
#include <math.h>
#define n  9999
#define mm 40000               // for star work
#define bok  48.

FILE   *fp,*fp1,*fp2,*fp3;
char   a[80],f1[60],b[80],h[36][80],c;
double xx[n],yy[n],zz[n];
char   ww[n];
double cx,cy,pi;
double r1[mm],d1[mm],r2[mm],d2[mm];
double m1[mm],m2[mm],ava[99],sig[99];
int    im1,im2,pd[99],ik[99];
double a8[8],b8[8];
char   head[72][80],w[mm];
float  gain,see;

outcat(l,ra,de,mag,nn)
int l,*nn; 
double *ra,*de,*mag;
{
  int i,j,k;
  char a[46],c[80],d[150];
  double x1,x2,x3,x4;
  fseek(fp,l*53,0); fread(a,1,45,fp); a[45]=0;
  printf("%s\n",a);
  fp1=fopen(a,"rb"); fread(head,72,80,fp1); fclose(fp1);
  k=indexpos(head,"A81     ",72);
  for(j=0;j<8;j++)sscanf(&head[k+j][10],"%lf",&a8[j]);
  xytoad(a8,b8);
  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.;

  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); fclose(fp2);
  sprintf(c,"sex %s",a); system(c);
  fp1=fopen("test.cat","r");
  j=0;
l50:
  fgets(d,150,fp1);  if(feof(fp1))goto l60;
  if(d[0]=='#')goto l50;
  sscanf(&d[48],"%lf %lf %lf %lf",&x1,&x2,&x3,&x4);
  if(x1>25.)goto l50;
  mag[j]=x1;
  xy_rad(a8[6],a8[7],x3,x4,&x1,&x2);
  if(x1<12.)x1+=24.; ra[j]=x1; de[j]=x2;
//  printf("%d %lf %lf %lf\n",j,mag[j],ra[j],de[j]);
  j++;
  goto l50;
l60:
  fclose(fp1);
  *nn=j;
}

work(ava,sig,m)
double *ava,*sig;
int *m;
{
  int i,j,k,kk;
  double x,y,x1,y1,z1,r,r0;
// ?1[] not change, ?2[] may change
  r0=0.5;
  kk=0;
  for(i=0;i<im2;i++){ w[i]=0;
    x=r2[i]; y=d2[i]; r=r0;
    for(j=0;j<im1;j++){
      x1=r1[j]-x;  x1*=cos(y*cy)*54000.;
      y1=d1[j]-y;  y1*=3600.;
      z1=sqrt(x1*x1+y1*y1);
      if(z1<r){ k=j; r=z1; }
    }
    if(r!=r0){ kk++; w[i]=1; m2[i]-=m1[k];}               // ccd_stand - CCDX
  }
  printf("position match: %d\n",kk);
l10:
  x=y=0.; k=0;
  for(i=0;i<im2;i++)if(w[i]){ k++; x+=m2[i]; } 
  if(k==0){ *m=0; *ava=99.; *sig=0; return; }
  x/=k;  
  for(i=0;i<im2;i++)if(w[i])y+=(m2[i]-x)*(m2[i]-x);  
  y/=k;
  printf("3sig match: %d %lf %lf\n",k,x,y);
  r0=3*y; j=0;
  for(i=0;i<im2;i++)if(w[i]){
    r=m2[i]-x; if(r<0.)r=-r;
    if(r>r0){ w[i]=0; j++; }
  }
  if(j && y>0.02)goto l10;
  *ava=x; *sig=y; *m=k;
}

main(int ac,char **av)
{
  int i,j,k,k1,k2,k3,ip=0,id;
  double x,y,x1,y1,z1;
  if(ac<2){ printf("\n\t star cal1 (2011)****\n");
            printf("\n\t Usage: starcal fit3.cal2,4,6,8,...\n\n");
            exit(0);
  }
		fp3=fopen("star.cal","w");
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  sprintf(b,"%s",av[1]);
  fp1=fopen(b,"r");
  if(fp1==0){ printf("%s not found !\n",b); exit(0); }
l10:
  fgets(a,80,fp1); if(feof(fp1))goto l20;
  x=0.; sscanf(a,"%s %lf",f1,&x);
  k=strlen(f1); c=f1[k-5];
  fp=fopen(f1,"rb");
  fread(h,36,80,fp); fclose(fp); 
  h[6][65]=h[6][78]=0;
  chs(&h[6][53],&x1);
  chs(&h[6][66],&y1);
  if(x1<12.)x1+=24.;
  xx[ip]=x1; yy[ip]=y1; zz[ip]=x; ww[ip]=c;
  ip++;
  goto l10;
l20:
  fclose(fp1);
  for(i=j=0;i<n;i++)if(zz[i]==0.)j++;
  printf("all ccd: %d  no_cal1: %d\n",ip,j);

  system("cp /vega2/rhbin/uefault.param default.param");
  system("cp /vega2/rhbin/uefault.conv default.conv");
  system("cp /vega2/rhbin/uefault.nnw default.nnw");
  fp=fopen(b,"rb");                // fit1.cal2...
  k1=k2=0;
  for(i=0;i<n;i++)if(zz[i]==0.){
    x=xx[i]; y=yy[i];
    id=0;
    for(j=0;j<n;j++)if(zz[j]!=0.){
      x1=xx[j]-x;  x1*=cos(y*cy)*900.;
      y1=yy[j]-y;  y1*=60.;
        
      z1=sqrt(x1*x1+y1*y1);
//if(z1<60.)printf("%f ",z1);
      if(z1>30.)continue;
      pd[id++]=j;
    }
    if(id==0){ k2++; x1=0.; goto l99; }
    printf("%d_%c:\n",i,ww[i]);
    outcat(i,r1,d1,m1,&im1);           // get mother's ra[].de[],mag[],
         printf("MMMMM***** %d\n",im1);
    for(j=0;j<id;j++){
      k=pd[j];
      outcat(k,r2,d2,m2,&im2);         // get each ra[n2],...
         printf("CCCCC***** %d\n",im2);
      work(&ava[j],&sig[j],&k3);       // ik[],  matched num      
      ava[j]+=zz[k];  ik[j]=k3;        // change to write down value.
      printf("%d %lf %lf\n",k3,ava[j],sig[j]);
      fprintf(fp3,"%d %lf %lf\n",k3,ava[j],sig[j]);
    }
// single ccd_match,  if k3 <20 stars, all k3<50, ignore
// sigma not count in,  we count k3, as weight
    k3=0; x1=0.;
    for(j=0;j<id;j++)if(ik[j]>20){ k3+=ik[j]; x1+=ava[j]*ik[j];}
    if(k3<50){ k2++; x1=0.; goto l99; }
    x1/=k3;
    k1++;
// x1= last ava
    fseek(fp,i*53,0); fread(a,1,45,fp); a[45]=0;
    printf("I will write: %s %lf\n",a,x1);
    fprintf(fp3,"%s %lf\n",a,x1);
    fp1=fopen(a,"rb+");
    if(fp1!=0){
      fread(head,72,80,fp1);
      for(j=0;j<80;j++)head[29][j]=32;
      sprintf(&head[29][0],"CALIBRAT=%21.3lf / star_cal1",x1);
      x1=pow(10,(x1-4.)*0.4);
      sprintf(&head[29][68],"W:%7.4lf",x1);
      for(j=32;j<80;j++)if(head[29][j]==0)head[29][j]=32;
      fseek(fp1,0,0); fwrite(head,72,80,fp1);  fclose(fp1);
    }
l99:  ;
  }
  printf("done: %d   none: %d\n",k1,k2);
  fclose(fp3);
}

chs(a,y)
char *a; double *y;
{
  double x;
  char b[20];
  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;
}

xy_rad(ac,dc,x,y,ra,de)
double ac,dc,x,y,*ra,*de;
{
  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];
  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);
  }
  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);
  *ra/=cx; *de/=cy;
}