#include <stdio.h>
#include <math.h>
#define  edge  600         // ~780
#define  err   1.0         // match star in 1 pixel. 1.3"
#define  m2    40000
#define  m1    m2*3
double a8[8],b8[8],pi,cx,cy;
char head[72][80],a[136],b[136],f1[80],d[80],c;
FILE *fp,*fp1;
double u[m1],v[m1],w[m1],s[m1];
float xx0[m2],yy0[m2],zz0[m2];
float u0[m2],v0[m2],w0[m2],w1[m2];
int m;
char a1[80],a2[80],mag2[20];

main(int ac, char **av)
{
  int i,j,k,jj,imag;
  double ra,de,x,y,z,t,x0,y0,z0,x1,y1;
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  if(ac<2){
    printf("\n\t ***** chang *.mag2 to X,Y file ******\n");
    printf("\n\tUsage: chxy t527i [!]               2011,3,3\n"); 
    printf("\n\toutput t527i.cat4, (new)t527i.coo (x,y adjusted)\n\n");
    exit(0);
  }
  findmag2(av[1]);
  sprintf(f1,"%s.fit",av[1]); c=f1[4];
  fp=fopen(f1,"rb"); if(fp==0){ printf("\n\t%s file not found!\n",f1);exit(0);}
  fread(head,72,80,fp); fclose(fp);
  k=indexpos(head,"A81     ",72);
  if(k==72){ printf("fits not coordination !\n"); exit(0); }
  for(i=0;i<8;i++)sscanf(&head[k+i][10],"%le",&a8[i]);
  xytoad(a8,b8);
  sprintf(f1,"%s.cat",av[1]);
  fp=fopen(f1,"r"); if(fp==0){ printf("4k_cat not found!\n");exit(0);}
  fgets(d,80,fp); m=0; ra=4096-edge;
l05:
  fgets(d,80,fp); if(feof(fp))goto l06;
  sscanf(&d[24],"%lf %lf %lf %lf",&x,&y,&z,&t);
  if(x<edge || y<edge)goto l05;
  if(x>ra   || y>ra)  goto l05;
  u[m]=x; v[m]=y; w[m]=z; s[m]=t; m++;
  goto l05;
l06:
  fclose(fp);
  j=k=jj=0; 
  fp=fopen(mag2,"r");
  fgets(a,136,fp);
  a[30]=0; if(ac==2)printf("%s\n",a);
  for(i=31;i<136;i++)if(a[i]==c)break; imag=i-4;
l10:
  fgets(a,136,fp); if(feof(fp))goto l20;
  fgets(b,136,fp);
  itohd(&a[7],&ra,&de);
  rad_xy(a8[6],a8[7],ra,de,&x,&y,0,0);
//  x-=0.09; y-=-0.13; test for o_filter
  a[29]=0; a[imag+6]=0; b[imag+6]=0;
  sscanf(&a[imag],"%lf",&z);
  de=err; w0[k]=w1[k]=0.;
  for(i=0;i<m;i++){
    x1=x-u[i]; y1=y-v[i];
    ra=x1*x1+y1*y1;
    if(ra<err){ de=ra; x0=x1; y0=y1; z0=z-w[i]; w0[k]=w[i]; w1[k]=s[i]; }
  }
  if(de==err)
    if(ac==2)
    printf("%s %8.2lf %7.2lf %s  %s\n",&a[6],x,y,&a[imag],&b[imag]);
  else { 
    if(ac==2)
    printf("%s %8.2lf %7.2lf %s  %s %5.2lf %5.2lf %7.3lf\n",
    &a[6],x,y,&a[imag],&b[imag],x0,y0,z0);
    xx0[j]=x0; yy0[j]=y0; if(z0>0.)zz0[jj++]=z0;
    j++;
  }
  u0[k]=x; v0[k]=y;
  k++;
  goto l10;
l20:
  fclose(fp);
  if(ac==2)
  printf("\t  2k.mag2            4k_xy           2k_mag       dx   dy    dm\n");
  printf("%d/%d      ",j,k);
  fmedian(xx0,j,&x0);
  fmedian(yy0,j,&y0);
  fmedian(zz0,jj,&z0);
  printf("file: %s    av_x0,y0: %6.2lf %6.2lf av_mag: %6.2lf\n",av[1],x0,y0,z0);
  for(i=0;i<k;i++)if(u0[i]!=0.){ u0[i]-=x0; v0[i]-=y0; }
  sprintf(f1,"%s.coo",av[1]);
  fp=fopen(f1,"r"); fgets(a1,80,fp); fgets(a2,80,fp); fclose(fp);
  fp=fopen(f1,"w"); fprintf(fp,"%s",a1); fprintf(fp,"%s\n",a2);
  for(i=0;i<k;i++)fprintf(fp,"%6d %8.3f %8.3f\n",i+1,u0[i],v0[i]);
  fclose(fp);
  sprintf(f1,"%s.cat4",av[1]);
  fp=fopen(f1,"w");
  for(i=0;i<k;i++)fprintf(fp,"%8.3f %8.3f\n",w0[i],w1[i]);
  fclose(fp);
}

findmag2(char *a)
{
  system("ls *.mag2 >mag2.tmp");
  fp=fopen("mag2.tmp","r"); if(fp==0)goto l01;
l00:
  fgets(mag2,40,fp); if(feof(fp))goto l01;
  if(mag2[6]=='k')goto l00;
  mag2[strlen(mag2)-1]=0;
  if(a[2]==mag2[2] && a[3]==mag2[3]) return;
l01:
  printf("%s?_??.mag2 not found !\n",a); exit(0); 
}
fmedian(a,n,xy)
float *a; int n; double *xy;
{
  int i,j;
  float z;
  for(i=0;i<n-1;i++)for(j=i+1;j<n;j++)if(a[i]>a[j]){
    z=a[i]; a[i]=a[j]; a[j]=z;
  }
  z=a[n/2]; if(n%2==0)z=(a[n/2]+a[n/2-1])*.5;
  *xy=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;
  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;}
}

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;
}
    
itohd(a,rc,dc)
char *a; double *rc,*dc;
{
  int i,j,k;
  int i1,i2,j1,j2;
  double x3,y3;
  char b[70];
  strcpy(b,a);
  for(i=0;i<24;i++)if(b[i]==':')b[i]=32;
  k=1; if(b[11]=='-'){ b[11]=32; k=-1; }
  sscanf(b,"%d %d %lf %d %d %lf",&i1,&i2,&x3,&j1,&j2,&y3);
//  printf("%2d %2d %5.2lf %2d %2d %4.1lf\n",i1,i2,x3,j1,j2,y3);
  x3=i1+i2/60.+x3/3600.;
  y3=j1+j2/60.+y3/3600.; if(k==-1)y3=-y3;
  *rc=x3;
  *dc=y3;
}