#include <stdio.h>
FILE *fp1,*fp2,*fp3;
char a[160],c1[14],c2[14];
char head[72][80];
double a8[8],b8[8];
double rc,dc,r1,d1,x,y, pi,cx,cy;

main(ac,av)
int ac; char *av[];
{
  int i,j,k;
  float skyadu,gain,rdno;

  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  if(ac<4){printf("\n\t **** ta02u(uband_shrink3)+TA02i.mag2 uedo.coo ****\n");
         printf("\n\tUsage: chmag22coo ta02u.fit Ta02i_15.mag2 ta02u.coo\n\n"); 
            exit(0); }

  fp1=fopen(av[1],"rb"); if(fp1==0){printf("%s not found!\n",av[1]); exit(0);}
  fp2=fopen(av[2],"r");  if(fp2==0){printf("%s not found!\n",av[2]); exit(0);}
  fp3=fopen(av[3],"w");  
  fread(head,72,80,fp1); fclose(fp1);
  k=indexpos(head,"A81     ",72);
  if(k==72){ printf("%s not coordination !\n",av[1]); exit(0); }
  for(i=0;i<8;i++)sscanf(&head[k+i][10],"%le",&a8[i]);
  xytoad(a8,b8);  rc=a8[6]; dc=a8[7]; 
  k=indexpos(head,"SKYADU  ",72);
  skyadu=99; if(k!=72)sscanf(&head[k][10],"%f",&skyadu);
  gain=2.; rdno=8.; 
  fprintf(fp3," NL   NX   NY  LOWBAD HIGHBAD  THRESH");
  fprintf(fp3,"     AP1  PH/ADU  RNOISE    FRAD\n");
  fprintf(fp3,"  1 4096 4096  %6.1f 60000.0   99.99    0.00",skyadu*0.8);
  fprintf(fp3,"%8.2f%8.2f    3.00\n\n",gain,rdno);
  j=0;
l10:
  fgets(a,160,fp2);
  fgets(a,160,fp2); if(feof(fp2))goto l20;
  sscanf(&a[7],"%s %s",c1,c2);
  chs(c1,&r1); chs(c2,&d1);
  rad_xy(rc,dc,r1,d1,&x,&y,1,0);
//    printf("%s %s  %8.2lf %8.2lf\n",c1,c2,x,y);
//  fprintf(fp3,"%6d %8.3lf %8.3lf %s %s\n",++j,x,y,c1,c2);
  fprintf(fp3,"%6d %8.3lf %8.3lf\n",++j,x,y);
  goto l10; 
l20:
  fclose(fp3);
  printf("%s produced!\n",av[3]);
}

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

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

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