56 double Interpola(
double vizq,
double vder,
double axid,
double xif) {
57 return vizq + ((vder - vizq) / axid) * xif;
60 double QuadraticEqP(
double A,
double B,
double C) {
61 return (-B + sqrt(B * B - 4 * A * C)) / 2 / A;
64 double QuadraticEqN(
double A,
double B,
double C) {
65 return (-B - sqrt(B * B - 4 * A * C)) / 2 / A;
68 float Max(
const double &a,
const float &b) {
69 return b > a ? (b) :
float(a);
72 float Max(
const float &a,
const double &b) {
73 return b > a ? float(b) : (a);
76 float Min(
const double &a,
const float &b) {
77 return b < a ? (b) :
float(a);
80 float Min(
const float &a,
const double &b) {
81 return b < a ? float(b) : (a);
84 float Sign(
const float &a,
const double &b) {
85 return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
88 float Sign(
const double &a,
const float &b) {
89 return (
float)(b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));
95 stPolar::stPolar(
double X,
double Y) {
96 Mod = sqrt(X * X + Y * Y);
100 void stPolar::operator()(
double X,
double Y) {
101 Mod = sqrt(X * X + Y * Y);
105 stRectan::stRectan() {
108 stRectan::stRectan(
double Ang,
double Mod) {
113 void stRectan::operator()(
double Ang,
double Mod) {
118 Base_interp::Base_interp() :
119 n(0), mm(0), jsav(0), cor(0), xx(), yy() {
122 Base_interp::Base_interp(
dVector &x,
const double *y,
int m) :
123 n(x.size()), mm(m), jsav(0), cor(0), xx(&x[0]), yy(y) {
124 dj = Max(1, (
int)
pow025((
double) n));
127 double Base_interp::interp(
double x) {
128 int jlo = cor ? hunt(x) : locate(x);
129 return rawinterp(jlo, x);
132 int Base_interp::locate(
const double x) {
133 int ju = 0, jm = 0, jl = 0;
134 if(n < 2 || mm < 2 || mm > n)
135 throw(
"locate size error");
136 bool ascnd = (xx[n - 1] >= xx[0]);
141 if((x >= xx[jm]) == ascnd) {
147 cor = abs(jl - jsav) > dj ? 0 : 1;
149 return Max(0, Min(n - mm, jl - ((mm - 2) >> 1)));
153 int Base_interp::hunt(
const double x) {
154 int jl = jsav, jm, ju, inc = 1;
155 if(n < 2 || mm < 2 || mm > n)
156 throw(
"locate size error");
157 bool ascnd = (xx[n - 1] >= xx[0]);
158 if(jl < 0 || jl > n - 1) {
162 if((x >= xx[jl]) == ascnd) {
168 }
else if((x < xx[ju]) == ascnd)
182 }
else if((x >= xx[jl]) == ascnd)
193 if((x >= xx[jm]) == ascnd) {
199 cor = abs(jl - jsav) > dj ? 0 : 1;
201 return Max(0, Min(n - mm, jl - ((mm - 2) >> 1)));
204 Linear_interp::Linear_interp() :
219 dj = Max(1, (
int)
pow025((
double) n));
222 double Linear_interp::rawinterp(
int j,
double x) {
223 if(xx[j] == xx[j + 1]) {
226 return yy[j] + ((x - xx[j]) / (xx[j + 1] - xx[j])) * (yy[j + 1] - yy[j]);
230 Hermite_interp::Hermite_interp() :
237 sety2(&xv[0], &yv[0]);
248 dj = Max(1, (
int)
pow025((
double) n));
249 sety2(&xv[0], &yv[0]);
252 void Hermite_interp::sety2(
const double *xv,
const double *yv) {
253 double DeltaK = 0., AlphaK = 0., BetaK = 0., TauK = 0.;
255 for(
int i = 1; i < n - 1; ++i) {
256 y2[i] = (yv[i] - yv[i - 1]) / 2. / (xv[i] - xv[i - 1]) + (yv[i + 1] - yv[i]) / 2. / (xv[i + 1] - xv[i]);
258 y2[0] = (yv[1] - yv[0]) / (xv[1] - xv[0]);
259 y2[n - 1] = (yv[n - 1] - yv[n - 2]) / (xv[n - 1] - xv[n - 2]);
261 for(
int i = 0; i < n - 1; i++) {
262 DeltaK = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]);
267 AlphaK = y2[i] / DeltaK;
268 BetaK = y2[i + 1] / DeltaK;
269 if(BetaK * BetaK + AlphaK * AlphaK > 9) {
270 TauK = 3 / sqrt(BetaK * BetaK + AlphaK * AlphaK);
271 y2[i] = TauK * AlphaK * DeltaK;
272 y2[i + 1] = TauK * BetaK * DeltaK;
278 double Hermite_interp::rawinterp(
int j,
double x) {
279 double ret_val = 0., h00 = 0., h10 = 0., h01 = 0., h11 = 0., t2 = 0., t3 = 0., t = 0., h = 0.;
284 }
else if(x >= xx[j + 1]) {
287 h = (xx[j + 1] - xx[j]);
291 h00 = 2 * t3 - 3 * t2 + 1;
292 h10 = t3 - 2 * t2 + t;
293 h01 = -2 * t3 + 3 * t2;
295 ret_val = h00 * yy[j] + h * h10 * y2[j] + h01 * yy[j + 1] + h * h11 * y2[j + 1];
300 Step_interp::Step_interp() :
315 dj = Max(1, (
int)
pow025((
double) n));
318 double Step_interp::rawinterp(
int j,
double x) {
319 if(xx[j] == xx[j + 1]) {
327 n(a.size()), lu(a), aref(a), indx(n) {
329 const double TINY = 1.0e-40;
330 int i = 0, imax = 0, j = 0, k = 0;
331 double big = 0., temp = 0.;
334 for(i = 0; i < n; i++) {
336 for(j = 0; j < n; j++)
337 if((temp = fabs(lu[i][j])) > big)
340 throw(
"Singular matrix in LUdcmp");
343 for(k = 0; k < n; k++) {
345 for(i = k; i < n; i++) {
346 temp = vv[i] * fabs(lu[i][k]);
353 for(j = 0; j < n; j++) {
355 lu[imax][j] = lu[k][j];
364 for(i = k + 1; i < n; i++) {
365 temp = lu[i][k] /= lu[k][k];
366 for(j = k + 1; j < n; j++)
367 lu[i][j] -= temp * lu[k][j];
374 int i, ii = 0, ip, j;
376 if(b.size() != n || x.size() != n)
377 throw(
"LUdcmp::solve bad sizes");
378 for(i = 0; i < n; i++)
380 for(i = 0; i < n; i++) {
385 for(j = ii - 1; j < i; j++)
386 sum -= lu[i][j] * x[j];
392 for(i = n - 1; i >= 0; i--) {
394 for(j = i + 1; j < n; j++)
395 sum -= lu[i][j] * x[j];
396 x[i] = sum / lu[i][i];