view release on metacpan or  search on metacpan

palsrc/palAopqk.c  view on Meta::CPAN

  /*  Breakpoint for fast/slow refraction algorithm:
   *  ZD greater than arctan(4), (see palRefco routine)
   *  or vector Z less than cosine(arctan(Z)) = 1/sqrt(17) */
  const double zbreak = 0.242535625;
  int i;

  double  sphi,cphi,st,v[3],xhd,yhd,zhd,diurab,f,
    xhdt,yhdt,zhdt,xaet,yaet,zaet,azobs,
    zdt,refa,refb,zdobs,dzd,dref,ce,
    xaeo,yaeo,zaeo,hmobs,dcobs,raobs;

  /*  sin, cos of latitude */
  sphi = aoprms[1];
  cphi = aoprms[2];

  /*  local apparent sidereal time */
  st = aoprms[13];

  /*  apparent ra,dec to cartesian -ha,dec */
  palDcs2c( rap-st, dap, v );
  xhd = v[0];
  yhd = v[1];
  zhd = v[2];

  /*  diurnal aberration */
  diurab = aoprms[3];
  f = (1.0-diurab*yhd);
  xhdt = f*xhd;
  yhdt = f*(yhd+diurab);
  zhdt = f*zhd;

  /*  cartesian -ha,dec to cartesian az,el (s=0,e=90) */
  xaet = sphi*xhdt-cphi*zhdt;
  yaet = yhdt;
  zaet = cphi*xhdt+sphi*zhdt;

  /*  azimuth (n=0,e=90) */
  if (xaet == 0.0 && yaet == 0.0) {
    azobs = 0.0;
  } else {
    azobs = atan2(yaet,-xaet);
  }

  /*  topocentric zenith distance */
  zdt = atan2(sqrt(xaet*xaet+yaet*yaet),zaet);

  /*
   *  refraction
   *  ---------- */

  /*  fast algorithm using two constant model */
  refa = aoprms[10];
  refb = aoprms[11];
  palRefz(zdt,refa,refb,&zdobs);

  /*  large zenith distance? */
  if (cos(zdobs) < zbreak) {

    /*     yes: use rigorous algorithm */

    /*     initialize loop (maximum of 10 iterations) */
    i = 1;
    dzd = 1.0e1;
    while (fabs(dzd) > 1e-10 && i <= 10) {

      /*        compute refraction using current estimate of observed zd */
      palRefro(zdobs,aoprms[4],aoprms[5],aoprms[6],
               aoprms[7],aoprms[8],aoprms[0],
               aoprms[9],1e-8,&dref);

      /*        remaining discrepancy */
      dzd = zdobs+dref-zdt;

      /*        update the estimate */
      zdobs = zdobs-dzd;

      /*        increment the iteration counter */
      i++;
    }
  }

  /*  to cartesian az/zd */
  ce = sin(zdobs);
  xaeo = -cos(azobs)*ce;
  yaeo = sin(azobs)*ce;
  zaeo = cos(zdobs);

  /*  cartesian az/zd to cartesian -ha,dec */
  v[0] = sphi*xaeo+cphi*zaeo;
  v[1] = yaeo;
  v[2] = -cphi*xaeo+sphi*zaeo;

  /*  to spherical -ha,dec */
  palDcc2s(v,&hmobs,&dcobs);

  /*  right ascension */
  raobs = palDranrm(st+hmobs);

  /*  return the results */
  *aob = azobs;
  *zob = zdobs;
  *hob = -hmobs;
  *dob = dcobs;
  *rob = raobs;

}