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