#include <stdio.h>
#include <math.h>
int    n,n2;                           // 2048,4096,   4064
double pi,cx,cy,sphi,cphi,s0,ut,cosdc;
double w_de, epoch;                 // output epoch
double a8[8],b8[6];                 // b8 is convert of a8
char   head[72][80],c1[14],c2[14];
FILE   *fp;
char   f1[60];
int    iy,im,id,itype,match;
int    iturn;                       // thick 2048 iturn=0 (include 1,3)
int    onum;                                    // thin  2048      =4
                                    // 4000,4096       =6 (include 7)
                                    // 4064            =8 ( 100707 define)
                                    // xuyi  4096      =9 ( not include here)
float  map[4096*4096];              // get skyadu, seeing
double rms,sky,seeing;
double xa,xd,xi,xn;                  // dummy
double oac,odc;
double oldac,olddc;

#define nn 29000
int    nstar,nu3;
double starx[nn],stary[nn],starv[nn], aa[nn],dd[nn];
short  ee[nn],pr2[nn],pd2[nn];                       // proper motion 
int    uno[nn];                               // number of u3
double a_c,d_c,up_a,dn_a,up_d,dn_d;  // always 2000.0,  for u3
double ha;
int    u3index[240][360];     
float    dra=0,dde=0;                  // for bok center adjust, in arcmin
int    ipic=0;
// struct { float x,y; short v; } v;

getstar(k)
int *k;
{
  char a[80];
  int i,j,ip=0;
  double z,z16=16.;                  // if more star, z16=16,15,14
  double edg;
  fp=fopen("default.sex","r");
  if(fp==0){
    system("cp /home/primefocus/jiang/vega2/rhbin/default.sex .");
    system("cp /home/primefocus/jiang/vega2/rhbin/default.conv .");
    system("cp /home/primefocus/jiang/vega2/rhbin/default.param .");
  } else fclose(fp);
  sprintf(a,"sex %s",f1);

  system(a);                         // max is nn
  edg=60.; if(n>=4000)edg*=2.;
l05:
  fp=fopen("test.cat","r");
  i=0; ip++;
l10:
  fgets(a,80,fp); if(feof(fp))goto l20;
  if(a[0]=='#')goto l10;
  sscanf(a,"%lf %lf %lf",&starv[i],&starx[i],&stary[i]);
  if(starv[i]>z16)goto l10;
  if(starx[i]<edg || starx[i]>n-edg )goto l10;
  if(stary[i]<edg || stary[i]>n-edg )goto l10;
  i++;  
  if(i<nn)goto l10;                 // nn=500
  else if(ip<3){ z16-=1.; fclose(fp);  goto l05; }
l20:
  if(i<50 && ip<2){ z16+=2.0; fclose(fp); goto l05; }
  fclose(fp);
  *k=i;
  for(i=0;i<*k-1;i++)for(j=i+1;j<*k;j++)if(starv[i]>starv[j]){
    z=starx[i]; starx[i]=starx[j]; starx[j]=z;
    z=stary[i]; stary[i]=stary[j]; stary[j]=z;
    z=starv[i]; starv[i]=starv[j]; starv[j]=z;
  }
// 2048_1
if(iturn<6)for(i=0;i<*k;i++){ starx[i]=n-starx[i]+1; stary[i]=n-stary[i]+1; }
}

double  fjd(int iy,int im,int id)
{
  int i,j;
  i=(int)((im+9)*0.09);
  i=(int)((i+iy-1900)*1.75+0.01);
  j=(int)(30.56*im);
  j=(iy-1950)*367+id+j-i;
  return(2433338.5+j);
}

getpar()
{
  int    i,j,ipos;
  double d1,d2,d,x,xk,x1;
  char   c[21];
  ipos=indexpos(head,"SKYADU  ",72);
  if(ipos!=72)sscanf(&head[ipos][12],"%lf",&sky); else sky=0.;
  ipos=indexpos(head,"DATE-OBS",72);
  strncpy(c,&head[ipos][10],20); for(i=0;i<21;i++)if(c[i]<'0'||c[i]>'9')c[i]=32;
  sscanf(c,"%d %d %d",&iy,&im,&id);
  if(iy<1900){ i=id; id=iy; iy=i+2000; if(iy>2050)iy-=100; }  // for SMT data
  ipos=indexpos(head,"TIME-OBS",72);
  if(ipos==72)ipos=indexpos(head,"TIME    ",72);              // for SMT time
  strncpy(c,&head[ipos][10],20); for(i=0;i<21;i++)if(c[i]<'.'||c[i]>'9')c[i]=32;
  sscanf(c,"%d %d %lf",&i,&j,&x);
  ipos=indexpos(head,"EXPTIME ",72);
  if(ipos==72)ipos=indexpos(head,"EXPOSURE",72);              // for SMT exp
  sscanf(&head[ipos][19],"%lf",&xk);
  ut=i+j/60.+x/3600.+xk/7200.;
  d1=fjd(1992,12,31);
  d2=(6+38/60.+40.1954/3600.)/24.;
  x1=-(111.+36.0/60.+1.6/3600.)/15.;       // 111:36:01.61 lamda of BOK
  if((n+n2)/2!=4064)x1=7+50/60.+18.344/3600.;     // SMT
  d=fjd(iy,im,id);
  d=(d-d1)*1.0027379093+d2;
  j=d; d-=j;
  s0=d*24.+x1+ut*1.0027379093;
  ipos=indexpos(head,"RA      ",72);
  strncpy(c,&head[ipos][10],20); for(i=0;i<21;i++)if(c[i]<'.'||c[i]>'9')c[i]=32;
  sscanf(c,"%d %d %lf",&i,&j,&x);
  d1=i+j/60.+x/3600.;
  ipos=indexpos(head,"DEC     ",72);
  strncpy(c,&head[ipos][10],20); xk=1.; for(i=0;i<21;i++)if(c[i]=='-')xk=-1.;
                                 for(i=0;i<21;i++)if(c[i]<'.'||c[i]>'9')c[i]=32;
  sscanf(c,"%d %d %lf",&i,&j,&x);
  d2=xk*(i+j/60.+x/3600.);
  ipos=indexpos(head,"EPOCH   ",72);
  if(ipos==72){ printf("I cannot found EPOCH !!!\n"); exit(0); }
  sscanf(&head[ipos][12],"%lf",&epoch);
  astprs2(d1,d2,epoch,&a_c,&d_c,2000.0);

    ipos=indexpos(head,"RA-OBS  ",72);
  strncpy(c,&head[ipos][10],20); for(i=0;i<21;i++)if(c[i]<'.'||c[i]>'9')c[i]=32;
    sscanf(c,"%d %d %lf",&i,&j,&x);
    oac=i+j/60.+x/3600.;
    ipos=indexpos(head,"DEC-OBS ",72);
    strncpy(c,&head[ipos][10],20); xk=1.; for(i=0;i<21;i++)if(c[i]=='-')xk=-1.;
                                 for(i=0;i<21;i++)if(c[i]<'.'||c[i]>'9')c[i]=32;
    sscanf(c,"%d %d %lf",&i,&j,&x);
    odc=xk*(i+j/60.+x/3600.);
  ipos=indexpos(head,"INSTRUME",72);
  iturn=0; if((n+n2)/2==4064)iturn=8;
  if(ipos!=72)for(i=11;i<20;i++){
    if(head[ipos][i]=='4')iturn=4;
    if(head[ipos][i]=='5')iturn=5;
    if(head[ipos][i]=='6' || head[ipos][i]=='7')iturn=6;
    if(head[ipos][i]=='8')iturn=8;
  }  

  if(iturn<6)w_de=n/2048.*58.;              // for old 2048, wjh
  if((n+n2)/2==4064)w_de=30.;
  if(iturn==6)w_de=n/4096.*90.;             // for E2V, wjh
  x=(31.+57./60.+46.5/3600.)*cy;            // 31:57:46.5 BOB latitude
  if((n+n2)/2!=4064)x=(40.+23./60.+36./3600.)*cy;  // SMT
  sphi=sin(x); cphi=cos(x);
}

//     rc,dc is 2000.  in hour,degree
//     aa,dd is 2000.0 + proper_motion(epoch)
fgsc(rc,dc,w_de)
double rc,dc,w_de;
{
  double hw_de,hw_al,x,y,z;
  int    i,j,k1,k2;
  short  s;
  fp=fopen("/home/primefocus/jiang/ucac3/u3index2.unf","rb");
  fread(u3index,360*240,4,fp); fclose(fp);
  hw_de=w_de/120.;
  up_d=dc+hw_de;  
  dn_d=dc-hw_de; 
  x=cos(up_d*cy);   y=cos(dn_d*cy);  z=x; if(z>y)z=y;
  hw_al=hw_de/z/15.; // if(hw_al>12.)hw_al=12.;
  up_a=rc+hw_al;      if(up_a>=24.)up_a-=24.;
  dn_a=rc-hw_al;      if(dn_a<  0.)dn_a+=24.;
  cosdc=cos(dc*cy); if(cosdc<.05)cosdc=.05;
  nu3=0;
  k1=(dn_d+90.)*2.+1.;  k2=(up_d+90.)*2.+0.9999; 
  for(i=k1;i<=k2;i++)cat_u3(i);  //  printf("%03d: u3_star--> %d\n",i,nu3); }
  for(i=0;i<nu3-1;i++)for(j=i+1;j<nu3;j++){
    x=ee[i]; if(x<0.)x=-x;
    y=ee[j]; if(y<0.)y=-y;
    if(x>y){
      x=aa[i]; aa[i]=aa[j]; aa[j]=x;
      x=dd[i]; dd[i]=dd[j]; dd[j]=x;
      s=ee[i]; ee[i]=ee[j]; ee[j]=s;
      s=pr2[i]; pr2[i]=pr2[j]; pr2[j]=s;
      s=pd2[i]; pd2[i]=pd2[j]; pd2[j]=s;
      k1=uno[i]; uno[i]=uno[j]; uno[j]=k1;
    }
  }
}

struct { int ra, dc; short mag, pr2, pd2; } u;
cat_u3(i360)
int i360;
{
  double x1,x2, ra,dc,mas,xmag;
  char   f1[20];
  int    i,j,k,i87;
  mas=3600000.;    i=i360-1;
  x1=dn_a; x2=up_a;  if(x2<x1)x2+=24.;
  j=x1*10.; j--; if(j<0)k=0; else  k=u3index[j][i];
  i87=0; // if(up_d>87. || dn_d<-87.){ i87=1;  k=0; }
  sprintf(f1,"/home/primefocus/jiang/ucac3/u%03d\0",i360);
  fp=fopen(f1,"rb");
l05:
  fseek(fp,k*14,0);
l10:
  fread(&u,1,14,fp);   if(feof(fp))goto l20;
  k++;
  ra=u.ra/mas/15.;
  if(ra>x2 && i87==0)goto l100;               // normal exit
  if(ra<x1 && i87==0)goto l10;
  dc=(u.dc-mas*90.)/mas;  if(dc<dn_d || dc>up_d)goto l10;
  ra+=u.pr2*.1*(epoch-2000.0)/mas/15./cosdc;  // ra,dc is 2000's data + pm
  dc+=u.pd2*.1*(epoch-2000.0)/mas;
  aa[nu3]=ra;
  dd[nu3]=dc;
  ee[nu3]=u.mag;
  pr2[nu3]=u.pr2;
  pd2[nu3]=u.pd2;
  uno[nu3]=i360*1000000+k;
  nu3++; if(nu3==nn)nu3--;                    // becouse of read more cat
  goto l10;
l20:
  x1-=24.; x2-=24.;  k=0;
  if(i87==0)goto l05;
l100:
  fclose(fp);
}

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

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

standc(rc,dc,ra,de,xi,xn)
double rc,dc,ra,de,*xi,*xn;
{
  double sd,cd,td,co,tmp;
  sd=sin(dc);
  cd=cos(dc);
  td=tan(de);
  co=cos(ra-rc);
  tmp=sd*td+cd*co;
  *xi=sin(ra-rc)/tmp;
  *xn=(cd*td-sd*co)/tmp;
}

astand(rc,dc,xi,xn,ra,de)
double rc,dc,xi,xn,*ra,*de;
{
  double cd,td;
  cd=cos(dc);
  td=tan(dc);
  *ra=atan(xi/cd/(1.-xn*td));
  *de=atan((xn+td)*cos(*ra)/(1.-xn*td));
  *ra+=rc;
  cd=(*de)-dc; if(cd<0.)cd=-cd;
  if(cd>2.){  *de=-(*de);  *ra+=pi; }
  if(*ra<0.)*ra+=(pi+pi);
}

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

main_job()
{
  char   c3[3];
  char   ns1[50],dis[60];
  double x9[99],y9[99];
  double sx[4],sy[4];
  int    i,j,k=0;
  float  x,y,u,v;
  int    ipp=0;

    a_c=oac+(dra/900.);
    d_c=odc+(dde/60.);
//    if(iturn==8){            // 4064 BOK
      a8[0]=-0.2193e-05;
      a8[1]= 0.2184e-07;
      a8[2]= a8[1];
      a8[3]=-a8[0];
      a8[4]= 0.4457e-02;
      a8[5]=-0.4475e-02;
//    }
    a8[6]=a_c*cx;
    a8[7]=d_c*cy;

    k=indexpos(head,"CCD_NO: ",72);
    k=head[k][29]-48;
      a8[0]=-0.2193e-05;
      a8[1]= 0.2184e-07;
      a8[2]= a8[1];
      a8[3]=-a8[0];
      a8[4]= 0.4457e-02;
      a8[5]=-0.4475e-02;
    u=0.018557/cos(d_c*cy);
    v=0.265689;
    x=1.33/3600./cos(d_c*cy);
    y=27./3600.;
    if(k==2){ a_c+=u; d_c-=v; }
    if(k==3){ a_c-=u; d_c+=v; }
    if(k==1){ a_c+=(u+x); d_c+=(v-y); }
    if(k==4){ a_c-=(u+x); d_c-=(v-y); }
    
    a8[6]=a_c*cx; a8[7]=d_c*cy;
    oldac=a8[6]; olddc=a8[7];
    itype=0;
  fgsc(a_c,d_c,w_de*1.1);
  if(itype==0){
l15:  
    ipp++;
    rms=55.; match=kgsc(99);
    if(match==0){
      shift_c();
      if(match<5){ if(ipp==2)return(0);  goto l15; }
    }
    a_c=a8[6]/cx; d_c=a8[7]/cy;
    fgsc(a_c,d_c,w_de);
  }
  rms=5.;  match=kgsc(0);
  if(match==0)return(0); 
  if(a8[6]<0.)a8[6]+=pi+pi;
  return(1);
}

double xx[nn],yy[nn],vv[nn],wa[nn],wd[nn],gr[nn],gd[nn],xixn[nn][2],xy[nn][2];
short  zz[nn];
char   xyw[nn];
int kgsc(kk)
int kk;
{
  int i,j,k,m,l,n12,nu,ns,iz;
  double x,z1,z2,z,sig,sigma;
  float  fpeek,fsigma;
  ns=nstar; nu=nu3; if(kk>0)ns=nu=kk;
  if(nstar<kk)ns=nstar; if(nu3<kk)nu=nu3;
//  printf("ccd_star:%d  u3_star:%d\n",ns,nu);
  for(i=0;i<ns;i++){ xx[i]=starx[i]; yy[i]=stary[i]; vv[i]=starv[i]; xyw[i]=0;}
  n12=10; if(ns<=kk)n12=6;
  xytoad(a8,b8);
  for(i=0;i<nu;i++){
    xa=aa[i]*cx; xd=dd[i]*cy;
    standc(a8[6],a8[7],xa,xd,&xi,&xn);
    wa[i]=b8[0]*xi+b8[2]*xn+b8[4];
    wd[i]=b8[1]*xi+b8[3]*xn+b8[5];
  }

  for(i=0;i<nu;i++){
    x=rms;
    for(j=0;j<ns;j++){
      z1=wa[i]-xx[j];
      z2=wd[i]-yy[j];
      z=sqrt(z1*z1+z2*z2);
      if(z<x){ k=j; x=z; }
    }
    if(x!=rms){ xyw[k]=1; gr[k]=aa[i]; gd[k]=dd[i]; }
  }
  j=0;
  for(i=0;i<ns;i++)if(xyw[i]){
    xy[j][0]=xx[j]=xx[i];
    xy[j][1]=yy[j]=yy[i];
             vv[j]=vv[i];
    gr[j]=gr[i]; gd[j]=gd[i];
    xyw[j]=1;
    xa=gr[j]*cx; xd=gd[j]*cy;
    standc(a8[6],a8[7],xa,xd,&xi,&xn);
    xixn[j][0]=xi;
    xixn[j][1]=xn;
    j++;   
  }
  l=j;
  k=0;
l151:
  if(l-k<n12)return(0);
  k++;
  i=3; plate_(xixn,xy,xyw,&l,a8,&i);
  iz=m=0; sig=z=0.;
  for(i=0;i<l;i++){
    xi=a8[0]*xy[i][0]+a8[2]*xy[i][1]+a8[4];
    xn=a8[1]*xy[i][0]+a8[3]*xy[i][1]+a8[5];
    astand(a8[6],a8[7],xi,xn,&xa,&xd);
    xa/=cx; xd/=cy;
    sigma=gr[i]-xa; if(sigma<0.)sigma=-sigma;
    if(sigma>20.)sigma-=24.;
    sigma=sigma*15.*cos(a8[7]);
    sigma=sigma*sigma+(gd[i]-xd)*(gd[i]-xd);
    sigma=sqrt(sigma)*3600.;
    if(xyw[i]==1){                           // xyw=0.1.-1
      zz[iz++]=sigma*10.+0.5;
      if(sigma>z){ z=sigma; j=i; }
      sig+=sigma; m++;
    }
  }
  x=sig/m;
//  if(x>6. && k>9)return(0);
  if(z>5.0){ xyw[j]=-xyw[j]; goto l151; }    // coord RMS precison   ???
  k=-k;
  z=n/2.; if(iturn<6)z=(n+1)/2.;         // 2048_4
  z2=n2/2.; if(iturn<6)z2=(n2+1)/2.;         // 2048_4
  xi=a8[0]*z+a8[2]*z2+a8[4];              // new center RA,DEC
  xn=a8[1]*z+a8[3]*z2+a8[5];              // center X,Y is n/2
  astand(a8[6],a8[7],xi,xn,&xa,&xd);     // due to history reason, =(n+1)/2
  a8[6]=xa; a8[7]=xd;
  for(i=0;i<l;i++){
    xa=gr[i]*cx;
    xd=gd[i]*cy;
    standc(a8[6],a8[7],xa,xd,&xi,&xn);
    xixn[i][0]=xi;
    xixn[i][1]=xn;
  }
  if(k<0)goto l151;
  if(iz>20){
    white_black(zz,1,iz,&fpeek,&fsigma);
    x=fpeek*0.1;
  } else return (0);
  printf("\nBy using %d stars to calculate 8 ceof, RMS:%6.2lf arcsec\n",iz,x);
  rms=x;
  return(iz);
}

double xmedian(a,m)
double *a; int m;
{
  int i,j;
  double x;
  for(i=0;i<m-1;i++)for(j=i+1;j<m;j++)if(a[i]>a[j]){
    x=a[i]; a[i]=a[j]; a[j]=x;
  }
  j=m/2;
  if(m%2)return(a[j]);
  else return((a[j-1]+a[j])/2.);
}

double adx[99],ady[99],czr1[99],czr2[99],stx[99],sty[99];
int    cz1[99],cz2[99];
short  w[199];
shift_c()              // use bright 30 object
{
  int    i,j,ip=0,m,n30;
  double r,rr;
  float peak,s;
  n30=99;
  for(i=0;i<n30;i++){                // 2048_7
    stx[i]=starx[i]; sty[i]=stary[i];
  }
l5:
  xytoad(a8,b8);
  for(i=0;i<n30;i++){
    xa=aa[i]*cx;
    xd=dd[i]*cy;
    standc(a8[6],a8[7],xa,xd,&xi,&xn);
    adx[i]=b8[0]*xi+b8[2]*xn+b8[4];
    ady[i]=b8[1]*xi+b8[3]*xn+b8[5];
  }
  m=0;
  for(i=0;i<n30;i++){
    r=500.;                         // find minum 
    cz1[i]=0;
    for(j=0;j<n30;j++){
      xa=adx[i]-stx[j]; xd=ady[i]-sty[j];
      rr=sqrt(xa*xa+xd*xd);
      if(rr<r){ r=rr; cz1[i]=j; czr1[i]=r; }
    }
    r=500.;                         // find 2nd minum
    cz2[i]=0;
    for(j=0;j<n30;j++){
      xa=adx[i]-stx[j]; xd=ady[i]-sty[j];
      rr=sqrt(xa*xa+xd*xd);
      if(rr<r && rr!=czr1[i]){ r=rr; cz2[i]=j; czr2[i]=r; }
    }
    if(cz1[i]){ m++; w[m]=czr1[i]*10.; }
    if(cz2[i]){ m++; w[m]=czr2[i]*10.; }
  }
  white_black(w,m,1,&peak,&s);      // use 1st,2nd minum to solve
//    for(i=0;i<m;i++)printf("%d ",w[i]);
  peak*=0.1;
//    printf("peak: %f\n",peak);             // should be a small value
  s=3.;
  m=0;
  for(i=0;i<n30;i++){
    if(cz1[i]==0 && cz2[i]==0)continue;
l70:
    if(cz1[i]==0){
      r=czr2[i]-peak; if(r<0.)r=-r;
      if(r>s)continue;
      j=cz2[i];
      m++;
      czr1[m]=adx[i]-stx[j];
      czr2[m]=ady[i]-sty[j];
      continue;
    }
    if(cz2[i]==0){
      r=czr1[i]-peak; if(r<0.)r=-r;
      if(r>s)continue;
      j=cz1[i];
      m++;
      czr1[m]=adx[i]-stx[j];
      czr2[m]=ady[i]-sty[j];
      continue;
    }
    xa=czr1[i]-peak; if(xa<0.)xa=-xa;
    xd=czr2[i]-peak; if(xd<0.)xd=-xd;
    if(xa>xd)cz1[i]=0;
    if(xd>xa)cz2[i]=0;
    goto l70;
  }
  if(m<5 && ip==0){
    ip=1; goto l5;
  }
  xa=xmedian(czr1,m);
  xd=xmedian(czr2,m);
  s=-sqrt(xa*xa+xd*xd)*a8[0]*10800/pi;
//  printf(" Auto shift %5.1lf arc_minutes\n",s);
  a8[6]+=xa*a8[0]/cos(a8[7]); 
  if(iturn!=8)a8[7]+=xd*a8[0];
  else a8[7]-=xd*a8[0];                   // BOK is -
}


int oa[1600],ob[1600],oc[1600],od[1600];
runnum()
{
  int i,j,k=0,m=2,l;
  for(i=0;i<80;i+=m)for(j=0;j<80;j+=m)
     { oa[k]=i; ob[k]=j; oc[k]=i+j; od[k]=k++; }
  for(i=0;i<k-1;i++)for(j=i+1;j<k;j++)if(oc[i]>oc[j]){
    l=oa[i]; oa[i]=oa[j]; oa[j]=l;
    l=ob[i]; ob[i]=ob[j]; ob[j]=l;
    l=oc[i]; oc[i]=oc[j]; oc[j]=l;
    l=od[i]; od[i]=od[j]; od[j]=l;
  }  
  return(k);
}

int oi[4]={1,1,-1,-1};
int oj[4]={1,-1,1,-1};
main(ac,av)
int ac; char *av[];
{
  int  i,j,k;
  double x,y,z;
  char a[80],c;
  pi=4.*atan(1.0); cx=pi/12.; cy=pi/180.;
  if(ac<2){ 
            printf("\n\t  ****** auto find coordination in 80arcmin\n");
            printf("\n\tusage: coord60 e*_?.fit\n\n");
            printf("\n\tbefore run: getcoord\n");
            printf("\n\tafter  run: getcoord\n\n");
            exit(0);
          }
//  check file
  strcpy(f1,av[1]);
  if(f1[0]!='e'){ printf("\n\tinput file 1ST must with char--e\n\n"); exit(0); }
//   j=nindex(f1,".");
  k=strlen(f1); c=k-5; if(c<0)c=0; j=0; for(i=c;i<k;i++)if(f1[i]=='.')j=i;  
  if(j==0) strcat(f1,".fit");
  fp=fopen(f1,"rb"); if(fp==0){ printf("\t%s file not found!\n",f1); exit(0); }
  fread(head,72,80,fp);
  k=indexpos(head,"NAXIS1  ",72);  sscanf(&head[k][24],"%d",&n); 
                                   sscanf(&head[k+1][24],"%d",&n2);
//  read rest of head
  k=indexpos(head,"END     ",72);
  if(k==72)fread(map,36,72,fp);
//  check itype
  i=indexpos(head,"A81     ",72);
  itype=1;     if(i==72)itype=0;          // not found A81
  if(itype==1 && ac>2)itype=2;
//  if(itype==0)newhead();                // produce new fits_head
  if(itype==1){ printf("\tcoordinating already done !\n"); exit(0); }
  getpar();                               // get s0,iy,im,id,ut,iturn,w_de,phi
  fread(map,n*n2,4,fp); swap4(map,n*n2*4); 
  fclose(fp);                             // for skyadu, seeing only

  getstar(&nstar);
  onum=runnum();
  for(i=1;i<onum;i++){
//    printf("%d_%d  ",oa[i],ob[i]);
    printf("."); fflush(stdout);
    for(j=0;j<4;j++){
      dra=oa[i]*oi[j];  dde=ob[i]*oj[j];
      if(main_job())goto l999;
    }
  }
l999:
  if(i==onum){ printf(" ********* Fail **********\n"); exit(0); }    

  xytoad(a8,b8);
//  a_c=a8[6]/cx; d_c=a8[7]/cy;  // 2000, center in HD 
  xa=n/2.; xd=n2/2.;
  xi=a8[0]*xa+a8[2]*xd+a8[4];
  xn=a8[1]*xa+a8[3]*xd+a8[5];
  astand(a8[6],a8[7],xi,xn,&xa,&xd);
  a_c=xa/cx; d_c=xd/cy;

    oldac=(a8[6]-oldac)/cy*60+dra;
    olddc=(a8[7]-olddc)/cy*60+dde;
    printf("shift:%6.1lf%6.1lf   (arcmin)\n",oldac,olddc);
    fp=fopen("shift.bok","w"); 
    fprintf(fp,"%6.1lf%6.1lf\n",oldac,olddc); fclose(fp); 
    printf("\n\t**********ok!  shift.bok  updated !\n");
    printf("\n\t\t\t  Now, you can run getcoord again! ***********\n\n");
}