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 char interpol_herm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/interpol_herm.C,v 1.13 2015/06/15 15:08:22 j_novak Exp $" ;
29 
30 /*
31  * $Id: interpol_herm.C,v 1.13 2015/06/15 15:08:22 j_novak Exp $
32  * $Log: interpol_herm.C,v $
33  * Revision 1.13 2015/06/15 15:08:22 j_novak
34  * New file interpol_bifluid for interpolation of 2-fluid EoSs
35  *
36  * Revision 1.12 2015/06/10 14:39:18 a_sourie
37  * New class Eos_bf_tabul for tabulated 2-fluid EoSs and associated functions for the computation of rotating stars with such EoSs.
38  *
39  * Revision 1.11 2015/01/27 14:22:38 j_novak
40  * New methods in Eos_tabul to correct for EoS themro consistency (optional).
41  *
42  * Revision 1.10 2014/10/13 08:53:32 j_novak
43  * Lorene classes and functions now belong to the namespace Lorene.
44  *
45  * Revision 1.9 2013/12/12 16:07:30 j_novak
46  * interpol_herm_2d outputs df/dx, used to get the magnetization.
47  *
48  * Revision 1.8 2012/09/04 14:53:28 j_novak
49  * Replacement of the FORTRAN version of huntm by a C one.
50  *
51  * Revision 1.7 2011/10/04 16:05:19 j_novak
52  * Update of Eos_mag class. Suppression of loge, re-definition of the derivatives
53  * and use of interpol_herm_2d.
54  *
55  * Revision 1.6 2011/10/03 13:44:45 j_novak
56  * Updated the y-derivative for the 2D version
57  *
58  * Revision 1.5 2011/09/27 15:38:11 j_novak
59  * New function for 2D interpolation added. The computation of 1st derivative is
60  * still missing.
61  *
62  * Revision 1.4 2003/11/21 16:14:51 m_bejger
63  * Added the linear interpolation
64  *
65  * Revision 1.3 2003/05/15 09:42:12 e_gourgoulhon
66  * Added the new function interpol_herm_der
67  *
68  * Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
69  * Modification of declaration of Fortran 77 prototypes for
70  * a better portability (in particular on IBM AIX systems):
71  * All Fortran subroutine names are now written F77_* and are
72  * defined in the new file C++/Include/proto_f77.h.
73  *
74  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
75  * LORENE
76  *
77  * Revision 2.0 2000/11/22 19:31:42 eric
78  * *** empty log message ***
79  *
80  *
81  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/interpol_herm.C,v 1.13 2015/06/15 15:08:22 j_novak Exp $
82  *
83  */
84 
85 // Headers Lorene
86 #include "tbl.h"
87 
88 namespace Lorene {
89 
90  //---------------------------------------------------------------
91  // Value bracketting in an ordered table (from Numerical Recipes)
92  //---------------------------------------------------------------
93  void huntm(const Tbl& xx, double& x, int& i_low) {
94 
95  assert (xx.get_etat() == ETATQCQ) ;
96  int nx = xx.get_taille() ;
97  bool ascend = ( xx(nx-1) > xx(0) ) ;
98  int i_hi ;
99  if ( (i_low < 0)||(i_low>=nx) ) {
100  i_low = -1 ;
101  i_hi = nx ;
102  }
103  else {
104  int inc = 1 ;
105  if ( (x >= xx(i_low)) == ascend ) {
106  if (i_low == nx -1) return ;
107  i_hi = i_low + 1 ;
108  while ( (x >= xx(i_hi)) == ascend ) {
109  i_low = i_hi ;
110  inc += inc ;
111  i_hi = i_low + inc ;
112  if (i_hi >= nx) {
113  i_hi = nx ;
114  break ;
115  }
116  }
117  } else {
118  if (i_low == 0) {
119  i_low = -1 ;
120  return ;
121  }
122  i_hi = i_low-- ;
123  while ( (x < xx(i_low)) == ascend ) {
124  i_hi = i_low ;
125  inc += inc ;
126  if ( inc >= i_hi ) {
127  i_low = 0 ;
128  break ;
129  }
130  else i_low = i_hi - inc ;
131  }
132  }
133  }
134  while ( (i_hi - i_low) > 1) {
135  int i_med = (i_hi + i_low) / 2 ;
136  if ( (x>=xx(i_med)) == ascend ) i_low = i_med ;
137  else i_hi = i_med ;
138  }
139  if (x == xx(nx-1)) i_low = nx-2 ;
140  if (x == xx(0)) i_low = 0 ;
141  return ;
142  }
143 
144  //---------------------
145  // Linear interpolation
146  //---------------------
147  void interpol_linear(const Tbl& xtab, const Tbl& ytab,
148  double x, int& i, double& y) {
149 
150  assert(ytab.dim == xtab.dim) ;
151  //assert(dytab.dim == xtab.dim) ;
152 
153  huntm(xtab, x, i) ;
154 
155  int i1 = i + 1 ;
156 
157  // double dx = xtab(i1) - xtab(i) ;
158  double y1 = ytab(i) ;
159  double y2 = ytab(i1) ;
160 
161  double x1 = xtab(i) ;
162  double x2 = xtab(i1) ;
163  double x12 = x1-x2 ;
164 
165  double a = (y1-y2)/x12 ;
166  double b = (x1*y2-y1*x2)/x12 ;
167 
168  y = x*a+b ;
169 
170  }
171 
172  //------------------------------------------------------------
173  // Cubic Hermite interpolation, returning the first derivative
174  //------------------------------------------------------------
175  void interpol_herm(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
176  double x, int& i, double& y, double& dy) {
177 
178  assert(ytab.dim == xtab.dim) ;
179  assert(dytab.dim == xtab.dim) ;
180 
181  huntm(xtab, x, i) ;
182 
183  int i1 = i + 1 ;
184 
185  double dx = xtab(i1) - xtab(i) ;
186 
187  double u = (x - xtab(i)) / dx ;
188  double u2 = u*u ;
189  double u3 = u2*u ;
190 
191  y = ytab(i) * ( 2.*u3 - 3.*u2 + 1.)
192  + ytab(i1) * ( 3.*u2 - 2.*u3)
193  + dytab(i) * dx * ( u3 - 2.*u2 + u )
194  - dytab(i1) * dx * ( u2 - u3 ) ;
195 
196  dy = 6. * ( ytab(i) / dx * ( u2 - u )
197  - ytab(i1) / dx * ( u2 - u ) )
198  + dytab(i) * ( 3.*u2 - 4.*u + 1. )
199  + dytab(i1) * ( 3.*u2 - 2.*u ) ;
200  }
201 
202 
203  //-------------------------------------------------------------
204  // Cubic Hermite interpolation, returning the second derivative
205  //-------------------------------------------------------------
206  void interpol_herm_der(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
207  double x, int& i, double& y, double& dy, double& ddy) {
208 
209  assert(ytab.dim == xtab.dim) ;
210  assert(dytab.dim == xtab.dim) ;
211 
212  huntm(xtab, x, i) ;
213 
214  // i-- ; // Fortran --> C
215 
216  int i1 = i + 1 ;
217 
218  double dx = xtab(i1) - xtab(i) ;
219 
220  double u = (x - xtab(i)) / dx ;
221  double u2 = u*u ;
222  double u3 = u2*u ;
223 
224  y = ytab(i) * ( 2.*u3 - 3.*u2 + 1.)
225  + ytab(i1) * ( 3.*u2 - 2.*u3)
226  + dytab(i) * dx * ( u3 - 2.*u2 + u )
227  - dytab(i1) * dx * ( u2 - u3 ) ;
228 
229  dy = 6. * ( ytab(i) - ytab(i1) ) * ( u2 - u ) / dx
230  + dytab(i) * ( 3.*u2 - 4.*u + 1. )
231  + dytab(i1) * ( 3.*u2 - 2.*u ) ;
232 
233  ddy = 6 * ( ( ytab(i) - ytab(i1) ) * ( 2.*u - 1. ) / dx
234  + dytab(i) * (6.*u - 4.)
235  + dytab(i1) * (6.*u - 2.) ) / dx ;
236 
237  }
238 
239  //----------------------------------------------
240  // Bi-cubic Hermite interpolation, for 2D arrays
241  //----------------------------------------------
242  void interpol_herm_2d(const Tbl& xtab, const Tbl& ytab, const Tbl& ftab,
243  const Tbl& dfdxtab, const Tbl& dfdytab, const Tbl&
244  d2fdxdytab, double x, double y, double& f, double&
245  dfdx, double& dfdy) {
246 
247  assert(ytab.dim == xtab.dim) ;
248  assert(ftab.dim == xtab.dim) ;
249  assert(dfdxtab.dim == xtab.dim) ;
250  assert(dfdytab.dim == xtab.dim) ;
251  assert(d2fdxdytab.dim == xtab.dim) ;
252 
253  int nbp1, nbp2;
254  nbp1 = xtab.get_dim(0);
255  nbp2 = xtab.get_dim(1);
256 
257  int i_near = 0 ;
258  int j_near = 0 ;
259 
260  while ((xtab(i_near,0) <= x) && (nbp2 > i_near)) {
261  i_near++;
262  }
263  if (i_near != 0) {
264  i_near-- ;
265  }
266  j_near = 0;
267  while ((ytab(i_near,j_near) < y) && (nbp1 > j_near)) {
268  j_near++ ;
269  }
270  if (j_near != 0) {
271  j_near-- ;
272  }
273 
274  int i1 = i_near+1 ; int j1 = j_near+1 ;
275 
276  double dx = xtab(i1, j_near) - xtab(i_near, j_near) ;
277  double dy = ytab(i_near, j1) - ytab(i_near, j_near) ;
278 
279  double u = (x - xtab(i_near, j_near)) / dx ;
280  double v = (y - ytab(i_near, j_near)) / dy ;
281 
282  double u2 = u*u ; double v2 = v*v ;
283  double u3 = u2*u ; double v3 = v2*v ;
284 
285  double psi0_u = 2.*u3 - 3.*u2 + 1. ;
286  double psi0_1mu = -2.*u3 + 3.*u2 ;
287  double psi1_u = u3 - 2.*u2 + u ;
288  double psi1_1mu = -u3 + u2 ;
289 
290  double psi0_v = 2.*v3 - 3.*v2 + 1. ;
291  double psi0_1mv = -2.*v3 + 3.*v2 ;
292  double psi1_v = v3 - 2.*v2 + v ;
293  double psi1_1mv = -v3 + v2 ;
294 
295  f = ftab(i_near, j_near) * psi0_u * psi0_v
296  + ftab(i1, j_near) * psi0_1mu * psi0_v
297  + ftab(i_near, j1) * psi0_u * psi0_1mv
298  + ftab(i1, j1) * psi0_1mu * psi0_1mv ;
299 
300  f += (dfdxtab(i_near, j_near) * psi1_u * psi0_v
301  - dfdxtab(i1, j_near) * psi1_1mu * psi0_v
302  + dfdxtab(i_near, j1) * psi1_u * psi0_1mv
303  - dfdxtab(i1, j1) * psi1_1mu * psi0_1mv) * dx ;
304 
305  f += (dfdytab(i_near, j_near) * psi0_u * psi1_v
306  + dfdytab(i1, j_near) * psi0_1mu * psi1_v
307  - dfdytab(i_near, j1) * psi0_u * psi1_1mv
308  - dfdytab(i1, j1) * psi0_1mu * psi1_1mv) * dy ;
309 
310  f += (d2fdxdytab(i_near, j_near) * psi1_u * psi1_v
311  - d2fdxdytab(i1, j_near) * psi1_1mu * psi1_v
312  - d2fdxdytab(i_near, j1) * psi1_u * psi1_1mv
313  + d2fdxdytab(i1, j1) * psi1_1mu * psi1_1mv) * dx * dy ;
314 
315  double dpsi0_u = 6.*(u2 - u) ;
316  double dpsi0_1mu = 6.*(u2 - u) ;
317  double dpsi1_u = 3.*u2 - 4.*u + 1. ;
318  double dpsi1_1mu = 3.*u2 - 2.*u ;
319 
320  dfdx = (ftab(i_near, j_near) * dpsi0_u * psi0_v
321  - ftab(i1, j_near) * dpsi0_1mu * psi0_v
322  + ftab(i_near, j1) * dpsi0_u * psi0_1mv
323  - ftab(i1, j1) * dpsi0_1mu * psi0_1mv ) / dx;
324 
325  dfdx += (dfdxtab(i_near, j_near) * dpsi1_u * psi0_v
326  + dfdxtab(i1, j_near) * dpsi1_1mu * psi0_v
327  + dfdxtab(i_near, j1) * dpsi1_u * psi0_1mv
328  + dfdxtab(i1, j1) * dpsi1_1mu * psi0_1mv) ;
329 
330  dfdx += (dfdytab(i_near, j_near) * dpsi0_u * psi1_v
331  - dfdytab(i1, j_near) * dpsi0_1mu * psi1_v
332  - dfdytab(i_near, j1) * dpsi0_u * psi1_1mv
333  + dfdytab(i1, j1) * dpsi0_1mu * psi1_1mv) * dy /dx ;
334 
335  dfdx += (d2fdxdytab(i_near, j_near) * dpsi1_u * psi1_v
336  + d2fdxdytab(i1, j_near) * dpsi1_1mu * psi1_v
337  - d2fdxdytab(i_near, j1) * dpsi1_u * psi1_1mv
338  - d2fdxdytab(i1, j1) * dpsi1_1mu * psi1_1mv) * dy ;
339 
340  double dpsi0_v = 6.*(v2 - v) ;
341  double dpsi0_1mv = 6.*(v2 - v) ;
342  double dpsi1_v = 3.*v2 - 4.*v + 1. ;
343  double dpsi1_1mv = 3.*v2 - 2.*v ;
344 
345  dfdy = (ftab(i_near, j_near) * psi0_u * dpsi0_v
346  + ftab(i1, j_near) * psi0_1mu * dpsi0_v
347  - ftab(i_near, j1) * psi0_u * dpsi0_1mv
348  - ftab(i1, j1) * psi0_1mu * dpsi0_1mv) / dy ;
349 
350  dfdy += (dfdxtab(i_near, j_near) * psi1_u * dpsi0_v
351  - dfdxtab(i1, j_near) * psi1_1mu * dpsi0_v
352  - dfdxtab(i_near, j1) * psi1_u * dpsi0_1mv
353  + dfdxtab(i1, j1) * psi1_1mu * dpsi0_1mv) * dx / dy ;
354 
355  dfdy += (dfdytab(i_near, j_near) * psi0_u * dpsi1_v
356  + dfdytab(i1, j_near) * psi0_1mu * dpsi1_v
357  + dfdytab(i_near, j1) * psi0_u * dpsi1_1mv
358  + dfdytab(i1, j1) * psi0_1mu * dpsi1_1mv) ;
359 
360  dfdy += (d2fdxdytab(i_near, j_near) * psi1_u * dpsi1_v
361  - d2fdxdytab(i1, j_near) * psi1_1mu * dpsi1_v
362  + d2fdxdytab(i_near, j1) * psi1_u * dpsi1_1mv
363  - d2fdxdytab(i1, j1) * psi1_1mu * dpsi1_1mv) * dx ;
364 
365  return ;
366  }
367 
368 
369 
370  void interpol_herm_2d_sans(const Tbl& xtab, const Tbl& ytab, const Tbl& ftab,
371  const Tbl& dfdxtab, const Tbl& dfdytab, double x,
372  double y, double& f, double& dfdx, double& dfdy) {
373 
374  assert(ytab.dim == xtab.dim) ;
375  assert(ftab.dim == xtab.dim) ;
376  assert(dfdxtab.dim == xtab.dim) ;
377  assert(dfdytab.dim == xtab.dim) ;
378 
379  int nbp1, nbp2;
380  nbp1 = xtab.get_dim(0);
381  nbp2 = xtab.get_dim(1);
382 
383  int i_near = 0 ;
384  int j_near = 0 ;
385 
386  while ((xtab(i_near,0) <= x) && (nbp2 > i_near)) {
387  i_near++;
388  }
389  if (i_near != 0) {
390  i_near-- ;
391  }
392  j_near = 0;
393  while ((ytab(i_near,j_near) < y) && (nbp1 > j_near)) {
394  j_near++ ;
395  }
396  if (j_near != 0) {
397  j_near-- ;
398  }
399 
400  int i1 = i_near+1 ; int j1 = j_near+1 ;
401 
402  double dx = xtab(i1, j_near) - xtab(i_near, j_near) ;
403  double dy = ytab(i_near, j1) - ytab(i_near, j_near) ;
404 
405  double u = (x - xtab(i_near, j_near)) / dx ;
406  double v = (y - ytab(i_near, j_near)) / dy ;
407 
408  double u2 = u*u ; double v2 = v*v ;
409  double u3 = u2*u ; double v3 = v2*v ;
410 
411  double psi0_u = 2.*u3 - 3.*u2 + 1. ;
412  double psi0_1mu = -2.*u3 + 3.*u2 ;
413  double psi1_u = u3 - 2.*u2 + u ;
414  double psi1_1mu = -u3 + u2 ;
415 
416  double psi0_v = 2.*v3 - 3.*v2 + 1. ;
417  double psi0_1mv = -2.*v3 + 3.*v2 ;
418  double psi1_v = v3 - 2.*v2 + v ;
419  double psi1_1mv = -v3 + v2 ;
420 
421  f = ftab(i_near, j_near) * psi0_u * psi0_v
422  + ftab(i1, j_near) * psi0_1mu * psi0_v
423  + ftab(i_near, j1) * psi0_u * psi0_1mv
424  + ftab(i1, j1) * psi0_1mu * psi0_1mv ;
425 
426  f += (dfdxtab(i_near, j_near) * psi1_u * psi0_v
427  - dfdxtab(i1, j_near) * psi1_1mu * psi0_v
428  + dfdxtab(i_near, j1) * psi1_u * psi0_1mv
429  - dfdxtab(i1, j1) * psi1_1mu * psi0_1mv) * dx ;
430 
431  f += (dfdytab(i_near, j_near) * psi0_u * psi1_v
432  + dfdytab(i1, j_near) * psi0_1mu * psi1_v
433  - dfdytab(i_near, j1) * psi0_u * psi1_1mv
434  - dfdytab(i1, j1) * psi0_1mu * psi1_1mv) * dy ;
435 
436  double dpsi0_u = 6.*(u2 - u) ;
437  double dpsi0_1mu = 6.*(u2 - u) ;
438  double dpsi1_u = 3.*u2 - 4.*u + 1. ;
439  double dpsi1_1mu = 3.*u2 - 2.*u ;
440 
441  dfdx = (ftab(i_near, j_near) * dpsi0_u * psi0_v
442  - ftab(i1, j_near) * dpsi0_1mu * psi0_v
443  + ftab(i_near, j1) * dpsi0_u * psi0_1mv
444  - ftab(i1, j1) * dpsi0_1mu * psi0_1mv ) / dx;
445 
446  dfdx += (dfdxtab(i_near, j_near) * dpsi1_u * psi0_v
447  + dfdxtab(i1, j_near) * dpsi1_1mu * psi0_v
448  + dfdxtab(i_near, j1) * dpsi1_u * psi0_1mv
449  + dfdxtab(i1, j1) * dpsi1_1mu * psi0_1mv) ;
450 
451  dfdx += (dfdytab(i_near, j_near) * dpsi0_u * psi1_v
452  - dfdytab(i1, j_near) * dpsi0_1mu * psi1_v
453  - dfdytab(i_near, j1) * dpsi0_u * psi1_1mv
454  + dfdytab(i1, j1) * dpsi0_1mu * psi1_1mv) * dy /dx ;
455 
456  double dpsi0_v = 6.*(v2 - v) ;
457  double dpsi0_1mv = 6.*(v2 - v) ;
458  double dpsi1_v = 3.*v2 - 4.*v + 1. ;
459  double dpsi1_1mv = 3.*v2 - 2.*v ;
460 
461  dfdy = (ftab(i_near, j_near) * psi0_u * dpsi0_v
462  + ftab(i1, j_near) * psi0_1mu * dpsi0_v
463  - ftab(i_near, j1) * psi0_u * dpsi0_1mv
464  - ftab(i1, j1) * psi0_1mu * dpsi0_1mv) / dy ;
465 
466  dfdy += (dfdxtab(i_near, j_near) * psi1_u * dpsi0_v
467  - dfdxtab(i1, j_near) * psi1_1mu * dpsi0_v
468  - dfdxtab(i_near, j1) * psi1_u * dpsi0_1mv
469  + dfdxtab(i1, j1) * psi1_1mu * dpsi0_1mv) * dx / dy ;
470 
471  dfdy += (dfdytab(i_near, j_near) * psi0_u * dpsi1_v
472  + dfdytab(i1, j_near) * psi0_1mu * dpsi1_v
473  + dfdytab(i_near, j1) * psi0_u * dpsi1_1mv
474  + dfdytab(i1, j1) * psi0_1mu * dpsi1_1mv) ;
475 
476  return ;
477  }
478 
479  //--------------------------------------------------------------------
480  // Quintic Hermite interpolation using data from the second derivative
481  //--------------------------------------------------------------------
482  void interpol_herm_2nd_der(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
483  const Tbl& d2ytab, double x, int& i, double& y,
484  double& dy) {
485 
486  assert(ytab.dim == xtab.dim) ;
487  assert(dytab.dim == xtab.dim) ;
488  assert(d2ytab.dim == xtab.dim) ;
489 
490  huntm(xtab, x, i) ;
491 
492  int i1 = i + 1 ;
493 
494  double dx = xtab(i1) - xtab(i) ;
495 
496  double u = (x - xtab(i)) / dx ;
497  double u2 = u*u ;
498  double u3 = u2*u ;
499  double u4 = u2*u2 ;
500  double u5 = u3*u2 ;
501 
502  double v = 1. - u ;
503  double v2 = v*v ;
504  double v3 = v2*v ;
505  double v4 = v2*v2 ;
506  double v5 = v3*v2 ;
507 
508  y = ytab(i) * ( -6.*u5 + 15.*u4 - 10.*u3 + 1. )
509  + ytab(i1) * ( -6.*v5 + 15.*v4 - 10.*v3 + 1. )
510  + dytab(i) * dx * ( -3.*u5 + 8.*u4 -6.*u3 + u )
511  - dytab(i1) * dx * ( -3.*v5 + 8.*v4 -6.*v3 + v )
512  + d2ytab(i) * dx*dx * ( -0.5*u5 + 1.5*u4 - 1.5*u3 + 0.5*u2 )
513  + d2ytab(i1) * dx*dx * ( -0.5*v5 + 1.5*v4 - 1.5*v3 + 0.5*v2 ) ;
514 
515  dy = 30.*( ytab(i) / dx * ( -u4 + 2.*u3 - u2 )
516  - ytab(i1) / dx * ( -v4 + 2.*v3 - v2 ) )
517  + dytab(i) * ( -15.*u4 + 32.*u3 - 18.*u2 + 1. )
518  + dytab(i1) * ( -15.*v4 + 32.*v3 - 18.*v2 + 1. )
519  + d2ytab(i) * dx * ( -2.5*u4 + 6.*u3 -4.5*u2 + u )
520  - d2ytab(i1) * dx * ( -2.5*v4 + 6.*v3 -4.5*v2 + v ) ;
521  }
522 
523 } // End of namespace Lorene
524 
Lorene
Lorene prototypes.
Definition: app_hor.h:64