view release on metacpan or  search on metacpan

palsrc/palAopqk.c  view on Meta::CPAN

*       programming in the calling application may allow the
*       problem to be reduced.  Prepare an alternative AOPRMS array,
*       computed for zero air-pressure;  this will disable the
*       refraction corrections and cause rapid execution.  Using
*       this AOPRMS array, a preliminary call to the present routine
*       will, depending on the application, produce a rough position
*       which may be enough to establish whether the full, slow
*       calculation (using the real AOPRMS array) is worthwhile.
*       For example, there would be no need for the full calculation
*       if the preliminary call had already established that the
*       source was well below the elevation limits for a particular
*       telescope.
*
*     - The azimuths etc produced by the present routine are with
*       respect to the celestial pole.  Corrections to the terrestrial
*       pole can be computed using palPolmo.

*  History:
*     2012-08-25 (TIMJ):
*        Initial version, copied from Fortran SLA
*        Adapted with permission from the Fortran SLALIB library.
*     {enter_further_changes_here}

*  Copyright:
*     Copyright (C) 2003 Rutherford Appleton Laboratory
*     Copyright (C) 2012 Science and Technology Facilities Council.
*     All Rights Reserved.

*  Licence:
*     This program is free software; you can redistribute it and/or
*     modify it under the terms of the GNU General Public License as
*     published by the Free Software Foundation; either version 3 of
*     the License, or (at your option) any later version.
*
*     This program is distributed in the hope that it will be
*     useful, but WITHOUT ANY WARRANTY; without even the implied
*     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
*     PURPOSE. See the GNU General Public License for more details.
*
*     You should have received a copy of the GNU General Public License
*     along with this program; if not, write to the Free Software
*     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
*     MA 02110-1301, USA.

*  Bugs:
*     {note_any_bugs_here}
*-
*/

#include <math.h>

#include "pal.h"

void palAopqk ( double rap, double dap, const double aoprms[14],
                double *aob, double *zob, double *hob,
                double *dob, double *rob ) {

  /*  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) {