#include <stdio.h>
#include <math.h>
#define n 100              // usl star number
#define bok 48.            // standard SCHMIDT system,type=1  
char   head[72][80],c[80],d[270],e[270];
FILE   *fp,*fp1;
double a8[8],b8[8],xx[n],yy[n],x,y,x3,x4,r,cx,cy,pi;  // position usl
int    n1,n2,i0,i1,i2,i3,i4,m,i20;
float  u1[n],u2[n],u3[n],u4[n];
float  v1[n],v2[n],v3[n],v4[n];
float  x1,x2,z1[4],z2[4];
char   c1[20],c2[20],c34=34;
float  s1,s2, see1,see2,flu1,flu2;
float  ss[8]={ 0., 2.26,1.92,1.89,1.92,2.71,2.61,2.13 };
float  ff[8]={0.,0.,1.514,1.542,1.201,1.629,1.615,1.179};
main(int ac,char **av)
{
  int i,j,k;
  if(ac<3){
    printf("\n\t****** test z_9,c_9 by using a6.fit *******\n");
    printf("\n\t\t\t\t\t2011,9,15");
    printf("\n\tUsage: testz_9 1 2  (1-8)\n");
    printf("\n\t1: usl");
    printf("\n\t2: a6.fit    (by using z_9)");
    printf("\n\t3: p5506_0053_1.fit");
    printf("\n\t4: p5506_0054_1.fit");
    printf("\n\t5: p5506_0056_4.fit");
    printf("\n\t6: p5508_0027_1.fit");
    printf("\n\t7: p5508_0028_1.fit");
    printf("\n\t8: p5508_0057_4.fit\n\n");
    exit(0);
  }
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  fp=fopen("a6.fit","rb"); if(fp==0){ printf("a6.fit not found!\n"); exit(0); }
  fread(head,72,80,fp); fclose(fp);
  sscanf(av[1],"%d",&i1);
  if(i1<1 || i1>8){ printf("par_1 out of range\n"); exit(0); } i1--;  
  sscanf(av[2],"%d",&i2);
  if(i2<1 || i2>8){ printf("par_1 out of range\n"); exit(0); } i2--; 
 
  fp=fopen("a6.all","r");
  fgets(e,270,fp);
  i=0;
l10:
  fgets(d,270,fp); if(feof(fp))goto l20;
  if(d[0]=='#')goto l10;
  sscanf(&d[i1*33+4],"%f %f %f %f",&u1[i],&u2[i],&u3[i],&u4[i]);
  sscanf(&d[i2*33+4],"%f %f %f %f",&v1[i],&v2[i],&v3[i],&v4[i]);
//  printf("%f %f %f -  %f %f %f\n",u1[i],u2[i],u3[i],v1[i],v2[i],v3[i]);  
  i++;
  goto l10;
l20:
  m=i; fclose(fp); 
  for(i=0;i<19;i++)c1[i]=e[i1*33+7+i]; c1[i]=0;
  for(i=0;i<19;i++)c2[i]=e[i2*33+7+i]; c2[i]=0;
//  printf("%s\n%s\n",c1,c2);
  see1=ss[i1]; see2=ss[i2];
  flu1=ff[i1]; flu2=ff[i2];

  k=indexpos(head,"A81     ",72);
  for(j=0;j<8;j++)sscanf(&head[k+j][10],"%lf",&a8[j]);
  xytoad(a8,b8);
  for(i=0;i<m;i++){
    x=u1[i]; y=u2[i];
    rad_xy(a8[6],a8[7],x,y,&x3,&x4,1,0);
    u1[i]=x3; u2[i]=x4;
    x=v1[i]; y=v2[i];
    rad_xy(a8[6],a8[7],x,y,&x3,&x4,1,0);
    v1[i]=x3; v2[i]=x4;
//    printf("%f %f - %f %f\n",u1[i],u2[i],v1[i],v2[i]);  
  }
  i0=0; i1=1; i2=2; i3=3; i4=4; i20=24;
  pgbegin_(&i0,"/xw",&i2,&i1,3);
  x1=960., x2=500.;        pgpap_(&x1,&x2);
  x1=14.; x2=23.;
  pgsci_(&i1); pgenv_(&x1,&x2,&x1,&x2,&i0,&i0);
  z1[0]=z2[0]=x1; z1[1]=z2[1]=x2;
  pgline_(&i2,z1,z2,&i1);
  s1=s2=0; pgsci_(&i3);
  for(j=k=0;j<m;j++){
    x1=u3[j]; x2=v3[j];
    pgpoint_(&i1,&x1,&x2,&i20);
    if(u3[j]>5. && v3[j]>5.){ k++;
      s1+=u4[j]; s2+=v4[j];
    }
  }
  s1/=k; s2/=k;
  sprintf(c,"average MAG_RMS: %5.3f",s1); 
  x1=18.; x2=16.; pgtext_(&x1,&x2,c,strlen(c));
  x2-=0.6; sprintf(c,"seeing: %6.2f",see1);
  if(see1>0.)pgtext_(&x1,&x2,c,strlen(c));
  x2-=0.6; sprintf(c,"fulx_r: %6.3f",flu1);
  if(flu1>0.)pgtext_(&x1,&x2,c,strlen(c));
  sprintf(c,"average MAG_RMS: %5.3f",s2); 
  x1=15.; x2=22.; pgtext_(&x1,&x2,c,strlen(c));
  x2-=0.6; sprintf(c,"seeing: %6.2f",see2);
  if(see2>0.)pgtext_(&x1,&x2,c,strlen(c));
  x2-=0.6; sprintf(c,"fulx_r: %6.3f",flu2);
  if(flu2>0.)pgtext_(&x1,&x2,c,strlen(c));
  sprintf(c,"all matched star: %d",k);
  pgsci_(&i1); pglabel_(c1,c2,c,19,19,strlen(c));

  x1=0.; x2=1536.;
  pgsci_(&i1); pgenv_(&x1,&x2,&x1,&x2,&i0,&i0);
  r=0.;
  for(j=k=0;j<m;j++){
    x1=u1[j]; x2=u2[j];
    pgsci_(&i3); pgpoint_(&i1,&x1,&x2,&i20);
    z1[0]=x1; z2[0]=x2;
    x3=x1-v1[j]; x4=x2-v2[j];
    z1[1]=x3*45.4+x1;
    z2[1]=x4*45.4+x2;
    if(v3[j]>5. && u3[j]>5.){
      pgsci_(&i2); pgline_(&i2,z1,z2,&i1);
      r=r+sqrt(x3*x3+x4*x4); k++;
    }
  }
  z1[0]=z1[1]=1300; z1[2]=z1[3]=1400;
  z2[0]=z2[3]=80;  z2[1]=z2[2]=50;
  pgsci_(&i3); pgline_(&i4,z1,z2,&i1);
  x1=1320; x2=80; sprintf(c,"1%c",c34); pgtext_(&x1,&x2,c,2);
  r=r/k*0.454;
  sprintf(c,"position RMS: %6.2lf%c",r,c34);
  pgsci_(&i1); pglabel_("","",c,0,0,strlen(c));
  pgend_();
}

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