-
Notifications
You must be signed in to change notification settings - Fork 12
/
mathutil.c
136 lines (124 loc) · 3.41 KB
/
mathutil.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
/*=======================================================================
solunar
mathutil.c
Functions for finding maxima, minima, axis crossings, etc.
(c)2005-2012 Kevin Boone
=======================================================================*/
#include <math.h>
#include "mathutil.h"
void mathutil_get_maxima (double *x, double *y, int npoints, double *mins,
int maxmins, int *nmins)
{
*nmins = 0;
int i;
for (i = 1; i < npoints - 1 && *nmins < maxmins; i++)
{
double y0 = y[i - 1];
double y1 = y[i];
double y2 = y[i + 1];
if (y0 < y1 && y2 < y1)
{
double A = (0.5 * (y0 + y2)) - y1;
double B = (0.5 * (y2 - y0));
double xExtreme = -B / (2.0 * A);
double xguess = x[i-1] + ((1 + xExtreme) * (x[i] - x[i-1]));
mins[*nmins] = xguess;
(*nmins)++;
}
}
}
void mathutil_get_minima (double *x, double *y, int npoints, double *mins,
int maxmins, int *nmins)
{
*nmins = 0;
int i;
for (i = 1; i < npoints - 1 && *nmins < maxmins; i++)
{
double y0 = y[i - 1];
double y1 = y[i];
double y2 = y[i + 1];
if (y0 > y1 && y2 > y1)
{
double A = (0.5 * (y0 + y2)) - y1;
double B = (0.5 * (y2 - y0));
double xExtreme = -B / (2.0 * A);
double xguess = x[i-1] + ((1 + xExtreme) * (x[i] - x[i-1]));
mins[*nmins] = xguess;
(*nmins)++;
}
}
}
// Find zero-crossings using muller's method
void mathutil_get_positive_axis_crossings (double *x, double *y,
int npoints, double *mins, int maxmins, int *nmins)
{
*nmins = 0;
int i;
for (i = 1; i < npoints - 1 && *nmins < maxmins; i++)
{
double y0 = y[i - 1];
double y1 = y[i];
double y2 = y[i + 1];
if (y0 < 0 && y1 <= 0 && y2 > 0)
{
double root1 = 0.0;
double root2 = 0.0;
double A = (0.5 * (y0 + y2)) - y1;
double B = (0.5 * (y2 - y0));
double C = y1;
double xExtreme = -B / (2.0 * A);
double discriminant = (B * B) - 4.0 * A * C;
if (discriminant >= 0.0)
{
double DX = 0.5 * sqrt(discriminant) / fabs(A);
root1 = xExtreme - DX;
root2 = xExtreme + DX;
double bisect = 0;
if (fabs (root2) < 1.0)
bisect = root2;
else if (fabs (root1) < 1.0)
bisect = root1;
double xguess = x[i] + (bisect * (x[i] - x[i-1]));
mins[*nmins] = xguess;
(*nmins)++;
}
}
}
}
// Find zero-crossings using muller's method
void mathutil_get_negative_axis_crossings (double *x, double *y,
int npoints, double *mins, int maxmins, int *nmins)
{
*nmins = 0;
int i;
for (i = 1; i < npoints - 1 && *nmins < maxmins; i++)
{
double y0 = y[i - 1];
double y1 = y[i];
double y2 = y[i + 1];
if (y0 > 0 && y1 >= 0 && y2 < 0)
{
double root1 = 0.0;
double root2 = 0.0;
double A = (0.5 * (y0 + y2)) - y1;
double B = (0.5 * (y2 - y0));
double C = y1;
double xExtreme = -B / (2.0 * A);
double discriminant = (B * B) - 4.0 * A * C;
if (discriminant >= 0.0)
{
double DX = 0.5 * sqrt(discriminant) / fabs(A);
root1 = xExtreme - DX;
root2 = xExtreme + DX;
double bisect = 0;
if (fabs (root2) < 1.0)
bisect = root2;
else if (fabs (root1) < 1.0)
bisect = root1;
double xguess = x[i] + (bisect * (x[i] - x[i-1]));
mins[*nmins] = xguess;
(*nmins)++;
}
}
}
}