#include <stdio.h>
#include <math.h>
int    n1,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    ik;       // INSTRUME  [29],1-4// 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
char   pgmark;

#define nn 50000
int    nstar,nu3;
double starx[nn],stary[nn],starv[nn], aa[nn],dd[nn];
short  ee[nn],pr2[nn],pd2[nn],nof[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][720];     
// struct { float x,y; short v; } v;
#define bok  48.
#define btype 2

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 /vega2/rhbin/default.sex .");
    system("cp /vega2/rhbin/default.conv .");
    system("cp /vega2/rhbin/default.param .");
  } else fclose(fp);
  sprintf(a,"sex %s",f1);

  system(a);                         // max is nn
  edg=6.; if(n1>=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]>n1-edg )goto l10;
  if(stary[i]<edg || stary[i]>n2-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;
  }
}

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,"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((n1+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,"INSTRUME",72);
  ik=head[ipos][29]-48;
  w_de=31.;
  x=(31.+57./60.+46.5/3600.)*cy;            // 31:57:46.5 BOB latitude
  if((n1+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("/PPMXL/ppmindex.unf","rb");
  fread(u3index,720*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.)*4.;  k2=(up_d+90.)*4.+0.9999; 
  for(i=k1;i<k2;i++){ cat_u3(i); printf("%03d: ppm_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,xa,xd,xmag;
  char   f1[30];
  int    i,j,k,i87;
  mas=3600000.;    i=i360;
  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; }
  if(i360<360)sprintf(f1,"/PPMXL/s%02d%c\0",89-i360/4,100-(i360%4));
  else     sprintf(f1,"/PPMXL/n%02d%c\0",(i360-360)/4, 97+(i360%4));
  printf("%s: ",&f1[7]);
  fp=fopen(f1,"rb");
  if(fp==0){ printf("datafile not found! %s\n",f1); exit(0); }
l05:
  fseek(fp,k*14,0);
l10:
  fread(&u,1,14,fp);   if(feof(fp))goto l20;
  k++;
  if(u.mag==0 || u.mag>19000)goto l10;
  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;
  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]=k;
  nu3++; if(nu3==nn)nu3--;
  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;
}

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;
  int    ipp=0;

// ******** auto match
  printf("wait, auto match stars. Now, pick out CCD stars\n");
//  getstar(nstar);  read in file uedo.coo from "sex"
  getstar(&nstar);
//for(i=0;i<nstar;i++)printf("%d %lf %lf %lf\n",i+1,starx[i],stary[i],starv[i]);
  if(itype){
    k=indexpos(head,"A81     ",72);
    for(i=0;i<8;i++)sscanf(&head[k+i][10],"%lf",&a8[i]);
  } 
  xytoad(a8,b8);
  a_c=n1/2; d_c=n2/2;  
  xy_rad(a8[6],a8[7],a_c,d_c,&xa,&xd,btype,1);
  fgsc(xa/cx,xd/cy,w_de);
// shift center
  fp=fopen("coord.cen","r"); fgets(dis,60,fp); fclose(fp);
  chs(dis,&a_c);  chs(&dis[14],&d_c);
//  a_c=49./60.+59.421/3600.;
//  d_c=42.+31./60.+18.25/3600.;
  ad_xy(a_c*cx,d_c*cy,&xa,&xd);
  b8[4]=xa; b8[5]=xd; xytoad(b8,a8);
  a8[6]=a_c*cx;       a8[7]=d_c*cy;  
  rms=100.;  match=kgsc(0);
  if(match==0){ printf("match fail (star too faint) !\n"); exit(0); }
}

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


ad_xy(ra,de,x,y)
double ra,de,*x,*y;
{
  double sd,cd,td,co,tmp, rc,dc,xi,xn;
  rc=a8[6]; dc=a8[7];
  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;
  *x=b8[0]*xi+b8[2]*xn+b8[4];
  *y=b8[1]*xi+b8[3]*xn+b8[5];
}


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  ppm_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++)rad_xy(a8[6],a8[7],aa[i],dd[i],&wa[i],&wd[i],btype,0);
/*
    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];
  }
*/
//  printf("gsc---------\n");
//  for(i=0;i<nu;i++)printf("%9.1f %9.1f\n",wa[i],wd[i]);
//  printf("ccd---------\n");
//  for(i=0;i<ns;i++)printf("%9.1f %9.1f\n",xx[i],yy[i]);
//  for(i=0;i<10;i++)printf("%f %f %f %f\n",xx[i],yy[i],wa[i],wd[i]);

  float x1,x2;

  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;
    rad_xy(a8[6],a8[7],gr[j],gd[j],&xi,&xn,btype,1);
    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++){
    xy_rad(a8[6],a8[7],xy[i][0],xy[i][1],&xa,&xd,btype,0);
/*
    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>0.5){ xyw[j]=-xyw[j]; goto l151; }    // coord RMS precison
  k=-k;
  for(i=0;i<l;i++){
    rad_xy(a8[6],a8[7],gr[i],gd[i],&xi,&xn,btype,1);
/*
    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;
  }
  printf("\tBy using %d stars to calculate 8 ceof, RMS:%6.2lf arcsec\n",iz,x);
  rms=x;

  if(pgmark!='!' && pgmark!='%'){
    i=0; j=1; pgbegin_(&i,"/xw",&j,&j,3);
    x1=0; x2=n2; pgenv_(&x1,&x2,&x1,&x2,&j,&j);
    for(i=0;i<nu;i++){
      k=1; if(xyw[i]==1)k=2; pgsci_(&k); k=22;
      x1=wa[i]; x2=wd[i]; pgpoint_(&j,&x1,&x2,&k);
    }
    k=3; pgsci_(&k); k=1;
    for(i=0;i<ns;i++){
      x1=xx[i]; x2=yy[i]; pgpoint_(&j,&x1,&x2,&k);
    }
    pgend_();  
  }

// for wcs
  if(kk==0){
    fp=fopen("wcs.coo","w");
    fprintf(fp,"# %d*%d\n",n1,n2);
    for(i=0;i<l;i++)
    if(xyw[i]==1)fprintf(fp,"%8.3lf %8.3lf %5.2lf\n",xx[i],yy[i],vv[i]);
    fclose(fp);
  }
  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[30],ady[30],czr1[30],czr2[30],stx[30],sty[30];
int    cz1[30],cz2[30];
short  w[60];
shift_c()              // use bright 30 object
{
  int    i,j,ip=0,m,n30;
  double r,rr;
  float peak,s;
  n30=30;
  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++)rad_xy(a8[6],a8[7],aa[i],dd[i],&adx[i],&ady[i],btype,0);
/*
    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]); 
  a8[7]-=xd*a8[0];                   // BOK is -
}

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\tusage: poord81 filename ! (jiang101006)\n");
            printf("\n\t 1+bok*(xi**2+xn**2) %f\n",bok);
            printf("\n\tuse file coord.cen, shift center\n\n"); 
            exit(0);
          }
//  check file
  strcpy(f1,av[1]);
//   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",&n1); 
                                   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(ac>2)pgmark=av[2][0];
//  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,w_de,phi
  fread(map,n1*n2,4,fp); swap4(map,n1*n2*4); 
  fclose(fp);                             // for skyadu, seeing only
  
  main_job();

// get sky, seeing       used starx,stary,nstar
  sky_seeing(map,&sky,&seeing);

  k=32;                       // clear head
  for(i=k;i<71;i++)for(j=0;j<80;j++)head[i][j]=32;
  k=indexpos(head,"CRVAL1  ",72);
  if(k!=72)for(i=k;i<k+4;i++)for(j=0;j<47;j++)head[i][j]=32;
  k=indexpos(head,"CTYPE1  ",72);
  if(k!=72)for(i=k;i<k+4;i++)for(j=0;j<50;j++)head[i][j]=32;
  put_a8();               // 49-58 write sky seeing,a8
  put_xmass();            // 15-17, 64 write center,h_angle,airmas
  put_moon();             // 67-71  phase, position, pos_angle
  put_ds9();              // 40-48
  
  for(i=0;i<72;i++)for(j=0;j<80;j++)if(head[i][j]==0)head[i][j]=32;
  printf(" re_write fits head !\n");  
  fp=fopen(f1,"rb+");
  fwrite(head,80,72,fp);
  fclose(fp);
/*
  if((n1+n2)/2!=4064)exit(0);
// wcs
  printf("wait...doing WCS... ");
  sprintf(a,"(imwcs -c ucac3 -v -d wcs.coo -h 200 -w %s>wcs.res)>&wcs.log",f1);
  system(a);
  printf("ok !\n");
// produce u3.cat
  fp=fopen("u33.cat","w");
  astprs2(a8[6]/cx,a8[7]/cy,2000.0,&x,&y,epoch);
  toms2(x,c1,1); toms2(y,c2,0);
  fprintf(fp,"#center:%s %s (%6.1lf) %s  pm(mas)/year\n#",c1,c2,epoch,f1);
  for(i=0;i<66;i++)fprintf(fp,"-"); fprintf(fp,"\n");
  for(i=0;i<nu3;i++){
    astprs2(aa[i],dd[i],2000.0,&x,&y,epoch);
    toms2(x,c1,1); toms2(y,c2,0);
    z=ee[i]*0.001; c='s'; if(z<0.){ z=-z; c='g'; }
    j=uno[i]; k=j/1000000; j-=k*1000000;
    fprintf(fp,"%s %s (%6.1lf) %03d-%06d %c %6.3lf%6.1f%6.1f\n",
        c1,c2,epoch,k,j,c,z,pr2[i]*0.1,pd2[i]*0.1);
  } fclose(fp);
*/
}

sky_seeing(map,sky,seeing)
float map[n2][n1];
double *sky,*seeing;
{
  int ix,iy, i,j,k,m, b[113];
  short d[300*300];
  float z,s1,s2;
                                 m=0; 
  *sky=0.;                    // use 4 block to get av_sky
  for(ix=n1/4;ix<n1;ix+=n1/2)for(iy=n1/4;iy<n1;iy+=n1/2){
  k=0;for(i=ix-150;i<ix+150;i++)for(j=iy-150;j<iy+150;j++)d[k++]=map[i][j]/2.;
    white_black(d,1,k,&z,&s1);   if(z==0.)m++;
    *sky+=z;
  } *sky/=2.;                    if(m)*sky*=4.;

  for(m=0;m<nstar;m++){
    ix=starx[m]+0.5; iy=stary[m]+0.5;
    for(j=0,i=ix-14;i<=ix+14;i++,j+=4){
    z=0.; for(k=-2;k<=2;k++)z+=map[iy+k][i]-*sky; b[j]=z; }
    for(j=0;j<112;j+=4)b[j+2]=(b[j]+b[j+4])/2.+0.5;
    for(j=0;j<112;j+=2)b[j+1]=(b[j]+b[j+2])/2.+0.5;
    for(j=0;j<112;j++)if(b[j]<0)b[j]=0;
    histat(b,&z,&s1,112);
    for(j=0,i=iy-14;i<=iy+14;i++,j+=4){
    z=0.; for(k=-2;k<=2;k++)z+=map[i][ix+k]-*sky; b[j]=z; }
    for(j=0;j<112;j+=4)b[j+2]=(b[j]+b[j+4])/2.+0.5;
    for(j=0;j<112;j+=2)b[j+1]=(b[j]+b[j+2])/2.+0.5;
    for(j=0;j<112;j++)if(b[j]<0)b[j]=0;
    histat(b,&z,&s2,112);
    d[m]=(s1+s2)/8.*2.355*100.+0.5;          // expand 100
  }
  k=0; for(i=0;i<m;i++)if(d[i]>200 && d[i]<999)d[k++]=d[i];
  white_black(d,1,k,&z,&s1); z*=0.01;        // shrink 100
  *seeing=z*0.45;
  printf("SKY:%9.2lf  seeing:%6.2lf   ",*sky,*seeing);
}

double sunlamda(jd)                   // in degree
double jd;
{                                     // sun lamda in degree, beta==0
  double y,g,x;
  int j;
  x=jd-2451545.;
  y=280.460+0.9856474*x;
  g=(357.528+0.9856993*x)/360.;
  j=g; g-=j;
  x=g*360.*cy;
  x=(y+1.915*sin(x)+0.02*sin(x+x))/360.;
  j=x; x-=j;
  return(x*360.);
}

getmxy(phase,pl,pb,xx)
double *phase,*pl,*pb,*xx;
{
  double a1,d1,a2,d2,x,y,jd;
  astprs2(a_c,d_c,2000.,&xa,&xd,epoch);
  ad_azl(xa,xd,&a1,&d1);
  jd=fjd(iy,im,id)+ut/24.;
  moon2(jd,&a2,&d2);                          // get monn l,b
  lb_ad(a2,d2,&x,&y);                         // get moon alpha delta
  a2-=sunlamda(jd); if(a2<0.)a2+=360.; 
  *phase=a2/360.*29.53;                       // get moon phase
  ad_azl(x,y,&a2,&d2);                        // get moon A, H
  *pl=a2;
  *pb=d2;
  x=sin(d1*cy)*sin(d2*cy)+cos(d1*cy)*cos(d2*cy)*cos(a1*cy-a2*cy);
  y=acos(x);
  *xx=y/cy;
}

azl_ad(a,h,alpha,delta)                     // dd->hd
double a,h,*alpha,*delta;
{
  double aa,hh,sina,cosa,sinh,cosh,sinde,cosde,sint,cost,t;
  aa=a*cy;
  hh=h*cy;
  sina=sin(aa);
  cosa=cos(aa);
  sinh=sin(hh);
  cosh=cos(hh);
  sinde=sphi*sinh-cphi*cosa*cosh;
  cosde=sqrt(1.-sinde*sinde);
  sint=cosh*sina/cosde;
  cost=(cosa*cosh+cphi*sinde)/(sphi*cosde);
  *delta=atan(sinde/cosde)/cy;
  t=atan2(sint,cost);
  *alpha=s0-t/cx;   
  if(*alpha< 0.)*alpha+=24.;
  if(*alpha>24.)*alpha-=24.;
}

ad_azl(alpha,delta,a,h)     // hd,dd
double alpha,delta,*a,*h;
{
  double hour,delt,s_delta,c_delta,s_hour,c_hour,sinh,cosh,sinA,cosA;
  hour=(s0-alpha)*cx; delt=delta*cy;
  s_delta=sin(delt);  c_delta=cos(delt);
  s_hour=sin(hour);   c_hour=cos(hour);
  sinh=sphi*s_delta+cphi*c_delta*c_hour;
  cosh=sqrt(1.-sinh*sinh);
  sinA=c_delta*s_hour/cosh;
  cosA=(sphi*c_delta*c_hour-cphi*s_delta)/cosh;
  *h=atan(sinh/cosh)/cy;
  *a=atan2(sinA,cosA);
  *a/=cy;                if(*a<0.)*a+=360.;
}

put_moon()            // 66-71
{
  int i,k;
  double phase,pl,pb,xx;
  k=indexpos(head,"A81     ",72)+10;
  getmxy(&phase,&pl,&pb,&xx);
  sprintf(&head[k++],
       "MPHASE  = %20.1lf / Define the length of Moon_Month 29.53d",phase);
  sprintf(&head[k++],
       "MAZIMUTH= %20.2lf / Moon Azimuth is measured from south thr. west",pl);
  sprintf(&head[k++],"MALITIUD= %20.2lf / both units in degrees",pb);
  sprintf(&head[k++],
       "MANGLE  = %20.2lf / Position angle of Moon to CCD_Field center",xx);
  for(i=0;i<80;i++)head[k][i]=32;
}
        

moon2(jd,pl,pb)              //in degree
double jd,*pl,*pb;
{
  double t,r,x,x1,x2,x3,x4,x5,x6;
  int j;
  t=(jd-2451545.0)/36525.0/360.;
  x=481267.833*t; j=x; x-=j;
  r=218.32+x*360.;
  x=477198.85*t;  j=x; x-=j; x1=(134.9+x*360.)*cy;
  x=413335.38*t;  j=x; x-=j; x2=(259.2-x*360.)*cy;   //-
  x=890534.23*t;  j=x; x-=j; x3=(235.7+x*360.)*cy;
  x=954397.70*t;  j=x; x-=j; x4=(269.9+x*360.)*cy;
  x= 35999.05*t;  j=x; x-=j; x5=(357.5+x*360.)*cy;
  x=966404.05*t;  j=x; x-=j; x6=(186.6+x*360.)*cy;
  r+=6.29*sin(x1)-1.27*sin(x2)+0.66*sin(x3)
    +0.21*sin(x4)-0.19*sin(x5)-0.11*sin(x6);
  r/=360.; j=r; r-=j;
  *pl=r*360.;
  x=483202.03*t;  j=x; x-=j; x1=( 93.3+x*360.)*cy;
  x=960400.87*t;  j=x; x-=j; x2=(228.2+x*360.)*cy;
  x=  6003.18*t;  j=x; x-=j; x3=(318.3+x*360.)*cy;
  x=407332.20*t;  j=x; x-=j; x4=(217.6-x*360.)*cy;    //-
  r=5.13*sin(x1)+0.28*sin(x2)-0.28*sin(x3)-0.17*sin(x4);
  r/=360.; j=r; r-=j;
  r*=360.; if(r>90.)r-=360.;
  *pb=r;
}

lb_ad(pl,pb,alpha,delta)       // dd-->hd
double pl,pb,*alpha,*delta;
{
  double xl,xm,xn;
  xl=cos(pb*cy)*cos(pl*cy);
  xm=0.9175*cos(pb*cy)*sin(pl*cy)-0.3978*sin(pb*cy);
  xn=0.3978*cos(pb*cy)*sin(pl*cy)+0.9175*sin(pb*cy);
  *alpha=atan2(xm,xl)/cy/15.;
  if(*alpha<0.)*alpha+=24.;
  *delta=asin(xn)/cy;
}

ad_lb2(al,de,pl,pb)
double al,de,*pl,*pb;
{
  double xl,xm,xn;
  xl=cos(de*cy)*cos(al*cx);
  xm=0.9175*cos(de*cy)*sin(al*cx)+0.3978*sin(de*cy);
  xn=-.3978*cos(de*cy)*sin(al*cx)+0.9175*sin(de*cy);
  *pl=atan2(xm,xl)/cy;  if(*pl<0.)*pl+=360.;
  *pb=asin(xn)/cy;
}

put_ds9()                        
{
  int k=48;
  double x,z0,z1,z2,z3;
  astprs2(a_c,d_c,2000.0,&xa,&xd,epoch);
  x=(a8[1]+a8[2])/(a8[0]-a8[3])/cy;       // BOK
  z0=n1/2.;  
  z1=n2/2.; 
  sprintf(&head[k++],"CTYPE1  = 'RA---ARC'           / DS9 parameters");
  sprintf(&head[k++],"CTYPE2  = 'DEC--ARC'           /");
  sprintf(&head[k++],"CROTA2  =  %16.8le    /",x);
  sprintf(&head[k++],"CRVAL1  =  %16.8le    /",xa*15.);
  sprintf(&head[k++],"CRVAL2  =  %16.8le    /",xd);
  sprintf(&head[k++],"CDELT1  =  %16.8le    /",a8[0]/cy);
  sprintf(&head[k++],"CDELT2  =  %16.8le    /",a8[3]/cy);
  xytoad(a8,b8);      ad_xy(a_c*cx,d_c*cy,&xa,&xd);
  sprintf(&head[k++],"CRPIX1  =  %16.8le    /",xa);
  sprintf(&head[k++],"CRPIX2  =  %16.8le    /",xd);
}

//c this subroutine written by burstein in FORTRAN.
double redden(al,b)
double al,b;
{
  int ired[1200],ihi[201],ilin,kz,kq,iebv,kw,kx,ky;
  float xrad,alin,bb,arec,alc,bmv,ahi,abv;
  FILE *fp10;
  int  ip11=451200;
  int  ip12=451200;
  int  ip13=168840;

  xrad=57.2958;
  alin=al/0.3+0.5;
  ilin=alin+0.51;
  bb=b; if(bb<0.)bb=-bb;
  if(bb<10.)return(-0.99);
  fp10=fopen("/UCAC3/burstein.red","rb");
  if(fp10==0){ printf("burstein.red not found !\n"); exit(0); }
  if(b>0.)goto l18;
  if(b+62.<0.)goto l13;
  arec=-(b+10.)/0.6+1.;
  kz=arec+0.51;
  fseek(fp10,(kz-1)*4800,0); fread(ired,4,1200,fp10);   
  kq=kz-1;
  iebv=ired[ilin-1];
  goto l19;
l13:
  alc=al/xrad;
  arec=101.+sin(alc)*(90.+b)/0.3;
  alin=101.+cos(alc)*(90.+b)/0.3;
  kw=arec+0.51;
  ilin=alin+0.51;
  fseek(fp10,(kw-1)*804+ip11+ip12+ip13,0); fread(ihi,4,201,fp10);
  ahi=ihi[ilin-1];
  bmv=-0.0372 + 0.357*ahi/10000.;
  iebv=bmv*1000.+0.5;
  kq=kw-1;
  goto l19;
l18:
  if(b-62.>0)goto l16;
  arec=(b-10.)/0.6+1.;
  ky=arec+0.51;
  fseek(fp10,(ky-1)*4800+ip11,0); fread(ired,4,1200,fp10);
  kq=ky-1;
  iebv=ired[ilin-1];
  goto l19;
l16:
  alc=al/xrad;
  arec = 101. + sin(alc)*(90.-b)/0.3;
  alin = 101. + cos(alc)*(90.-b)/0.3;
  kx=arec+0.51;
  ilin=alin+0.51;
  fseek(fp10,(kx-1)*804+ip11+ip12,0); fread(ihi,4,201,fp10);
  ahi=ihi[ilin-1];
  bmv = -0.0372 + 0.357*ahi/10000.;
  iebv=bmv*1000.+0.5;
  kq=kx-1;
l19:
  fclose(fp10);
  abv=iebv;
  return(abv*4./1000. + 0.005);
}

put_a8()      // 49-58
{
  int i,j,k;
  j=k=32;
  toms2(a_c,c1,1); toms2(d_c,c2,0);
  sprintf(&head[k++],"SKYADU  =         %12.2lf / (per pixel)",sky);
  sprintf(&head[k++],"SEEING  =         %12.2lf / (arcsec)",seeing);
  i=indexpos(head,"VOLT1   ",72);
  if(i!=72){
    sprintf(&head[i+1],"VOLT2   =       %14.2lf / SKY ADU/PIXEL         ",sky);
    sprintf(&head[i+2],"VOLT3   =       %14.2lf / SEEING  (arcsec)",seeing);
  }
  sprintf(&head[k++],"A81     = %20.12le / 6 plate coef.",a8[0]);
  sprintf(&head[k++],"A82     = %20.12le / matched star:%6d",a8[1],match);
  sprintf(&head[k++],"A83     = %20.12le / RMS (arcsec):%6.2lf",a8[2],rms);
  sprintf(&head[k++],"A84     = %20.12le / total ppmxl: %6d",a8[3],nu3);
  sprintf(&head[k++],"A85     = %20.12le /",a8[4]);
  sprintf(&head[k++],"A86     = %20.12le /",a8[5]);
  sprintf(&head[k++],"A87     = %20.12le / p_c (2000.0) %s",a8[6],c1);
  sprintf(&head[k++],"A88     = %20.12le / in rad.      %s",a8[7],c2);
  ad_lb2(a_c,d_c,&xa,&xd);
  sprintf(&head[j++][56],"l_sun = %9.2lf",xa);
  sprintf(&head[j++][56],"b_sun = %9.2lf",xd);
  astprs2(a_c,d_c,2000.,&xa,&xd,1950.);
  galactic(xa*cx,xd*cy,&xi,&xn);
  sprintf(&head[j++][56],"l_gal = %9.2lf",xi);
  sprintf(&head[j++][56],"b_gal = %9.2lf",xn);
  xa=redden(xi,xn);
  sprintf(&head[k++],"EXTINCT = %20.3lf / BH Bmag, EXTINCTION = 4*E(B-V)",xa);
  astprs2(a_c,d_c,2000.,&xa,&xd,epoch);
  ad_azl(xa,xd,&xi,&xn);                     // hd,dd, 
  sprintf(&head[j++][56],"Azimuth = %7.2lf",xi);
  sprintf(&head[j++][56],"Altitude= %7.2lf",xn);
}

double xmass(angle,delta)
double angle,delta;
{
  double x,z;
  x=sin(delta)*sphi+cos(delta)*cos(angle)*cphi;
//       z=90.-atan2(x,sqrt(1.-x*x))*180./pi      ! Zenth
  x=1./x-1.;
  z=1.+x-x*(0.0018167+x*(0.002875+x*0.0008083));
  return(z);
}

put_xmass()                                    // 15,16,17,64
{
  int i,j,k;
  double xk,x,xmm,xmb,xme;
  j=(iy+((im-1)*30.4+id)/365.)*10.+0.5;
  epoch=j*0.1;
  astprs2(a_c,d_c,2000.,&xa,&xd,epoch);
  toms2(xa,c1,1);
  toms2(xd,c2,0);
  k=indexpos(head,"RA      ",72);
  sprintf(&head[k][10],"'     %s'",c1);
  sprintf(&head[k][49],"(%6.1lf)",epoch);
  k=indexpos(head,"EPOCH   ",72);
  if(k!=72)sprintf(&head[k][23]," %6.1lf",epoch); head[k][10]=32;
  k=indexpos(head,"DEC     ",72);
  sprintf(&head[k][10],"'      %s'",c2);
  sprintf(&head[k][49],"(%6.1lf)",epoch);
  ha=s0-xa; while(ha>=24.)ha-=24.;
            while(ha<0.)  ha+=24.;
  toms2(ha,c1,1);
  k=indexpos(head,"HA      ",72);
  sprintf(&head[k][10],"'     %s'",c1);
  xmm=xmass(ha*cx,xd*cy);                            // middle time
  k=indexpos(head,"EXPTIME ",72);
  if(k==72)k=indexpos(head,"EXPOSURE",72);
  sscanf(&head[k][19],"%d",&xk);                     // exp time
  xmb=xmass((ha-xk/7200.)*cx,xd*cy);                 // begin time
  xme=xmass((ha+xk/7200.)*cx,xd*cy);                 // end time
  x=(xmb+4.*xmm+xme)/6.;
  k=indexpos(head,"A81     ",72)+9;
  sprintf(&head[k],"AIRMASS = %20.3lf /",x);
}
galactic(alpha,delta,xl,xb)
double alpha,delta,*xl,*xb;
{
  double a0,theta,xl0,x1,x2;
  a0=282.25*cy;
  theta=62.6*cy;
  xl0=33.0*cy;
  x1=cos(delta)*cos(alpha-a0);
  x2=cos(delta)*sin(alpha-a0)*cos(theta)+sin(delta)*sin(theta);
  *xl=atan2(x2,x1)+xl0;
  if(*xl<0.)*xl+=2.0*pi;
  *xb=asin(sin(delta)*cos(theta)-cos(delta)*sin(alpha-a0)*sin(theta));
  *xl/=cy;  *xb/=cy;
}

/*
        subroutine newhead(head,nh)
        character*80 head(1),ohead(108)
        do 10 i=1,nh
10      ohead(i)=head(i)
        do 20 i=6,72
        do 20 j=1,80
20      head(i)(j:j)=char(32)
c        call puthead( 6,head,ohead,"CCDNAME ",nh)
        call puthead( 7,head,ohead,"OBJECT  ",nh)
        call puthead( 8,head,ohead,"OBSERVER",nh)
        call puthead( 9,head,ohead,"INSTRUME",nh)
        call puthead(10,head,ohead,"TELESCOP",nh)
        call puthead(11,head,ohead,"SITELONG",nh)
        call puthead(12,head,ohead,"SITELAT ",nh)
        call puthead(13,head,ohead,"SITEELEV",nh)
        call puthead(14,head,ohead,"FILTER  ",nh)
        call puthead(15,head,ohead,"GAIN8   ",nh)
        call puthead(16,head,ohead,"RDNOISE8",nh)
        call puthead(17,head,ohead,"DARKCUR ",nh)
        call puthead(18,head,ohead,"CAMTEMP ",nh)
        call puthead(19,head,ohead,"DEWTEMP ",nh)
        call puthead(20,head,ohead,"DATE-OBS",nh)
        call puthead(21,head,ohead,"TIME-OBS",nh)
        call puthead(22,head,ohead,"EXPTIME ",nh)
        call puthead(23,head,ohead,"RA      ",nh)
        call puthead(24,head,ohead,"DEC     ",nh)
        call puthead(25,head,ohead,"EPOCH   ",nh)
        call puthead(26,head,ohead,"HA      ",nh)
        head(11)(12:23)="111:36:01.6'"
        head(12)(12:23)=" 31:57:46.5'"
        head(13)(12:23)="       2071'"
        call puthead(72,head,ohead,"END     ",nh)
        end

        subroutine puthead(k,head,ohead,c8,nh)
        character*80 head(1),ohead(1),c8*8
        ipos=indexpos(ohead,c8,nh)
        head(k)=ohead(ipos)
        end
*/

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;
/*                         // same discrube formula as following 3 lines
  double sd,cd,td,co;      // china bai_ke_quan_shu (astronmy) p.552
  sd=sin(dc);  cd=cos(dc);  td=tan(der);
  co=cos(rar-ac);  tmp=sd*td+cd*co;
  xi=sin(rar-ac)/tmp;
  xn=(cd*td-sd*co)/tmp;
*/
  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);
    tmp=atan(tmp)/tmp;
    xi*=tmp; xn*=tmp;
  }
  if(type==2){
    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;}
}

xy_rad(ac,dc,x,y,ra,de,type,k)
double ac,dc,x,y,*ra,*de;  int type,k;
{
  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];
  if(type==1){
    tmp=sqrt(xi*xi+xn*xn);
    tmp=tan(tmp)/tmp;
    xi*=tmp;  xn*=tmp;
  }
  if(type==2){
    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);
//                                 printf("%lf\n",tmp);
    }
    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);
  if(k==0){ *ra/=cx; *de/=cy; }
}