OpenWAM
Math_wam.cpp
Go to the documentation of this file.
1 /* --------------------------------------------------------------------------------*\
2 ==========================|
3  \\ /\ /\ // O pen | OpenWAM: The Open Source 1D Gas-Dynamic Code
4  \\ | X | // W ave |
5  \\ \/_\/ // A ction | CMT-Motores Termicos / Universidad Politecnica Valencia
6  \\/ \// M odel |
7  ----------------------------------------------------------------------------------
8  License
9 
10  This file is part of OpenWAM.
11 
12  OpenWAM is free software: you can redistribute it and/or modify
13  it under the terms of the GNU General Public License as published by
14  the Free Software Foundation, either version 3 of the License, or
15  (at your option) any later version.
16 
17  OpenWAM is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with OpenWAM. If not, see <http://www.gnu.org/licenses/>.
24 
25 
26  \*-------------------------------------------------------------------------------- */
27 
53 // ---------------------------------------------------------------------------
54 #include "Math_wam.h"
55 
56 double Interpola(double vizq, double vder, double axid, double xif) {
57  return vizq + ((vder - vizq) / axid) * xif;
58 }
59 
60 double QuadraticEqP(double A, double B, double C) {
61  return (-B + sqrt(B * B - 4 * A * C)) / 2 / A;
62 }
63 
64 double QuadraticEqN(double A, double B, double C) {
65  return (-B - sqrt(B * B - 4 * A * C)) / 2 / A;
66 }
67 
68 float Max(const double &a, const float &b) {
69  return b > a ? (b) : float(a);
70 }
71 
72 float Max(const float &a, const double &b) {
73  return b > a ? float(b) : (a);
74 }
75 
76 float Min(const double &a, const float &b) {
77  return b < a ? (b) : float(a);
78 }
79 
80 float Min(const float &a, const double &b) {
81  return b < a ? float(b) : (a);
82 }
83 
84 float Sign(const float &a, const double &b) {
85  return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
86 }
87 
88 float Sign(const double &a, const float &b) {
89  return (float)(b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));
90 }
91 
92 stPolar::stPolar() {
93 }
94 
95 stPolar::stPolar(double X, double Y) {
96  Mod = sqrt(X * X + Y * Y);
97  Ang = atan(Y / X);
98 }
99 
100 void stPolar::operator()(double X, double Y) {
101  Mod = sqrt(X * X + Y * Y);
102  Ang = atan(Y / X);
103 }
104 
105 stRectan::stRectan() {
106 }
107 
108 stRectan::stRectan(double Ang, double Mod) {
109  X = Mod * cos(Ang);
110  Y = Mod * sin(Ang);
111 }
112 
113 void stRectan::operator()(double Ang, double Mod) {
114  X = Mod * cos(Ang);
115  Y = Mod * sin(Ang);
116 }
117 
118 Base_interp::Base_interp() :
119  n(0), mm(0), jsav(0), cor(0), xx(), yy() {
120 }
121 
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));
125 }
126 
127 double Base_interp::interp(double x) {
128  int jlo = cor ? hunt(x) : locate(x);
129  return rawinterp(jlo, x);
130 }
131 
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]);
137  jl = 0;
138  ju = n - 1;
139  while(ju - jl > 1) {
140  jm = (ju + jl) >> 1;
141  if((x >= xx[jm]) == ascnd) {
142  jl = jm;
143  } else {
144  ju = jm;
145  }
146  }
147  cor = abs(jl - jsav) > dj ? 0 : 1;
148  jsav = jl;
149  return Max(0, Min(n - mm, jl - ((mm - 2) >> 1)));
150 
151 }
152 
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) {
159  jl = 0;
160  ju = n - 1;
161  } else {
162  if((x >= xx[jl]) == ascnd) {
163  for(;;) {
164  ju = jl + inc;
165  if(ju >= n) {
166  ju = n - 1;
167  break;
168  } else if((x < xx[ju]) == ascnd)
169  break;
170  else {
171  jl = ju;
172  inc += inc;
173  }
174  }
175  } else {
176  ju = jl;
177  for(;;) {
178  jl = jl - inc;
179  if(jl <= 0) {
180  jl = 0;
181  break;
182  } else if((x >= xx[jl]) == ascnd)
183  break;
184  else {
185  ju = jl;
186  inc += inc;
187  }
188  }
189  }
190  }
191  while(ju - jl > 1) {
192  jm = (ju + jl) >> 1;
193  if((x >= xx[jm]) == ascnd) {
194  jl = jm;
195  } else {
196  ju = jm;
197  }
198  }
199  cor = abs(jl - jsav) > dj ? 0 : 1;
200  jsav = jl;
201  return Max(0, Min(n - mm, jl - ((mm - 2) >> 1)));
202 }
203 
204 Linear_interp::Linear_interp() :
205  Base_interp() {
206 }
207 
208 Linear_interp::Linear_interp(dVector &xv, dVector &yv) :
209  Base_interp(xv, &yv[0], 2) {
210 }
211 
212 void Linear_interp::operator()(dVector & xv, dVector & yv) {
213  xx = &xv[0];
214  yy = &yv[0];
215  n = xv.size();
216  mm = 2;
217  jsav = 0;
218  cor = 0;
219  dj = Max(1, (int) pow025((double) n));
220 }
221 
222 double Linear_interp::rawinterp(int j, double x) {
223  if(xx[j] == xx[j + 1]) {
224  return yy[j];
225  } else {
226  return yy[j] + ((x - xx[j]) / (xx[j + 1] - xx[j])) * (yy[j + 1] - yy[j]);
227  }
228 }
229 
230 Hermite_interp::Hermite_interp() :
231  Base_interp(), y2() {
232 }
233 ;
234 
235 Hermite_interp::Hermite_interp(dVector &xv, dVector &yv) :
236  Base_interp(xv, &yv[0], 2), y2(xv.size()) {
237  sety2(&xv[0], &yv[0]);
238 }
239 
240 void Hermite_interp::operator()(dVector & xv, dVector & yv) {
241  xx = &xv[0];
242  yy = &yv[0];
243  n = xv.size();
244  mm = 2;
245  jsav = 0;
246  cor = 0;
247  y2.resize(n);
248  dj = Max(1, (int) pow025((double) n));
249  sety2(&xv[0], &yv[0]);
250 }
251 
252 void Hermite_interp::sety2(const double *xv, const double *yv) {
253  double DeltaK = 0., AlphaK = 0., BetaK = 0., TauK = 0.;
254 
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]);
257  }
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]);
260 
261  for(int i = 0; i < n - 1; i++) {
262  DeltaK = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]);
263  if(DeltaK == 0) {
264  y2[i] = 0;
265  y2[i + 1] = 0;
266  } else {
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;
273  }
274  }
275  }
276 }
277 
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.;
280  // int k=0;
281 
282  if(x <= xx[j]) {
283  ret_val = yy[j];
284  } else if(x >= xx[j + 1]) {
285  ret_val = yy[j + 1];
286  } else {
287  h = (xx[j + 1] - xx[j]);
288  t = (x - xx[j]) / h;
289  t2 = t * t;
290  t3 = t2 * t;
291  h00 = 2 * t3 - 3 * t2 + 1;
292  h10 = t3 - 2 * t2 + t;
293  h01 = -2 * t3 + 3 * t2;
294  h11 = t3 - t2;
295  ret_val = h00 * yy[j] + h * h10 * y2[j] + h01 * yy[j + 1] + h * h11 * y2[j + 1];
296  }
297  return ret_val;
298 }
299 
300 Step_interp::Step_interp() :
301  Base_interp() {
302 }
303 
304 Step_interp::Step_interp(dVector &xv, dVector &yv) :
305  Base_interp(xv, &yv[0], 2) {
306 }
307 
308 void Step_interp::operator()(dVector & xv, dVector & yv) {
309  xx = &xv[0];
310  yy = &yv[0];
311  n = xv.size();
312  mm = 2;
313  jsav = 0;
314  cor = 0;
315  dj = Max(1, (int) pow025((double) n));
316 }
317 
318 double Step_interp::rawinterp(int j, double x) {
319  if(xx[j] == xx[j + 1]) {
320  return yy[j];
321  } else {
322  return yy[j];
323  }
324 }
325 
326 LUdcmp::LUdcmp(dMatrix &a) :
327  n(a.size()), lu(a), aref(a), indx(n) {
328 
329  const double TINY = 1.0e-40;
330  int i = 0, imax = 0, j = 0, k = 0;
331  double big = 0., temp = 0.;
332  dVector vv(n);
333  d = 1.0;
334  for(i = 0; i < n; i++) {
335  big = 0.0;
336  for(j = 0; j < n; j++)
337  if((temp = fabs(lu[i][j])) > big)
338  big = temp;
339  if(big == 0.0)
340  throw("Singular matrix in LUdcmp");
341  vv[i] = 1.0 / big;
342  }
343  for(k = 0; k < n; k++) {
344  big = 0.0;
345  for(i = k; i < n; i++) {
346  temp = vv[i] * fabs(lu[i][k]);
347  if(temp > big) {
348  big = temp;
349  imax = i;
350  }
351  }
352  if(k != imax) {
353  for(j = 0; j < n; j++) {
354  temp = lu[imax][j];
355  lu[imax][j] = lu[k][j];
356  lu[k][j] = temp;
357  }
358  d = -d;
359  vv[imax] = vv[k];
360  }
361  indx[k] = imax;
362  if(lu[k][k] == 0.0)
363  lu[k][k] = TINY;
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];
368 
369  }
370  }
371 }
372 
373 void LUdcmp::solve(dVector &b, dVector &x) {
374  int i, ii = 0, ip, j;
375  double sum = 0.;
376  if(b.size() != n || x.size() != n)
377  throw("LUdcmp::solve bad sizes");
378  for(i = 0; i < n; i++)
379  x[i] = b[i];
380  for(i = 0; i < n; i++) {
381  ip = indx[i];
382  sum = x[ip];
383  x[ip] = x[i];
384  if(ii != 0)
385  for(j = ii - 1; j < i; j++)
386  sum -= lu[i][j] * x[j];
387  else if(sum != 0.0)
388  ii = i + 1;
389  x[i] = sum;
390 
391  }
392  for(i = n - 1; i >= 0; i--) {
393  sum = x[i];
394  for(j = i + 1; j < n; j++)
395  sum -= lu[i][j] * x[j];
396  x[i] = sum / lu[i][i];
397  }
398 
399 }
400 
pow025
T pow025(T x)
Returns x to the power of 0.25.
Definition: Math_wam.h:127
dMatrix
std::vector< std::vector< double > > dMatrix
2-dimensional double matrix
Definition: Math_wam.h:71
Base_interp
Definition: Math_wam.h:285
Math_wam.h
dVector
std::vector< double > dVector
Double vector.
Definition: Math_wam.h:70