1 /* This file is part of the GNU plotutils package. */
2
3 /*
4 * Copyright (C) 1982-1994, Nicholas B. Tufillaro. All rights reserved.
5 *
6 * GNU enhancements Copyright (C) 1996, 1997, 1998, 1999, 2005, 2008, Free
7 * Software Foundation, Inc.
8 */
9
10 /*
11 * find maximum errors
12 *
13 */
14
15 #include "sys-defines.h"
16 #include "ode.h"
17 #include "extern.h"
18
19 static double ssemax, abemax, acemax;
20 static char *ssenam, *abenam, *acenam;
21
22 void
maxerr(void)23 maxerr (void)
24 {
25 struct sym *sp, *dq;
26
27 dq = symtab->sy_link;
28 ssemax = abemax = acemax = 0.0;
29 for (sp = dq; sp != NULL; sp = sp->sy_link)
30 {
31 if (ssemax < sp->sy_sserr)
32 {
33 ssemax = sp->sy_sserr;
34 ssenam = sp->sy_name;
35 }
36 if (abemax < sp->sy_aberr)
37 {
38 abemax = sp->sy_aberr;
39 abenam = sp->sy_name;
40 }
41 if (acmax < sp->sy_acerr)
42 {
43 acemax = sp->sy_acerr;
44 acenam = sp->sy_name;
45 }
46 }
47 }
48
49 bool
hierror(void)50 hierror (void) /* not enough accuracy */
51 {
52 double t = symtab->sy_val[0];
53
54 if (t + tstep == t)
55 {
56 fprintf (stderr, "%s: %s\n", progname, "step size below lower limit");
57 longjmp (mark, 1);
58 }
59 if (ssemax <= ssmax && abemax <= abmax && acemax <= acmax)
60 return false;
61 if (fabs(tstep) >= fabs(hmin))
62 return true;
63 if (sflag)
64 return false;
65 if (ssemax > ssmax)
66 fprintf (stderr,
67 "%s: relative error limit exceeded while calculating %.*s'\n",
68 progname, NAMMAX, ssenam);
69 else if (abemax > abmax)
70 fprintf (stderr,
71 "%s: absolute error limit exceeded while calculating %.*s'\n",
72 progname, NAMMAX, abenam);
73 else if (acemax > acmax)
74 fprintf (stderr,
75 "%s: accumulated error limit exceeded while calculating %.*s'\n",
76 progname, NAMMAX, acenam);
77 longjmp (mark, 1);
78
79 /* doesn't return, but must keep unintelligent compilers happy */
80 return false;
81 }
82
83 bool
lowerror(void)84 lowerror (void) /* more than enough accuracy */
85 {
86 if (ssemax < ssmin || abemax < abmin)
87 if (fabs(tstep) <= fabs(hmax))
88 return true;
89 return false;
90 }
91
92 /*
93 * interpolate to tstop in Runge-Kutta routines
94 */
95 #define PASTSTOP(stepvar) (t + 0.9375*stepvar > tstop && \
96 t + 0.0625*stepvar < tstop)
97 #define BEFORESTOP(stepvar) (t + 0.9375*stepvar < tstop && \
98 t + 0.0625*stepvar > tstop)
99
100 bool
intpr(double t)101 intpr (double t)
102 {
103 if (tstep > 0)
104 if (!PASTSTOP(tstep))
105 return false;
106 if (tstep < 0)
107 if (!BEFORESTOP(tstep))
108 return false;
109 if (tstep > 0)
110 while (PASTSTOP(tstep))
111 tstep = HALF * tstep;
112 if (tstep < 0)
113 while (BEFORESTOP(tstep))
114 tstep = HALF * tstep;
115 return true;
116 }
117