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
28char 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
88namespace 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 prototypes.
Definition app_hor.h:64