#include <stdio.h>
#include <math.h>
FILE *fp2,*fp3,*fp4;
char a[130];
char c1[14],c2[14],head[72][80];
double a8[8],b8[8],pi,cx,cy;
#define bok 48.

main(ac,av)
int ac; char *av[];
{
  int   i,j,k;
  double sharp,x,y,dx,dy,bx,by,r,s,alpha,delta;
  char  f1[60],f2[60],f3[60],f4[60];
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
 
  if(ac<2){ printf("\n\t ********lastcat ******\n");
            printf("\n\tUsage: lastcat uedo [1.0]\n\n");
            exit(0);
          }
  strcpy(f1,av[1]);
  k=nindex(f1,"."); if(k>0)f1[k]=0;
  strcpy(f2,f1);  strcpy(f3,f1); strcpy(f4,f1); 
  strcat(f2,".fit"); strcat(f3,".ext"); strcat(f4,".cat");
  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); }
  sharp=1.0; if(ac>2)sscanf(av[2],"%lf",&sharp);
  fread(head,72,80,fp2); fclose(fp2);
  k=indexpos(head,"A81     ",72);
  if(k==72){ printf("%s not coordinating!\n",f2); exit(0); }
  for(i=0;i<8;i++)sscanf(&head[k+i][10],"%lf",&a8[i]);
  fp4=fopen("sort.par","r"); for(i=0;i<3;i++)fgets(a,130,fp4); fclose(fp4);
  sscanf(a,"%d",&j);
  fp4=fopen(f4,"w");
  puthead();
  k=0;
l10:
  fgets(a,130,fp3); if(feof(fp3))goto l20;
  sscanf(a,"%lf %lf %lf %lf %lf %lf %lf %lf",&bx,&by,&dy,&dx,&r,&r,&r,&s);
  x=dx; y=dy; if(s>sharp){ x=bx; y=by; }
  xy_rad(a8[6],a8[7],y,x,&alpha,&delta);
  toms2(alpha,c1,1); toms2(delta,c2,0);
  if(s>sharp){ sprintf(&a[16],"%8.2lf%8.2lf",y,x); for(i=32;i<53;i++)a[i]=32; }
  k++;
  sscanf(&a[107],"%d",&i); a[106]=0; if(j==0)i=k;
  fprintf(fp4,"%4d%s %s%s\n",i,c1,c2,&a[16]);
  goto l10; 
l20:
  fclose(fp3);
  printf("%d  %s produced!\n",k,f4);
}

puthead()
{
  FILE *fp;
  int  i,j;
  fp=fopen("headcat.par","r");
  if(fp==0){ printf(" Cannot found headcat.par!\n"); exit(0); }
l10:
  fgets(a,130,fp); if(feof(fp))goto l20;
  if(a[0]<=32)goto l10;
  if(a[0]=='#'){ fputs(a,fp4); goto l10; }
  for(i=0;i<72;i++){
    for(j=0;j<strlen(a)-1;j++)if(head[i][j]!=a[j])goto l15;
    fputc('#',fp4); fwrite(&head[i][0],1,79,fp4); fprintf(fp4,"\n");
    goto l10;
l15:;
  }
  goto l10;  
l20:
  fclose(fp);
}

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

toms2(aa,c,k)
double aa;           // input
char c[14];          // output
int k;               // input
/*
 hour or degree to char_line
 if k=1 (hour case) in **:**:**.***
    k=0 (degree       -**:**:**.**
*/
{
  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]=':';
  }
}

rad_xy(ac,dc,ra,de,x,y)
double ac,dc,ra,de,*x,*y;
{
  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);
  tmp=1.+bok*(xi*xi+xn*xn);
  xi*=tmp; xn*=tmp;
  *x=b8[0]*xi+b8[2]*xn+b8[4];
  *y=b8[1]*xi+b8[3]*xn+b8[5];
}

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