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 $" ;
93 void huntm(
const Tbl& xx,
double& x,
int& i_low) {
95 assert (xx.get_etat() == ETATQCQ) ;
96 int nx = xx.get_taille() ;
97 bool ascend = ( xx(nx-1) > xx(0) ) ;
99 if ( (i_low < 0)||(i_low>=nx) ) {
105 if ( (x >= xx(i_low)) == ascend ) {
106 if (i_low == nx -1) return ;
108 while ( (x >= xx(i_hi)) == ascend ) {
123 while ( (x < xx(i_low)) == ascend ) {
130 else i_low = i_hi - inc ;
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 ;
139 if (x == xx(nx-1)) i_low = nx-2 ;
140 if (x == xx(0)) i_low = 0 ;
147 void interpol_linear(
const Tbl& xtab,
const Tbl& ytab,
148 double x,
int& i,
double& y) {
150 assert(ytab.dim == xtab.dim) ;
158 double y1 = ytab(i) ;
159 double y2 = ytab(i1) ;
161 double x1 = xtab(i) ;
162 double x2 = xtab(i1) ;
165 double a = (y1-y2)/x12 ;
166 double b = (x1*y2-y1*x2)/x12 ;
175 void interpol_herm(
const Tbl& xtab,
const Tbl& ytab,
const Tbl& dytab,
176 double x,
int& i,
double& y,
double& dy) {
178 assert(ytab.dim == xtab.dim) ;
179 assert(dytab.dim == xtab.dim) ;
185 double dx = xtab(i1) - xtab(i) ;
187 double u = (x - xtab(i)) / dx ;
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 ) ;
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 ) ;
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) {
209 assert(ytab.dim == xtab.dim) ;
210 assert(dytab.dim == xtab.dim) ;
218 double dx = xtab(i1) - xtab(i) ;
220 double u = (x - xtab(i)) / dx ;
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 ) ;
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 ) ;
233 ddy = 6 * ( ( ytab(i) - ytab(i1) ) * ( 2.*u - 1. ) / dx
234 + dytab(i) * (6.*u - 4.)
235 + dytab(i1) * (6.*u - 2.) ) / dx ;
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) {
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) ;
254 nbp1 = xtab.get_dim(0);
255 nbp2 = xtab.get_dim(1);
260 while ((xtab(i_near,0) <= x) && (nbp2 > i_near)) {
267 while ((ytab(i_near,j_near) < y) && (nbp1 > j_near)) {
274 int i1 = i_near+1 ;
int j1 = j_near+1 ;
276 double dx = xtab(i1, j_near) - xtab(i_near, j_near) ;
277 double dy = ytab(i_near, j1) - ytab(i_near, j_near) ;
279 double u = (x - xtab(i_near, j_near)) / dx ;
280 double v = (y - ytab(i_near, j_near)) / dy ;
282 double u2 = u*u ;
double v2 = v*v ;
283 double u3 = u2*u ;
double v3 = v2*v ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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 ;
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) ;
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 ;
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) {
374 assert(ytab.dim == xtab.dim) ;
375 assert(ftab.dim == xtab.dim) ;
376 assert(dfdxtab.dim == xtab.dim) ;
377 assert(dfdytab.dim == xtab.dim) ;
380 nbp1 = xtab.get_dim(0);
381 nbp2 = xtab.get_dim(1);
386 while ((xtab(i_near,0) <= x) && (nbp2 > i_near)) {
393 while ((ytab(i_near,j_near) < y) && (nbp1 > j_near)) {
400 int i1 = i_near+1 ;
int j1 = j_near+1 ;
402 double dx = xtab(i1, j_near) - xtab(i_near, j_near) ;
403 double dy = ytab(i_near, j1) - ytab(i_near, j_near) ;
405 double u = (x - xtab(i_near, j_near)) / dx ;
406 double v = (y - ytab(i_near, j_near)) / dy ;
408 double u2 = u*u ;
double v2 = v*v ;
409 double u3 = u2*u ;
double v3 = v2*v ;
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 ;
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 ;
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 ;
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 ;
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 ;
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 ;
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;
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) ;
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 ;
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 ;
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 ;
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 ;
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) ;
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,
486 assert(ytab.dim == xtab.dim) ;
487 assert(dytab.dim == xtab.dim) ;
488 assert(d2ytab.dim == xtab.dim) ;
494 double dx = xtab(i1) - xtab(i) ;
496 double u = (x - xtab(i)) / dx ;
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 ) ;
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 ) ;