#include <stdio.h>
FILE *fp;
char a[80],f1[60],f2[20],c1[14],c2[14],f10;
char head[72][80],c=37;
double u1,v1,x1,y1,x[4],y[4],ra,de,epoch;
int  w[4],p4=2,m1,m2,irm=1;
main(int ac,char **av)
{
  int i,j,k;
  if(ac<2){ printf("\n\t******* get real coordination ********\n");
            printf("\n\t                       jiang 2011,9,20\n"); 
            printf("\n\tUsage: getcoord d5534.0045.fits [!]\n");
            printf("\n\t! means do 4_ccds coord. else do 2,3 only(save time)");
            printf("\n\t(need program: sex  support)\n");
            printf("\n\tnext may: look4 e5534.0045");
            printf("\n\tnext may: c40 e5534.0045*.fit [!]");
            printf("\n\t------------\n\tif pointer error ~20arcmin");
            printf("\n\tlookall, another window run u3bok");
            printf("\n\tmove cursor in u3bok, center fit");
            printf("\n\t------------\n\tif pointer error ~80arcmin"); 
            printf("\n\tcoord60 e5534_?.fit");
            printf("\n\t------------");
            printf("\n\trun getcoord again\n\n");
            exit(0);
  }
  m1=1; m2=3; // do ccd_2, ccd_3 only, reason: save time, 
              // if want 4 ccd all work then set m1=0; m2=4;
  if(ac>2 && av[2][0]=='!'){ m1=0; m2=4; }
  for(i=0;i<4;i++)w[i]=0;
  fp=fopen(av[1],"rb");
  if(fp==0){ printf("\n\tinput file not found!\n"); exit(0); }    
  fclose(fp);
  strcpy(f1,av[1]); k=strlen(f1); f1[k-5]=0; 
// if e*_1.fit existed, not do d4,e4;
  sprintf(f2,"%s_1.fit",f1); f2[0]='e';
//  fp=fopen(f2,"rb"); if(fp){ fclose(fp); irm=0; goto l10; }  
//
  sprintf(a,"d4 %s",av[1]); system(a); // for display reason, produce all 4 e*
  for(i=0;i<4;i++){ sprintf(a,"e4 %s_%d.fit",f1,i+1); system(a); }
l10:
  f10=f1[0]; f1[0]='e';
  for(i=m1;i<m2;i++){
    sprintf(a,"coord8 %s_%d %c",f1,i+1,c); system(a); }
  for(i=m1;i<m2;i++){         
    sprintf(a,"%s_%d.fit",f1,i+1);
    fp=fopen(a,"rb"); if(fp==0)continue;
    fread(head,72,80,fp); fclose(fp);
//  write down getcoord.par
    if(i==m1){
      fp=fopen("getcoord.par","w");
      fprintf(fp,"%s\n",a);
      head[5][23]=head[6][23]=0;
      fprintf(fp,"%s %s\n",&head[5][11],&head[6][11]);  
      fclose(fp);
    }
    k=indexpos(head,"A87     ",72);
    if(k==72)continue; 
    w[i]=1;
    chs(&head[k][46],&x[i]);    head[k][67]=0;
    chs(&head[k+1][46],&y[i]);  head[k+1][67]=0;
                                head[k-7][30]=0;
    printf("%d: %s %s  seeing:%s\n",
      i+1,&head[k][46],&head[k+1][46],&head[k-7][24]);
  }
  if(((w[0]&w[3])|(w[1]&w[2]))==0){
    printf("Sorry, I cannot work out!! \t run coord60 pls!\n\n"); exit(0);
  }
  if((w[0]&w[3])==0)w[0]=w[3]=0;  
  if((w[1]&w[2])==0)w[1]=w[2]=0;  
  x1=y1=0.;  k=0;                   
  for(i=0;i<4;i++){         
    if(w[i]){ x1+=x[i]; y1+=y[i]; k++; }
  }
  x1/=k; y1/=k;
  if(w[1]&w[2]){ if(x[2]-x[1]>20.){ x1+=12.; if(x1>=24.)x1-=24.; } }
  else         { if(x[3]-x[0]>20.){ x1+=12.; if(x1>=24.)x1-=24.; } }
  toms1(x1,c1,1); toms1(y1,c2,0);
  for(i=0;i<60;i++)printf("-"); printf("\n");
  head[5][23]=head[6][23]=0;
  chs(&head[5][11],&u1);
  chs(&head[6][11],&v1);  
  k=indexpos(head,"DATE-OBS",72);
  for(i=0;i<30;i++){ c=head[k][i+8]; if(c<'0' || c>'9')c=32; a[i]=c; }
  sscanf(a,"%d %d %d",&i,&j,&k); epoch=i+(j-1)/12.+(k-1)/30./12.;
  k=epoch*10+0.5; epoch=k*.1;
  astprs2(u1,v1,2000.,&ra,&de,epoch);  toms1(ra,c1,1); toms1(de,c2,0);
  printf("OBS:  %s %s (2000.0)  %s %s (%6.1lf)\n",
  &head[5][11],&head[6][11],c1,c2,epoch);
  toms1(x1,c1,1); toms1(y1,c2,0);
  printf("CAL: %s  %s  (2000.0)",c1,c2);
  astprs2(x1,y1,2000.,&ra,&de,epoch);  toms1(ra,c1,1); toms1(de,c2,0);
  printf("  %s %s (%6.1lf)\n",c1,c2,epoch);

  x1-=u1; y1-=v1;
  x1*=3600.; y1*=3600.;  x1=-x1; y1=-y1;
  printf(" %12.2lf (%6.1lf  %6.1lf\n",x1,x1*15.,y1);
  sprintf(a,"rm %c*.fit",f10);
  if(irm)system(a);
//  system("rm wcs.coo");
  system("rm test.cat");
//  system("rm shift.bok");
}

chs(a,y)
char *a; double *y;
{
  double x;
  char b[20];
  int i,j,k;
  strncpy(b,a,20); 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;
}

toms1(aa,c,k)
double aa;           // input
char c[14];          // output
int k;               // input
{
  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:%5.2f",i,j,x);
  if(k==0)sprintf(&c[1],"%2.2d:%2.2d:%4.1f",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]=':';
  }
}

astprs2(ra1, dc1, ep1, ra2, dc2, ep2)   // RA in hour, DEC in degree
double  ra1,dc1,ep1,*ra2,*dc2,ep2;      // only ra2,dc2   is output
{
  double r0[3],r1[3],p[3][3],arc;
  double r2,d2;

  arc=45./atan(1.);
  *ra2=ra1; *dc2=dc1;
  if(ep1 == ep2)return;
  r2 =ra1*15./arc; d2 =dc1/arc;
  r0[0]=cos(r2)*cos(d2); r0[1]=sin(r2)*cos(d2); r0[2]=sin(d2);
  if(ep1 != 2000.){
    astrox(ep1, p);
    r1[0] = p[0][0] * r0[0] + p[0][1] * r0[1] + p[0][2] * r0[2];
    r1[1] = p[1][0] * r0[0] + p[1][1] * r0[1] + p[1][2] * r0[2];
    r1[2] = p[2][0] * r0[0] + p[2][1] * r0[1] + p[2][2] * r0[2];
    r0[0] = r1[0]; r0[1] = r1[1]; r0[2] = r1[2];
  }
  if(ep2 != 2000.){
    astrox(ep2, p);
    r1[0] = p[0][0] * r0[0] + p[1][0] * r0[1] + p[2][0] * r0[2];
    r1[1] = p[0][1] * r0[0] + p[1][1] * r0[1] + p[2][1] * r0[2];
    r1[2] = p[0][2] * r0[0] + p[1][2] * r0[1] + p[2][2] * r0[2];
    r0[0] = r1[0];    r0[1] = r1[1];    r0[2] = r1[2];
  }
  *ra2  = atan2(r0[1], r0[0])/15.*arc;
  *dc2 = asin(r0[2])*arc;
  if(*ra2<0)*ra2+=24.;
}

astrox (epoch, p)
double epoch,p[3][3];
{
  double t,a,b,c,ca,cb,cc,sa,sb,sc,arc;
  arc=45./atan(1.);
  astjuy(epoch,&t);
  t = (t - 2451545.0) / 36525.;
  a = t * (0.6406161 + t * (0.0000839 + t * 0.0000050));
  b = t * (0.6406161 + t * (0.0003041 + t * 0.0000051));
  c = t * (0.5567530 - t * (0.0001185 + t * 0.0000116));
  ca = cos (a/arc);
  sa = sin (a/arc);
  cb = cos (b/arc);
  sb = sin (b/arc);
  cc = cos (c/arc);
  sc = sin (c/arc);
  p[0][0] = ca * cb * cc - sa * sb;
  p[1][0] = -sa * cb * cc - ca * sb;
  p[2][0] = -cb * sc;
  p[0][1] = ca * sb * cc + sa * cb;
  p[1][1] = -sa * sb * cc + ca * cb;
  p[2][1] = -sb * sc;
  p[0][2] = ca * sc;
  p[1][2] = -sa * sc;
  p[2][2] = cc;
}

astjuy (epoch,t)
double epoch,*t;
{
  double jd;
  int year,centuy;
  year = epoch - 1;
  centuy = year / 100;
  jd =1721425.5+365.*year-centuy+ year/4 + centuy/4;
  year=epoch;
  *t= jd + (epoch - year) * 365.25;
}