home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / ephem421.zip / SRCH.C < prev    next >
C/C++ Source or Header  |  1990-09-13  |  8KB  |  348 lines
  1. /* this file contains functions to support iterative ephem searches.
  2.  * we support several kinds of searching and solving algorithms.
  3.  * values used in the evaluations come from the field logging flog.c system.
  4.  * the expressions being evaluated are compiled and executed from compiler.c.
  5.  */
  6.  
  7. #include <stdio.h>
  8. #include <math.h>
  9. #include "screen.h"
  10.  
  11. extern char *strcpy();
  12.  
  13. static int (*srch_f)();
  14. static int srch_tmscalled;
  15. static char expbuf[NC];        /* [0] == '\0' when expression is invalid */
  16. static double tmlimit = 1./60.;    /* search accuracy, in hrs; def is one minute */
  17.  
  18.  
  19. srch_setup()
  20. {
  21.     int srch_minmax(), srch_solve0(), srch_binary();
  22.     static char *chcs[] = {
  23.         "Find extreme", "Find 0", "Binary", "New function", "Accuracy",
  24.         "Stop"
  25.     };
  26.     static int fn;    /* start with 0, then remember for next time */
  27.  
  28.     /* let op select algorithm, edit, set accuracy
  29.      * or stop if currently searching
  30.      * algorithms require a function.
  31.      */
  32.     ask:
  33.     switch (popup(chcs, fn, srch_f ? 6 : 5)) {
  34.     case 0: fn = 0;
  35.         if (expbuf[0] == '\0')
  36.             set_function();
  37.         srch_f = expbuf[0] ? srch_minmax : 0;
  38.         if (srch_f)
  39.             break;
  40.         else
  41.             goto ask;
  42.     case 1: fn = 1;
  43.         if (expbuf[0] == '\0')
  44.             set_function();
  45.         srch_f = expbuf[0] ? srch_solve0 : 0;
  46.         if (srch_f)
  47.             break;
  48.         else
  49.             goto ask;
  50.     case 2: fn = 2;
  51.         if (expbuf[0] == '\0')
  52.             set_function();
  53.         srch_f = expbuf[0] ? srch_binary : 0;
  54.         if (srch_f)
  55.             break;
  56.         else
  57.             goto ask;
  58.     case 3: fn = 3; srch_f = 0; set_function(); goto ask;
  59.     case 4: fn = 4; srch_f = 0; set_accuracy(); goto ask;
  60.     case 5: srch_f = 0; srch_prstate(0); return;
  61.     default: return;
  62.     }
  63.  
  64.     /* new search */
  65.     srch_tmscalled = 0;
  66.     srch_prstate (0);
  67. }
  68.  
  69. /* if searching is in effect call the search type function.
  70.  * it might modify *tmincp according to where it next wants to eval.
  71.  * (remember tminc is in hours, not days).
  72.  * if searching ends for any reason it is also turned off.
  73.  * also, flog the new value.
  74.  * return 0 if caller can continue or -1 if it is time to stop.
  75.  */
  76. srch_eval(mjd, tmincp)
  77. double mjd;
  78. double *tmincp;
  79. {
  80.     char errbuf[128];
  81.     int s;
  82.     double v;
  83.  
  84.     if (!srch_f)
  85.         return (0);
  86.  
  87.     if (execute_expr (&v, errbuf) < 0) {
  88.         srch_f = 0;
  89.         f_msg (errbuf);
  90.     } else {
  91.         s = (*srch_f)(mjd, v, tmincp);
  92.         if (s < 0)
  93.         srch_f = 0;
  94.         (void) flog_log (R_SRCH, C_SRCH, v, "");
  95.         srch_tmscalled++;
  96.     }
  97.  
  98.     srch_prstate (0);
  99.     return (s);
  100. }
  101.  
  102. /* print state of searching. */
  103. srch_prstate (force)
  104. int force;
  105. {
  106.     int srch_minmax(), srch_solve0(), srch_binary();
  107.     static (*last)();
  108.  
  109.     if (force || srch_f != last) {
  110.         f_string (R_SRCH, C_SRCHV,
  111.             srch_f == srch_minmax   ? "Extrema" :
  112.             srch_f == srch_solve0   ? " Find 0" :
  113.             srch_f == srch_binary ?   " Binary" :
  114.                           "    off");
  115.         last = srch_f;
  116.     }
  117. }
  118.  
  119. srch_ison()
  120. {
  121.     return (srch_f != 0);
  122. }
  123.  
  124. /* display current expression. then if type in at least one char make it the
  125.  * current expression IF it compiles ok.
  126.  * TODO: editing?
  127.  */
  128. static
  129. set_function()
  130. {
  131.     static char prompt[] = "Function: ";
  132.     char newexp[NC];
  133.     int s;
  134.  
  135.     f_prompt (prompt);
  136.     (void) fputs (expbuf, stdout);
  137.     c_pos (R_PROMPT, sizeof(prompt));
  138.  
  139.     s = read_line (newexp, PW-sizeof(prompt));
  140.     if (s >= 0) {
  141.         char errbuf[NC];
  142.         if (s > 0 && compile_expr (newexp, errbuf) < 0)
  143.         f_msg (errbuf);
  144.         else
  145.         (void) strcpy (expbuf, newexp);
  146.     }
  147. }
  148.  
  149. static
  150. set_accuracy()
  151. {
  152.     static char p[] = "Desired accuracy (         hrs): ";
  153.     int hrs, mins, secs;
  154.     char buf[NC];
  155.  
  156.     f_prompt (p);
  157.     f_time (R_PROMPT, C_PROMPT+18, tmlimit); /* place in blank spot */
  158.     c_pos (R_PROMPT, sizeof(p));
  159.     if (read_line (buf, PW-sizeof(p)) > 0) {
  160.         f_dec_sexsign (tmlimit, &hrs, &mins, &secs);
  161.         f_sscansex (buf, &hrs, &mins, &secs);
  162.         sex_dec (hrs, mins, secs, &tmlimit);
  163.     }
  164. }
  165.  
  166. /* use successive paraboloidal fits to find when expression is at a
  167.  * local minimum or maximum.
  168.  */
  169. static
  170. srch_minmax(mjd, v, tmincp)
  171. double mjd;
  172. double v;
  173. double *tmincp;
  174. {
  175.     static double base;        /* for better stability */
  176.     static double x_1, x_2, x_3;    /* keep in increasing order */
  177.     static double y_1, y_2, y_3;
  178.     double xm, a, b;
  179.  
  180.     if (srch_tmscalled == 0) {
  181.         base = mjd;
  182.         x_1 = 0.0;
  183.         y_1 = v;
  184.         return (0);
  185.     }
  186.     mjd -= base;
  187.     if (srch_tmscalled == 1) {
  188.         /* put in one of first two slots */
  189.         if (mjd < x_1) {
  190.             x_2 = x_1;  y_2 = y_1;
  191.         x_1 = mjd; y_1 = v;
  192.         } else {
  193.         x_2 = mjd; y_2 = v;
  194.         }
  195.         return (0);
  196.     }
  197.     if (srch_tmscalled == 2 || fabs(mjd - x_1) < fabs(mjd - x_3)) {
  198.         /* closer to x_1 so discard x_3.
  199.          * or if it's our third value we know to "discard" x_3.
  200.          */
  201.         if (mjd > x_2) {
  202.         x_3 = mjd; y_3 = v;
  203.         } else {
  204.         x_3 = x_2;  y_3 = y_2;
  205.         if (mjd > x_1) {
  206.             x_2 = mjd; y_2 = v;
  207.         } else {
  208.             x_2 = x_1;  y_2 = y_1;
  209.             x_1 = mjd; y_1 = v;
  210.         }
  211.         }
  212.         if (srch_tmscalled == 2)
  213.         return (0);
  214.     } else {
  215.         /* closer to x_3 so discard x_1 */
  216.         if (mjd < x_2) {
  217.         x_1 = mjd;  y_1 = v;
  218.         } else {
  219.         x_1 =  x_2;  y_1 = y_2;
  220.         if (mjd < x_3) {
  221.             x_2 = mjd; y_2 = v;
  222.         } else {
  223.             x_2 =  x_3; y_2 = y_3;
  224.             x_3 = mjd; y_3 = v;
  225.         }
  226.         }
  227.     }
  228.  
  229. #ifdef TRACEMM
  230.     { char buf[NC];
  231.       sprintf (buf, "x_1=%g y_1=%g x_2=%g y_2=%g x_3=%g y_3=%g",
  232.                         x_1, y_1, x_2, y_2, x_3, y_3);
  233.       f_msg (buf);
  234.     }
  235. #endif
  236.     a = y_1*(x_2-x_3) - y_2*(x_1-x_3) + y_3*(x_1-x_2);
  237.     if (fabs(a) < 1e-10) {
  238.         /* near-0 zero denominator, ie, curve is pretty flat here,
  239.          * so assume we are done enough.
  240.          * signal this by forcing a 0 tminc.
  241.          */
  242.         *tmincp = 0.0;
  243.         return (-1);
  244.     }
  245.     b = (x_1*x_1)*(y_2-y_3) - (x_2*x_2)*(y_1-y_3) + (x_3*x_3)*(y_1-y_2);
  246.     xm = -b/(2.0*a);
  247.     *tmincp = (xm - mjd)*24.0;
  248.     return (fabs (*tmincp) < tmlimit ? -1 : 0);
  249. }
  250.  
  251. /* use secant method to solve for time when expression passes through 0.
  252.  */
  253. static
  254. srch_solve0(mjd, v, tmincp)
  255. double mjd;
  256. double v;
  257. double *tmincp;
  258. {
  259.     static double x0, x_1;    /* x(n-1) and x(n) */
  260.     static double y_0, y_1;    /* y(n-1) and y(n) */
  261.     double x_2;        /* x(n+1) */
  262.     double df;        /* y(n) - y(n-1) */
  263.  
  264.     switch (srch_tmscalled) {
  265.     case 0: x0 = mjd; y_0 = v; return(0);
  266.     case 1: x_1 = mjd; y_1 = v; break;
  267.     default: x0 = x_1; y_0 = y_1; x_1 = mjd; y_1 = v; break;
  268.     }
  269.  
  270.     df = y_1 - y_0;
  271.     if (fabs(df) < 1e-10) {
  272.         /* near-0 zero denominator, ie, curve is pretty flat here,
  273.          * so assume we are done enough.
  274.          * signal this by forcing a 0 tminc.
  275.          */
  276.         *tmincp = 0.0;
  277.         return (-1);
  278.     }
  279.     x_2 = x_1 - y_1*(x_1-x0)/df;
  280.     *tmincp = (x_2 - mjd)*24.0;
  281.     return (fabs (*tmincp) < tmlimit ? -1 : 0);
  282. }
  283.  
  284. /* binary search for time when expression changes from its initial state.
  285.  * if the change is outside the initial tminc range, then keep searching in that
  286.  *    direction by tminc first before starting to divide down.
  287.  */
  288. static
  289. srch_binary(mjd, v, tmincp)
  290. double mjd;
  291. double v;
  292. double *tmincp;
  293. {
  294.     static double lb, ub;        /* lower and upper bound */
  295.     static int initial_state;
  296.     int this_state = v >= 0.5;
  297.  
  298. #define    FLUNDEF    -9e10
  299.  
  300.     if (srch_tmscalled == 0) {
  301.         if (*tmincp >= 0.0) {
  302.         /* going forwards in time so first mjd is lb and no ub yet */
  303.         lb = mjd;
  304.         ub = FLUNDEF;
  305.         } else {
  306.         /* going backwards in time so first mjd is ub and no lb yet */
  307.         ub = mjd;
  308.         lb = FLUNDEF;
  309.         }
  310.         initial_state = this_state;
  311.         return (0);
  312.     }
  313.  
  314.     if (ub != FLUNDEF && lb != FLUNDEF) {
  315.         if (this_state == initial_state)
  316.         lb = mjd;
  317.         else
  318.         ub = mjd;
  319.         *tmincp = ((lb + ub)/2.0 - mjd)*24.0;
  320. #ifdef TRACEBIN
  321.         { char buf[NC];
  322.           sprintf (buf, "lb=%g ub=%g tminc=%g mjd=%g is=%d ts=%d",
  323.                 lb, ub, *tmincp, mjd, initial_state, this_state);
  324.           f_msg (buf);
  325.         }
  326. #endif
  327.         /* signal to stop if asking for time change less than TMLIMIT */
  328.         return (fabs (*tmincp) < tmlimit ? -1 : 0);
  329.     } else if (this_state != initial_state) {
  330.         /* gone past; turn around half way */
  331.         if (*tmincp >= 0.0)
  332.         ub = mjd;
  333.         else
  334.         lb = mjd;
  335.         *tmincp /= -2.0;
  336.         return (0);
  337.     } else {
  338.         /* just keep going, looking for first state change but we keep
  339.          * learning the lower (or upper, if going backwards) bound.
  340.          */
  341.         if (*tmincp >= 0.0)
  342.         lb = mjd;
  343.         else
  344.         ub = mjd;
  345.         return (0);
  346.     }
  347. }
  348.