#include <stdio.h>
#include <math.h>
char c1[14],c2[14],head[72][80];
double a8[8],pi,cx,cy;
FILE *fp,*fp1,*fp2;
int  n1,n2;
#define bok 48.

main()
{
  char a[50],f1[50];
  double x,y,alpha,delta;
  double a1,a2,a3,a4,d1,d2,d3,d4;
  int i,k;
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  system("ls p*.fit >v_1.tmp");
  fp1=fopen("v_1.tmp","r");
  fp2=fopen("v_1.cat","w");
l10:
  fgets(a,50,fp1); if(feof(fp1))goto l20;
  strcpy(f1,a); k=strlen(f1); f1[k-1]=0;
  fp=fopen(f1,"rb"); 
  fread(head,72,80,fp); fclose(fp);
  k=indexpos(head,"NAXIS1  ",72);
  sscanf(&head[k++][23],"%d",&n1);
  sscanf(&head[k++][23],"%d",&n2);
  if(n1!=4096 || n2!=4032){ printf("\nnot BOK_ccd !\n\n"); exit(0); }
  k=indexpos(head,"A81     ",72);
  if(k==72){ printf("\nfits_file not coordinated !\n"); exit(0); }
  for(i=0;i<8;i++)sscanf(&head[k+i][10],"%le",&a8[i]);
  k=indexpos(head,"CCD_NO: ",72);
  sscanf(&head[k][23],"%d",&i);

  x=y=0;     if(i==4)x=n1/2;
  xy_rad(a8[6],a8[7],x,y,&a1,&d1);
  x=0; y=n2; if(i==4){ x=n1/2; y=n2/2; }
  xy_rad(a8[6],a8[7],x,y,&a2,&d2);
  x=n1; y=0;   
  xy_rad(a8[6],a8[7],x,y,&a3,&d3);
  x=n1; y=n2; if(i==4)y=n2/2;
  xy_rad(a8[6],a8[7],x,y,&a4,&d4);
  x=n1/2; y=n2/2;
  xy_rad(a8[6],a8[7],x,y,&alpha,&delta);
  toms2(alpha,c1,1);  toms2(delta,c2,0);
//  printf("%s%s %s %10.7lf %10.7lf %10.7lf %10.7lf\n",f1,c1,c2,a3,a2,d1,d4);
  fprintf(fp2,"%s%s %s %10.7lf %10.7lf %10.7lf %10.7lf\n",f1,c1,c2,a3,a2,d1,d4);
  goto l10;
l20:
  fclose(fp1); fclose(fp2);
  system("rm v_1.tmp");
}

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

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