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 
103 namespace 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
Lorene prototypes.
Definition: app_hor.h:67
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::extract_row
Tbl extract_row(const Tbl &xytab, const Tbl &xtab, double xx)
Extraction of a row of a 2D Tbl
Definition: interpol_herm.C:631
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Lorene::Tbl::get_ndim
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
Definition: tbl.h:400
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
Lorene::extract_column
Tbl extract_column(const Tbl &, const Tbl &, double)
Extraction of a column of a 2D Tbl
Definition: interpol_herm.C:615
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
Lorene::Tbl
Basic array class.
Definition: tbl.h:161