#include <stdio.h>
#include <math.h>
double a8[8],b8[8],pi,cx,cy;
FILE   *fp,*fp1,*fp2;
char   head[72][80],d[150],f1[50],f2[50],c1[14],c2[14];
int    n1,n2,m=19;
int    ft[19]={10,11,12,13,16,25,32,33,34,35,36,37,38,39,40,41,42,43,44,45};

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

#define bok  48.
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; }
}

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]=':';
  }
}

main(int ac, char **av)
{
  int    i,j,k,l;
  char   c[80];
  float  see,gain;
  double x3,x4,rc,dc,alpha,delta;
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  if(ac<2){ printf("\n\t** UBAND shrink3 (4k) ---> p*.cat (bertiine)**\n");
            printf("\n\tUsage: u4k2cat ta02u.fit\n\n"); exit(0); }
  system("cp /vega2/rhbin/uefault.param default.param");
  system("cp /vega2/rhbin/uefault.conv default.conv");
  system("cp /vega2/rhbin/uefault.nnw default.nnw");

  for(i=1;i<ac;i++){
    strcpy(f1,av[i]); printf("%4d: %s ",i,f1);
    k=strlen(f1); if(f1[k-3]!='f' || f1[k-2]!='i' || f1[k-1]!='t')goto l99;
    fp =fopen(f1,"rb"); fread(head,72,80,fp); fclose(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);
    if(n1!=n2 && n1!=4096){ printf(" not a 4k file !\n"); continue; }
    rc=a8[6]; dc=a8[7];
    see=2.0; gain=1.8;
    k=indexpos(head,"SEEING  ",72);
    if(k!=72) sscanf(&head[k][15],"%f",&see);
    k=indexpos(head,"GAIN    ",72); 
    if(k!=72)sscanf(&head[k][15],"%f",&gain);

    fp1=fopen("/vega2/rhbin/default.blc","r");
    fp2=fopen("default.sex","w");
l30:
    fgets(c,80,fp1); if(feof(fp1))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,fp2); goto l30;
l40:
    fclose(fp1); fclose(fp2);
    sprintf(c,"sex %s",f1); system(c);                          

    strcpy(f2,f1);  k=strlen(f2);  f2[k-3]='c';  f2[k-2]='a'; 
    fp2=fopen(f2,"w");
    fprintf(fp2,"# %s           / SEX using /vega2/rhbin/default.blc\n",f1);
    for(j=0;j<m;j++){ k=ft[j];
      for(l=0;l<79;l++)c[l]=head[k][l]; c[l]=0;
      fprintf(fp2,"# %s\n",c);
    }
    fprintf(fp2,"#-------------------------------------------------------\n");
    fp1=fopen("test.cat","r");
    for(j=0;j<15;j++){ fgets(d,150,fp1); fputs(d,fp2); }
    fprintf(fp2,"#  16 ra  (2000.0)\n");       
    fprintf(fp2,"#  17 dec (2000.0)\n");       
    fprintf(fp2,"#-------------------------------------------------------\n");
l50:
    fgets(d,150,fp1);  if(feof(fp1))goto l60;
    sscanf(&d[67],"%lf %lf",&x3,&x4); d[134]=0; 
    xy_rad(rc,dc,x3,x4,&alpha,&delta,1,0);
    toms2(alpha,c1,1); toms2(delta,c2,0);
    fprintf(fp2,"   %s %s %s\n",d,c1,c2);
    goto l50;
l60:
    fclose(fp1); fclose(fp2);
    printf("ok!");
l99: 
    printf("\n");
  }
}