LORENE
interpol_herm.C
1/*
2 * Hermite interpolation functions.
3 *
4 */
5
6/*
7 * Copyright (c) 2000-2002 Eric Gourgoulhon
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation; either version 2 of the License, or
14 * (at your option) any later version.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27
28
29
30/*
31 * $Id: interpol_herm.C,v 1.18 2022/03/22 13:19:53 g_servignat Exp $
32 * $Log: interpol_herm.C,v $
33 * Revision 1.18 2022/03/22 13:19:53 g_servignat
34 * Modified row and column extraction of 2D Tbl's
35 *
36 * Revision 1.17 2022/03/01 10:03:38 g_servignat
37 * Added subtbl extraction for interpol_linear_2D purposes
38 *
39 * Revision 1.16 2022/02/25 10:45:21 g_servignat
40 * Added 2D linear interpolation
41 *
42 * Revision 1.15 2021/05/14 15:39:23 g_servignat
43 * Added sound speed computation from enthalpy to Eos class and tabulated+polytropic derived classes
44 *
45 * Revision 1.14 2016/12/05 16:18:11 j_novak
46 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
47 *
48 * Revision 1.13 2015/06/15 15:08:22 j_novak
49 * New file interpol_bifluid for interpolation of 2-fluid EoSs
50 *
51 * Revision 1.12 2015/06/10 14:39:18 a_sourie
52 * New class Eos_bf_tabul for tabulated 2-fluid EoSs and associated functions for the computation of rotating stars with such EoSs.
53 *
54 * Revision 1.11 2015/01/27 14:22:38 j_novak
55 * New methods in Eos_tabul to correct for EoS themro consistency (optional).
56 *
57 * Revision 1.10 2014/10/13 08:53:32 j_novak
58 * Lorene classes and functions now belong to the namespace Lorene.
59 *
60 * Revision 1.9 2013/12/12 16:07:30 j_novak
61 * interpol_herm_2d outputs df/dx, used to get the magnetization.
62 *
63 * Revision 1.8 2012/09/04 14:53:28 j_novak
64 * Replacement of the FORTRAN version of huntm by a C one.
65 *
66 * Revision 1.7 2011/10/04 16:05:19 j_novak
67 * Update of Eos_mag class. Suppression of loge, re-definition of the derivatives
68 * and use of interpol_herm_2d.
69 *
70 * Revision 1.6 2011/10/03 13:44:45 j_novak
71 * Updated the y-derivative for the 2D version
72 *
73 * Revision 1.5 2011/09/27 15:38:11 j_novak
74 * New function for 2D interpolation added. The computation of 1st derivative is
75 * still missing.
76 *
77 * Revision 1.4 2003/11/21 16:14:51 m_bejger
78 * Added the linear interpolation
79 *
80 * Revision 1.3 2003/05/15 09:42:12 e_gourgoulhon
81 * Added the new function interpol_herm_der
82 *
83 * Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
84 * Modification of declaration of Fortran 77 prototypes for
85 * a better portability (in particular on IBM AIX systems):
86 * All Fortran subroutine names are now written F77_* and are
87 * defined in the new file C++/Include/proto_f77.h.
88 *
89 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
90 * LORENE
91 *
92 * Revision 2.0 2000/11/22 19:31:42 eric
93 * *** empty log message ***
94 *
95 *
96 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/interpol_herm.C,v 1.18 2022/03/22 13:19:53 g_servignat Exp $
97 *
98 */
99
100// Headers Lorene
101#include "tbl.h"
102
103namespace Lorene {
104
105 //---------------------------------------------------------------
106 // Value bracketting in an ordered table (from Numerical Recipes)
107 //---------------------------------------------------------------
108 void huntm(const Tbl& xx, double& x, int& i_low) {
109
110 assert (xx.get_etat() == ETATQCQ) ;
111 int nx = xx.get_taille() ;
112 bool ascend = ( xx(nx-1) > xx(0) ) ;
113 int i_hi ;
114 if ( (i_low < 0)||(i_low>=nx) ) {
115 i_low = -1 ;
116 i_hi = nx ;
117 }
118 else {
119 int inc = 1 ;
120 if ( (x >= xx(i_low)) == ascend ) {
121 if (i_low == nx -1) return ;
122 i_hi = i_low + 1 ;
123 while ( (x >= xx(i_hi)) == ascend ) {
124 i_low = i_hi ;
125 inc += inc ;
126 i_hi = i_low + inc ;
127 if (i_hi >= nx) {
128 i_hi = nx ;
129 break ;
130 }
131 }
132 } else {
133 if (i_low == 0) {
134 i_low = -1 ;
135 return ;
136 }
137 i_hi = i_low-- ;
138 while ( (x < xx(i_low)) == ascend ) {
139 i_hi = i_low ;
140 inc += inc ;
141 if ( inc >= i_hi ) {
142 i_low = 0 ;
143 break ;
144 }
145 else i_low = i_hi - inc ;
146 }
147 }
148 }
149 while ( (i_hi - i_low) > 1) {
150 int i_med = (i_hi + i_low) / 2 ;
151 if ( (x>=xx(i_med)) == ascend ) i_low = i_med ;
152 else i_hi = i_med ;
153 }
154 if (x == xx(nx-1)) i_low = nx-2 ;
155 if (x == xx(0)) i_low = 0 ;
156 return ;
157 }
158
159 //---------------------
160 // Linear interpolation
161 //---------------------
162 void interpol_linear(const Tbl& xtab, const Tbl& ytab,
163 double x, int& i, double& y) {
164
165 assert(ytab.dim == xtab.dim) ;
166 //assert(dytab.dim == xtab.dim) ;
167
168 huntm(xtab, x, i) ;
169
170 int i1 = i + 1 ;
171
172 // double dx = xtab(i1) - xtab(i) ;
173 double y1 = ytab(i) ;
174 double y2 = ytab(i1) ;
175
176 double x1 = xtab(i) ;
177 double x2 = xtab(i1) ;
178 double x12 = x1-x2 ;
179
180 double a = (y1-y2)/x12 ;
181 double b = (x1*y2-y1*x2)/x12 ;
182
183 y = x*a+b ;
184
185 }
186
187 //---------------------
188 // Linear interpolation
189 //---------------------
190 void interpol_linear_2d(const Tbl& x1tab, const Tbl& x2tab, const Tbl& ytab,
191 double x1, double x2, int& i, int& j, double& y) {
192
193 assert(ytab.dim.ndim == 2) ;
194 assert(x1tab.dim.ndim == 1) ;
195 assert(x2tab.dim.ndim == 1) ;
196 assert(ytab.dim.dim[1] == x2tab.dim.dim[0]) ;
197 assert(ytab.dim.dim[0] == x1tab.dim.dim[0]) ;
198
199 huntm(x1tab, x1, i) ;
200 huntm(x2tab, x2, j) ;
201
202 int i1 = i + 1 ;
203 int j1 = j + 1 ;
204
205 // double dx = xtab(i1) - xtab(i) ;
206 double y11 = ytab(j,i) ;
207 double y21 = ytab(j1,i) ;
208 double y12 = ytab(j,i1) ;
209
210 double x11 = x1tab(i) ;
211 double x12 = x1tab(i1) ;
212 double x21 = x2tab(j) ;
213 double x22 = x2tab(j1) ;
214 double x1diff = x12-x11 ;
215 double x2diff = x22-x21 ;
216
217 double a = (y21-y11)/x2diff ;
218 double b = (y12-y11)/x1diff ;
219
220 y = y11 + (x1-x11)*b + (x2-x21)*a ;
221
222 }
223
224 //------------------------
225 // Quadratic interpolation
226 //------------------------
227 void interpol_quad(const Tbl& xtab, const Tbl& ytab,
228 double x, int& i, double& y) {
229
230 assert(ytab.dim == xtab.dim) ;
231 //assert(dytab.dim == xtab.dim) ;
232
233 huntm(xtab, x, i) ;
234
235 int i1 = i - 1 ;
236 int i2 = i + 1 ;
237
238 double y0=ytab(i1) ;
239 double y1=ytab(i) ;
240 double y2=ytab(i2) ;
241
242 double x0=xtab(i1) ;
243 double x1=xtab(i) ;
244 double x2=xtab(i2) ;
245 double x01=x0-x1 ;
246 double x02=x0-x2 ;
247 double x12=x1-x2 ;
248
249 y = y0*(x-x1)*(x-x2)/(x01*x02) - y1*(x-x0)*(x-x2)/(x01*x12) + y2*(x-x0)*(x-x1)/(x02*x12) ;
250
251 }
252
253 //------------------------------------------------------------
254 // Cubic Hermite interpolation, returning the first derivative
255 //------------------------------------------------------------
256 void interpol_herm(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
257 double x, int& i, double& y, double& dy) {
258
259 assert(ytab.dim == xtab.dim) ;
260 assert(dytab.dim == xtab.dim) ;
261
262 huntm(xtab, x, i) ;
263
264 int i1 = i + 1 ;
265
266 double dx = xtab(i1) - xtab(i) ;
267
268 double u = (x - xtab(i)) / dx ;
269 double u2 = u*u ;
270 double u3 = u2*u ;
271
272 y = ytab(i) * ( 2.*u3 - 3.*u2 + 1.)
273 + ytab(i1) * ( 3.*u2 - 2.*u3)
274 + dytab(i) * dx * ( u3 - 2.*u2 + u )
275 - dytab(i1) * dx * ( u2 - u3 ) ;
276
277 dy = 6. * ( ytab(i) / dx * ( u2 - u )
278 - ytab(i1) / dx * ( u2 - u ) )
279 + dytab(i) * ( 3.*u2 - 4.*u + 1. )
280 + dytab(i1) * ( 3.*u2 - 2.*u ) ;
281 }
282
283
284 //-------------------------------------------------------------
285 // Cubic Hermite interpolation, returning the second derivative
286 //-------------------------------------------------------------
287 void interpol_herm_der(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
288 double x, int& i, double& y, double& dy, double& ddy) {
289
290 assert(ytab.dim == xtab.dim) ;
291 assert(dytab.dim == xtab.dim) ;
292
293 huntm(xtab, x, i) ;
294
295 // i-- ; // Fortran --> C
296
297 int i1 = i + 1 ;
298
299 double dx = xtab(i1) - xtab(i) ;
300
301 double u = (x - xtab(i)) / dx ;
302 double u2 = u*u ;
303 double u3 = u2*u ;
304
305 y = ytab(i) * ( 2.*u3 - 3.*u2 + 1.)
306 + ytab(i1) * ( 3.*u2 - 2.*u3)
307 + dytab(i) * dx * ( u3 - 2.*u2 + u )
308 - dytab(i1) * dx * ( u2 - u3 ) ;
309
310 dy = 6. * ( ytab(i) - ytab(i1) ) * ( u2 - u ) / dx
311 + dytab(i) * ( 3.*u2 - 4.*u + 1. )
312 + dytab(i1) * ( 3.*u2 - 2.*u ) ;
313
314 ddy = 6 * ( ( ytab(i) - ytab(i1) ) * ( 2.*u - 1. ) / dx
315 + dytab(i) * (6.*u - 4.)
316 + dytab(i1) * (6.*u - 2.) ) / dx ;
317
318 }
319
320 //----------------------------------------------
321 // Bi-cubic Hermite interpolation, for 2D arrays
322 //----------------------------------------------
323 void interpol_herm_2d(const Tbl& xtab, const Tbl& ytab, const Tbl& ftab,
324 const Tbl& dfdxtab, const Tbl& dfdytab, const Tbl&
325 d2fdxdytab, double x, double y, double& f, double&
326 dfdx, double& dfdy) {
327
328 assert(ytab.dim == xtab.dim) ;
329 assert(ftab.dim == xtab.dim) ;
330 assert(dfdxtab.dim == xtab.dim) ;
331 assert(dfdytab.dim == xtab.dim) ;
332 assert(d2fdxdytab.dim == xtab.dim) ;
333
334 int nbp1, nbp2;
335 nbp1 = xtab.get_dim(0);
336 nbp2 = xtab.get_dim(1);
337
338 int i_near = 0 ;
339 int j_near = 0 ;
340
341 while ((xtab(i_near,0) <= x) && (nbp2 > i_near)) {
342 i_near++;
343 }
344 if (i_near != 0) {
345 i_near-- ;
346 }
347 j_near = 0;
348 while ((ytab(i_near,j_near) < y) && (nbp1 > j_near)) {
349 j_near++ ;
350 }
351 if (j_near != 0) {
352 j_near-- ;
353 }
354
355 int i1 = i_near+1 ; int j1 = j_near+1 ;
356
357 double dx = xtab(i1, j_near) - xtab(i_near, j_near) ;
358 double dy = ytab(i_near, j1) - ytab(i_near, j_near) ;
359
360 double u = (x - xtab(i_near, j_near)) / dx ;
361 double v = (y - ytab(i_near, j_near)) / dy ;
362
363 double u2 = u*u ; double v2 = v*v ;
364 double u3 = u2*u ; double v3 = v2*v ;
365
366 double psi0_u = 2.*u3 - 3.*u2 + 1. ;
367 double psi0_1mu = -2.*u3 + 3.*u2 ;
368 double psi1_u = u3 - 2.*u2 + u ;
369 double psi1_1mu = -u3 + u2 ;
370
371 double psi0_v = 2.*v3 - 3.*v2 + 1. ;
372 double psi0_1mv = -2.*v3 + 3.*v2 ;
373 double psi1_v = v3 - 2.*v2 + v ;
374 double psi1_1mv = -v3 + v2 ;
375
376 f = ftab(i_near, j_near) * psi0_u * psi0_v
377 + ftab(i1, j_near) * psi0_1mu * psi0_v
378 + ftab(i_near, j1) * psi0_u * psi0_1mv
379 + ftab(i1, j1) * psi0_1mu * psi0_1mv ;
380
381 f += (dfdxtab(i_near, j_near) * psi1_u * psi0_v
382 - dfdxtab(i1, j_near) * psi1_1mu * psi0_v
383 + dfdxtab(i_near, j1) * psi1_u * psi0_1mv
384 - dfdxtab(i1, j1) * psi1_1mu * psi0_1mv) * dx ;
385
386 f += (dfdytab(i_near, j_near) * psi0_u * psi1_v
387 + dfdytab(i1, j_near) * psi0_1mu * psi1_v
388 - dfdytab(i_near, j1) * psi0_u * psi1_1mv
389 - dfdytab(i1, j1) * psi0_1mu * psi1_1mv) * dy ;
390
391 f += (d2fdxdytab(i_near, j_near) * psi1_u * psi1_v
392 - d2fdxdytab(i1, j_near) * psi1_1mu * psi1_v
393 - d2fdxdytab(i_near, j1) * psi1_u * psi1_1mv
394 + d2fdxdytab(i1, j1) * psi1_1mu * psi1_1mv) * dx * dy ;
395
396 double dpsi0_u = 6.*(u2 - u) ;
397 double dpsi0_1mu = 6.*(u2 - u) ;
398 double dpsi1_u = 3.*u2 - 4.*u + 1. ;
399 double dpsi1_1mu = 3.*u2 - 2.*u ;
400
401 dfdx = (ftab(i_near, j_near) * dpsi0_u * psi0_v
402 - ftab(i1, j_near) * dpsi0_1mu * psi0_v
403 + ftab(i_near, j1) * dpsi0_u * psi0_1mv
404 - ftab(i1, j1) * dpsi0_1mu * psi0_1mv ) / dx;
405
406 dfdx += (dfdxtab(i_near, j_near) * dpsi1_u * psi0_v
407 + dfdxtab(i1, j_near) * dpsi1_1mu * psi0_v
408 + dfdxtab(i_near, j1) * dpsi1_u * psi0_1mv
409 + dfdxtab(i1, j1) * dpsi1_1mu * psi0_1mv) ;
410
411 dfdx += (dfdytab(i_near, j_near) * dpsi0_u * psi1_v
412 - dfdytab(i1, j_near) * dpsi0_1mu * psi1_v
413 - dfdytab(i_near, j1) * dpsi0_u * psi1_1mv
414 + dfdytab(i1, j1) * dpsi0_1mu * psi1_1mv) * dy /dx ;
415
416 dfdx += (d2fdxdytab(i_near, j_near) * dpsi1_u * psi1_v
417 + d2fdxdytab(i1, j_near) * dpsi1_1mu * psi1_v
418 - d2fdxdytab(i_near, j1) * dpsi1_u * psi1_1mv
419 - d2fdxdytab(i1, j1) * dpsi1_1mu * psi1_1mv) * dy ;
420
421 double dpsi0_v = 6.*(v2 - v) ;
422 double dpsi0_1mv = 6.*(v2 - v) ;
423 double dpsi1_v = 3.*v2 - 4.*v + 1. ;
424 double dpsi1_1mv = 3.*v2 - 2.*v ;
425
426 dfdy = (ftab(i_near, j_near) * psi0_u * dpsi0_v
427 + ftab(i1, j_near) * psi0_1mu * dpsi0_v
428 - ftab(i_near, j1) * psi0_u * dpsi0_1mv
429 - ftab(i1, j1) * psi0_1mu * dpsi0_1mv) / dy ;
430
431 dfdy += (dfdxtab(i_near, j_near) * psi1_u * dpsi0_v
432 - dfdxtab(i1, j_near) * psi1_1mu * dpsi0_v
433 - dfdxtab(i_near, j1) * psi1_u * dpsi0_1mv
434 + dfdxtab(i1, j1) * psi1_1mu * dpsi0_1mv) * dx / dy ;
435
436 dfdy += (dfdytab(i_near, j_near) * psi0_u * dpsi1_v
437 + dfdytab(i1, j_near) * psi0_1mu * dpsi1_v
438 + dfdytab(i_near, j1) * psi0_u * dpsi1_1mv
439 + dfdytab(i1, j1) * psi0_1mu * dpsi1_1mv) ;
440
441 dfdy += (d2fdxdytab(i_near, j_near) * psi1_u * dpsi1_v
442 - d2fdxdytab(i1, j_near) * psi1_1mu * dpsi1_v
443 + d2fdxdytab(i_near, j1) * psi1_u * dpsi1_1mv
444 - d2fdxdytab(i1, j1) * psi1_1mu * dpsi1_1mv) * dx ;
445
446 return ;
447 }
448
449
450
451 void interpol_herm_2d_sans(const Tbl& xtab, const Tbl& ytab, const Tbl& ftab,
452 const Tbl& dfdxtab, const Tbl& dfdytab, double x,
453 double y, double& f, double& dfdx, double& dfdy) {
454
455 assert(ytab.dim == xtab.dim) ;
456 assert(ftab.dim == xtab.dim) ;
457 assert(dfdxtab.dim == xtab.dim) ;
458 assert(dfdytab.dim == xtab.dim) ;
459
460 int nbp1, nbp2;
461 nbp1 = xtab.get_dim(0);
462 nbp2 = xtab.get_dim(1);
463
464 int i_near = 0 ;
465 int j_near = 0 ;
466
467 while ((xtab(i_near,0) <= x) && (nbp2 > i_near)) {
468 i_near++;
469 }
470 if (i_near != 0) {
471 i_near-- ;
472 }
473 j_near = 0;
474 while ((ytab(i_near,j_near) < y) && (nbp1 > j_near)) {
475 j_near++ ;
476 }
477 if (j_near != 0) {
478 j_near-- ;
479 }
480
481 int i1 = i_near+1 ; int j1 = j_near+1 ;
482
483 double dx = xtab(i1, j_near) - xtab(i_near, j_near) ;
484 double dy = ytab(i_near, j1) - ytab(i_near, j_near) ;
485
486 double u = (x - xtab(i_near, j_near)) / dx ;
487 double v = (y - ytab(i_near, j_near)) / dy ;
488
489 double u2 = u*u ; double v2 = v*v ;
490 double u3 = u2*u ; double v3 = v2*v ;
491
492 double psi0_u = 2.*u3 - 3.*u2 + 1. ;
493 double psi0_1mu = -2.*u3 + 3.*u2 ;
494 double psi1_u = u3 - 2.*u2 + u ;
495 double psi1_1mu = -u3 + u2 ;
496
497 double psi0_v = 2.*v3 - 3.*v2 + 1. ;
498 double psi0_1mv = -2.*v3 + 3.*v2 ;
499 double psi1_v = v3 - 2.*v2 + v ;
500 double psi1_1mv = -v3 + v2 ;
501
502 f = ftab(i_near, j_near) * psi0_u * psi0_v
503 + ftab(i1, j_near) * psi0_1mu * psi0_v
504 + ftab(i_near, j1) * psi0_u * psi0_1mv
505 + ftab(i1, j1) * psi0_1mu * psi0_1mv ;
506
507 f += (dfdxtab(i_near, j_near) * psi1_u * psi0_v
508 - dfdxtab(i1, j_near) * psi1_1mu * psi0_v
509 + dfdxtab(i_near, j1) * psi1_u * psi0_1mv
510 - dfdxtab(i1, j1) * psi1_1mu * psi0_1mv) * dx ;
511
512 f += (dfdytab(i_near, j_near) * psi0_u * psi1_v
513 + dfdytab(i1, j_near) * psi0_1mu * psi1_v
514 - dfdytab(i_near, j1) * psi0_u * psi1_1mv
515 - dfdytab(i1, j1) * psi0_1mu * psi1_1mv) * dy ;
516
517 double dpsi0_u = 6.*(u2 - u) ;
518 double dpsi0_1mu = 6.*(u2 - u) ;
519 double dpsi1_u = 3.*u2 - 4.*u + 1. ;
520 double dpsi1_1mu = 3.*u2 - 2.*u ;
521
522 dfdx = (ftab(i_near, j_near) * dpsi0_u * psi0_v
523 - ftab(i1, j_near) * dpsi0_1mu * psi0_v
524 + ftab(i_near, j1) * dpsi0_u * psi0_1mv
525 - ftab(i1, j1) * dpsi0_1mu * psi0_1mv ) / dx;
526
527 dfdx += (dfdxtab(i_near, j_near) * dpsi1_u * psi0_v
528 + dfdxtab(i1, j_near) * dpsi1_1mu * psi0_v
529 + dfdxtab(i_near, j1) * dpsi1_u * psi0_1mv
530 + dfdxtab(i1, j1) * dpsi1_1mu * psi0_1mv) ;
531
532 dfdx += (dfdytab(i_near, j_near) * dpsi0_u * psi1_v
533 - dfdytab(i1, j_near) * dpsi0_1mu * psi1_v
534 - dfdytab(i_near, j1) * dpsi0_u * psi1_1mv
535 + dfdytab(i1, j1) * dpsi0_1mu * psi1_1mv) * dy /dx ;
536
537 double dpsi0_v = 6.*(v2 - v) ;
538 double dpsi0_1mv = 6.*(v2 - v) ;
539 double dpsi1_v = 3.*v2 - 4.*v + 1. ;
540 double dpsi1_1mv = 3.*v2 - 2.*v ;
541
542 dfdy = (ftab(i_near, j_near) * psi0_u * dpsi0_v
543 + ftab(i1, j_near) * psi0_1mu * dpsi0_v
544 - ftab(i_near, j1) * psi0_u * dpsi0_1mv
545 - ftab(i1, j1) * psi0_1mu * dpsi0_1mv) / dy ;
546
547 dfdy += (dfdxtab(i_near, j_near) * psi1_u * dpsi0_v
548 - dfdxtab(i1, j_near) * psi1_1mu * dpsi0_v
549 - dfdxtab(i_near, j1) * psi1_u * dpsi0_1mv
550 + dfdxtab(i1, j1) * psi1_1mu * dpsi0_1mv) * dx / dy ;
551
552 dfdy += (dfdytab(i_near, j_near) * psi0_u * dpsi1_v
553 + dfdytab(i1, j_near) * psi0_1mu * dpsi1_v
554 + dfdytab(i_near, j1) * psi0_u * dpsi1_1mv
555 + dfdytab(i1, j1) * psi0_1mu * dpsi1_1mv) ;
556
557 return ;
558 }
559
560 //--------------------------------------------------------------------
561 // Quintic Hermite interpolation using data from the second derivative
562 //--------------------------------------------------------------------
563 void interpol_herm_2nd_der(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
564 const Tbl& d2ytab, double x, int& i, double& y,
565 double& dy, double& d2y) {
566
567 assert(ytab.dim == xtab.dim) ;
568 assert(dytab.dim == xtab.dim) ;
569 assert(d2ytab.dim == xtab.dim) ;
570
571 huntm(xtab, x, i) ;
572
573 int i1 = i + 1 ;
574
575 double dx = xtab(i1) - xtab(i) ;
576
577 double u = (x - xtab(i)) / dx ;
578 double u2 = u*u ;
579 double u3 = u2*u ;
580 double u4 = u2*u2 ;
581 double u5 = u3*u2 ;
582
583 double v = 1. - u ;
584 double v2 = v*v ;
585 double v3 = v2*v ;
586 double v4 = v2*v2 ;
587 double v5 = v3*v2 ;
588
589 y = ytab(i) * ( -6.*u5 + 15.*u4 - 10.*u3 + 1. )
590 + ytab(i1) * ( -6.*v5 + 15.*v4 - 10.*v3 + 1. )
591 + dytab(i) * dx * ( -3.*u5 + 8.*u4 -6.*u3 + u )
592 - dytab(i1) * dx * ( -3.*v5 + 8.*v4 -6.*v3 + v )
593 + d2ytab(i) * dx*dx * ( -0.5*u5 + 1.5*u4 - 1.5*u3 + 0.5*u2 )
594 + d2ytab(i1) * dx*dx * ( -0.5*v5 + 1.5*v4 - 1.5*v3 + 0.5*v2 ) ;
595
596 dy = 30.*( ytab(i) / dx * ( -u4 + 2.*u3 - u2 )
597 - ytab(i1) / dx * ( -v4 + 2.*v3 - v2 ) )
598 + dytab(i) * ( -15.*u4 + 32.*u3 - 18.*u2 + 1. )
599 + dytab(i1) * ( -15.*v4 + 32.*v3 - 18.*v2 + 1. )
600 + d2ytab(i) * dx * ( -2.5*u4 + 6.*u3 -4.5*u2 + u )
601 - d2ytab(i1) * dx * ( -2.5*v4 + 6.*v3 -4.5*v2 + v ) ;
602
603 d2y = 60.*(ytab(i)/(dx*dx) * (-2.*u3 + 3*u2 - u)
604 +ytab(i1)/(dx*dx) *(-2.*v3 + 3*v2 - v))
605 + 12.*dytab(i)/dx *(-5.*u3 + 8.*u2 - 3.*u)
606 - 12.*dytab(i1)/dx *(-5.*v3 + 8.*v2 - 3.*v)
607 + d2ytab(i) *(-10.*u3 + 18.*u2 - 9.*u + 1.)
608 + d2ytab(i1) *(-10.*v3 + 18.*v2 - 9.*v + 1.) ;
609 }
610
611
615 Tbl extract_column(const Tbl& xytab, const Tbl& ytab, double yy){
616 assert(xytab.get_ndim() == 2) ;
617 assert(ytab.get_ndim() == 1) ;
618 int i_low = ytab.get_taille()/2;
619 huntm(ytab,yy,i_low) ;
620 int n_h = xytab.get_dim(0) ;
621 Tbl resu(n_h) ; resu.set_etat_qcq() ;
622 for (int i=0 ; i<n_h ; i++){
623 resu.set(i) = xytab(i_low,i) ;
624 }
625 return resu;
626 }
627
631 Tbl extract_row(const Tbl& xytab, const Tbl& xtab, double xx){
632 assert(xytab.get_ndim() == 2) ;
633 assert(xtab.get_ndim() == 1) ;
634 int i_low = xtab.get_taille()/2;
635 huntm(xtab,xx,i_low) ;
636 int n_h = xytab.get_dim(1) ;
637 Tbl resu(n_h) ; resu.set_etat_qcq() ;
638 for (int i=0 ; i<n_h ; i++){
639 resu.set(i) = xytab(i,i_low) ;
640 }
641 return resu;
642 }
643
644
645} // End of namespace Lorene
646
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim).
Definition tbl.h:400
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille).
Definition tbl.h:397
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i]).
Definition tbl.h:403
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::extract_row
Tbl extract_row(const Tbl &xytab, const Tbl &xtab, double xx)
Definition interpol_herm.C:631
Lorene::y
Coord y
y coordinate centered on the grid
Definition map.h:739
Lorene::extract_column
Tbl extract_column(const Tbl &, const Tbl &, double)
Definition interpol_herm.C:615
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738