home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part21 / refract.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  1.8 KB  |  67 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* correct the true altitude, ta, for refraction to the apparent altitude, aa,
  6.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  7.  * the temperature, tr, in degrees C.
  8.  */
  9. void
  10. refract (pr, tr, ta, aa)
  11. double pr, tr;
  12. double ta;
  13. double *aa;
  14. {
  15.     double r;    /* refraction correction*/
  16.  
  17.         if (ta >= degrad(15.)) {
  18.         /* model for altitudes at least 15 degrees above horizon */
  19.             r = 7.888888e-5*pr/((273+tr)*tan(ta));
  20.     } else if (ta > degrad(-5.)) {
  21.         /* hairier model for altitudes at least -5 and below 15 degrees */
  22.         double a, b, tadeg = raddeg(ta);
  23.         a = ((2e-5*tadeg+1.96e-2)*tadeg+1.594e-1)*pr;
  24.         b = (273+tr)*((8.45e-2*tadeg+5.05e-1)*tadeg+1);
  25.         r = degrad(a/b);
  26.     } else {
  27.         /* do nothing if more than 5 degrees below horizon.
  28.          */
  29.         r = 0;
  30.     }
  31.  
  32.     *aa  =  ta + r;
  33. }
  34.  
  35. /* correct the apparent altitude, aa, for refraction to the true altitude, ta,
  36.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  37.  * the temperature, tr, in degrees C.
  38.  */
  39. void
  40. unrefract (pr, tr, aa, ta)
  41. double pr, tr;
  42. double aa;
  43. double *ta;
  44. {
  45.     double err;
  46.     double appar;
  47.     double true;
  48.  
  49.     /* iterative solution: search for the true that refracts to the
  50.      *   given apparent.
  51.      * since refract() is discontinuous at -5 degrees, there is a range
  52.      *   of apparent altitudes between about -4.5 and -5 degrees that are
  53.      *   not invertable (the graph of ap vs. true has a vertical step at
  54.      *   true = -5). thus, the iteration just oscillates if it gets into
  55.      *   this region. if this happens the iteration is forced to abort.
  56.      *   of course, this makes unrefract() discontinuous too.
  57.      */
  58.     true = aa;
  59.     do {
  60.         refract (pr, tr, true, &appar);
  61.         err = appar - aa;
  62.         true -= err;
  63.     } while (fabs(err) >= 1e-6 && true > degrad(-5));
  64.  
  65.     *ta = true;
  66. }
  67.