#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(int ac, char **av)
{
  char a[60],f1[60];
  double x,y,alpha,delta;
  double a1,a2,a3,a4,d1,d2,d3,d4;
  int i,k;
  int idate;
  if(ac<2){ printf("\n\t**** produce fitma.pos9 ***\n\n"); exit(0); }
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  fp1=fopen("ma.cal1","r");
  fp2=fopen("fitma.pos9","w");
l10:
  fgets(a,60,fp1); if(feof(fp1))goto l20;
  if(a[0]=='#')goto l10;
  sscanf(a,"%s",f1);
  idate=0; sscanf(&a[1],"%d",&k); if(k>5700)idate=1;
  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 ! %s\n",a); goto l10; }
  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(idate==0 && i==4)x=n1/2;
  xy_rad(a8[6],a8[7],x,y,&a1,&d1);
  x=0; y=n2; if(idate==0 && 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(idate==0 && 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);
  head[29][77]=0;
  fprintf(fp2,"%s%11.7lf%11.7lf%11.7lf%11.7lf%11.7lf%11.7lf%11.7lf%11.7lf %s\n",
          f1,a1,d1,a2,d2,a3,d3,a4,d4,&head[29][70]);
  goto l10;
l20:
  fclose(fp1); fclose(fp2);
}

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