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