1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21!/*---------------------------------------------------------------------------- 22! This module contains a data type (spline_t) to contain 1D functions, 23! along with a series of procedures to manage them (to define them, to obtain 24! its values, to operate on them, etc). The internal representation of the 25! functions is done through cubic splines, handled by the GSL library. For 26! the user of the module, this internal representation is hidden; one just 27! works with what are called hereafter "spline functions". 28! 29! To define a function, one must supply a set {x(i),y(i)} of pairs of values 30! -- the abscissa and the value of the function. 31! 32! [*] DATA TYPE: 33! To define a spline function: 34! type(spline_t) :: f 35! 36! [1] INITIALIZATION: 37! Before using any function, one should initialize it: 38! 39! Interface: 40! subroutine spline_init(spl) 41! type(spline_t), intent(out) :: spl [or spl(:) or spl(:, :)] 42! end subroutine spline_init 43! 44! Usage: 45! call spline_init(f) 46! 47! [2] FINALIZATION: 48! To empty any allocated space, one should finalize the function: 49! 50! Interface: 51! subroutine spline_end(spl) 52! type(spline_t), intent(inout) :: spl [or spl(:) or spl(:, :)] 53! end subroutine spline_end 54! 55! Usage 56! type(spline_t) :: f 57! call spline_end(f) 58! 59! [3] TO DEFINE A FUNCTION: 60! To "fill" an initialized function f,use spline_fit 61! 62! Interface: 63! subroutine spline_fit(n, x, y, spl) 64! integer, intent(in) :: nrc 65! real(X), intent(in) :: x(n), y(n) 66! type(spline_t), intent(out) :: spl 67! end subroutine spline_fit 68! 69! (X may be 4 or eight, for single or double precision) 70! 71! Usage: 72! call spline_fit(n, x, y, f) 73! n is the number of values that are supplied, x the abscissas, and y 74! the value of the function to represent at each point. 75! 76! [4] FUNCTION VALUES: 77! To retrieve the value of a function at a given point: 78! 79! Interface: 80! 81! real(8) function spline_eval(spl, x) 82! type(spline_t), intent(in) :: spl 83! real(8), intent(in) :: x 84! end function spline_eval 85! 86! Usage: 87! To retrieve the value of function f at point x, and place it into 88! real value val: 89! val = spline_eval(f, x) 90! 91! [5] SUM: 92! If you have two defined functions, f and g, and an initialized function 93! h, you may sum f and g and place the result onto g. For this purpose, 94! the grid that defined the first operand, f, will be used -- the values 95! of g will be interpolated to get the sum and place it in h. 96! 97! Interface: 98! subroutine spline_sum(spl1, spl2, splsum) 99! type(spline_t), intent(in) :: spl1, spl2 100! type(spline_t), intent(out) :: splsum 101! end subroutine spline_sum 102! 103! Usage: 104! call spline_init(f) 105! call spline_init(g) 106! call spline_init(h) 107! call spline_fit(npointsf, xf, yf, f) 108! call spline_fit(npointsg, xg, yg, g) 109! call spline_sum(f, g, h) 110! 111! [6] MULTIPLICATION BY A SCALAR 112! You may multiply a given spline-represented spline by a real number: 113! 114! Interface: 115! subroutine spline_times(a, spl) 116! type(spline_t), intent(inout) :: spl 117! real(8), intent(in) :: a 118! end subroutine spline_times 119! 120! Usage: 121! call spline_init(f) 122! call spline_fit(npoints, x, y, f) ! Fill f with y values at x points 123! call spline_times(a, f) ! Now f contains a*y values at x points. 124! 125! [7] INTEGRAL: 126! Given a defined function, the function spline_integral returns its 127! integral. The interval of integration may or may not be supplied. 128! 129! Interface: 130! real(8) function spline_integral(spl [,a,b]) 131! type(spline_integral), intent(in) :: spl 132! real(8), intent(in), optional :: a, b 133! end function spline_integral 134! 135! [8] DOT PRODUCT: 136! Given two defined functions f and g, the function spline_dotp returns 137! the value of their dot-product: int {dx f(x)*g(x)}. The mesh used to do 138! so the mesh of the fist-function (note that as a result the definition 139! is no longer conmutative). 140! 141! Interface: 142! real(8) function spline_dotp(spl1, spl2) 143! type(spline_t), intent(in) :: spl1, spl2 144! end function spline_dotp 145! 146! Note: The following routines, spline_3dft, spline_cut and spline_filter, 147! assume that the spline functions are the radial part of a 3 dimensional function with 148! spherical symmetry. This is why the Fourier transform of F(\vec{r}) = f(r), is: 149! F(\vec{g}) = f(g) = \frac{4\pi}{g} \int_{0}^{\infty} { r*sin(g*r)*f(r) } 150! which coincides with the inverse Fourier transform, except that the inverse Fourier 151! transform should be multiplied by a (2*\pi)^{-3} factor. 152! 153! [9] FOURIER TRANSFORM: 154! If a spline function f is representing the radial part of a spherically 155! symmetric function F(\vec{r}), its Fourier transform is: 156! F(\vec{g}) = f(g) = \frac{4\pi}{g} \int_{0}^{\infty} { r*sin(g*r)*f(r) } 157! It is assumed that f is defined in some interval [0,rmax], and that it is 158! null at rmax and beyond. One may obtain f(g) by using spline_3dft. 159! The result is placed on the spline data type splw. This has to be initialized, 160! and may or may not be filled. If it is filled, the abscissas that it has 161! are used to define the function. If it is not filled, an equally spaced 162! grid is constructed to define the function, in the interval [0, gmax], where 163! gmax has to be supplied by the caller. 164! 165! Interface: 166! subroutine spline_3dft(spl, splw, gmax) 167! type(spline_t), intent(in) :: spl 168! type(spline_t), intent(inout) :: splw 169! real(8), intent(in), optional :: gmax 170! end subroutine spline_3dft 171! 172! [10] BESSEL TRANSFORM: 173! 174! 175! [11] CUTTING A FUNCTION: 176! spline_cut multiplies a given function by a cutoff-function, which 177! is defined to be one in [0, cutoff], and \exp\{-beta*(x/cutoff-1)^2\} 178! 179! Interface: 180! subroutine spline_cut(spl, cutoff, beta) 181! type(spline_t), intent(in) :: spl 182! real(8), intent(in) :: cutoff, beta 183! end subroutine spline_cut 184! 185! [13] PRINTING A FUNCTION: 186! It prints to a file the (x,y) values that were used to define a function. 187! The file is pointed to by its Fortran unit given by argument iunit. 188! 189! Interface: 190! subroutine spline_print(spl, iunit) 191! type(spline_t), intent(in) :: spl [ or spl(:) or spl(:, :)] 192! integer, intent(in) :: iunit 193! end subroutine spline_print 194! 195! [14] DIFFERENTIATE A FUNCTION: 196! 197!----------------------------------------------------------------------------*/! 198module splines_oct_m 199 use global_oct_m 200 use iso_c_binding 201 use loct_math_oct_m 202 use messages_oct_m 203 use parser_oct_m 204 use profiling_oct_m 205 206 implicit none 207 208 ! Define which routines can be seen from the outside. 209 private 210 public :: & 211 spline_t, & ! [*] 212 spline_init, & ! [1] 213 spline_end, & ! [2] 214 spline_fit, & ! [3] 215 spline_eval, & ! [4] 216 spline_eval_vec, & ! [4] 217 spline_sum, & ! [5] 218 spline_times, & ! [6] 219 spline_integral, & ! [7] 220 spline_dotp, & ! [8] 221 spline_3dft, & ! [9] 222 spline_besselft, & ! [10] 223 spline_cut, & ! [11] 224 spline_print, & ! [13] 225 spline_der, & 226 spline_der2, & 227 spline_div, & 228 spline_mult, & 229 spline_force_pos, & 230 spline_range_min, & 231 spline_range_max, & 232 spline_cutoff_radius 233 234 !> the basic spline datatype 235 type spline_t 236 private 237 real(8) :: x_limit(2) 238 type(c_ptr) :: spl, acc 239 end type spline_t 240 241 !> Both the filling of the function, and the retrieval of the values 242 !! may be done using single- or double-precision values. 243 interface spline_eval_vec 244 module procedure spline_eval8_array 245 module procedure spline_evalz_array 246 end interface spline_eval_vec 247 248 !> The integral may be done with or without integration limits, but 249 !! we want the interface to be common. 250 interface spline_integral 251 module procedure spline_integral_full 252 module procedure spline_integral_limits 253 end interface spline_integral 254 255 !> Some operations may be done for one spline-function, or for an array of them 256 interface spline_init 257 module procedure spline_init_0 258 module procedure spline_init_1 259 module procedure spline_init_2 260 end interface spline_init 261 262 interface spline_end 263 module procedure spline_end_0 264 module procedure spline_end_1 265 module procedure spline_end_2 266 end interface spline_end 267 268 interface spline_print 269 module procedure spline_print_0 270 module procedure spline_print_1 271 module procedure spline_print_2 272 end interface spline_print 273 274 ! These are interfaces to functions defined in oct_gsl_f.c, which in turn 275 ! take care of calling the GSL library. 276 interface 277 subroutine oct_spline_end(spl, acc) 278 use iso_c_binding 279 implicit none 280 type(c_ptr), intent(inout) :: spl, acc 281 end subroutine oct_spline_end 282 283 subroutine oct_spline_fit(nrc, x, y, spl, acc) 284 use iso_c_binding 285 implicit none 286 integer, intent(in) :: nrc 287 real(8), intent(in) :: x 288 real(8), intent(in) :: y 289 type(c_ptr), intent(inout) :: spl 290 type(c_ptr), intent(inout) :: acc 291 end subroutine oct_spline_fit 292 293 real(8) pure function oct_spline_eval(x, spl, acc) 294 use iso_c_binding 295 296 real(8), intent(in) :: x 297 type(c_ptr), intent(in) :: spl 298 type(c_ptr), intent(in) :: acc 299 end function oct_spline_eval 300 301 pure subroutine oct_spline_eval_array(nn, xf, spl, acc) 302 use iso_c_binding 303 304 integer, intent(in) :: nn 305 real(8), intent(inout) :: xf 306 type(c_ptr), intent(in) :: spl 307 type(c_ptr), intent(in) :: acc 308 end subroutine oct_spline_eval_array 309 310 pure subroutine oct_spline_eval_arrayz(nn, xf, spl, acc) 311 use iso_c_binding 312 313 integer, intent(in) :: nn 314 complex(8), intent(inout) :: xf 315 type(c_ptr), intent(in) :: spl 316 type(c_ptr), intent(in) :: acc 317 end subroutine oct_spline_eval_arrayz 318 319 real(8) pure function oct_spline_eval_der(x, spl, acc) 320 use iso_c_binding 321 322 real(8), intent(in) :: x 323 type(c_ptr), intent(in) :: spl 324 type(c_ptr), intent(in) :: acc 325 end function oct_spline_eval_der 326 327 real(8) pure function oct_spline_eval_der2(x, spl, acc) 328 use iso_c_binding 329 330 real(8), intent(in) :: x 331 type(c_ptr), intent(in) :: spl 332 type(c_ptr), intent(in) :: acc 333 end function oct_spline_eval_der2 334 335 integer pure function oct_spline_npoints(spl, acc) 336 use iso_c_binding 337 338 type(c_ptr), intent(in) :: spl 339 type(c_ptr), intent(in) :: acc 340 end function oct_spline_npoints 341 342 pure subroutine oct_spline_x(spl, acc, x) 343 use iso_c_binding 344 345 type(c_ptr), intent(in) :: spl 346 type(c_ptr), intent(in) :: acc 347 real(8), intent(out) :: x 348 end subroutine oct_spline_x 349 350 subroutine oct_spline_y(spl, acc, y) 351 use iso_c_binding 352 implicit none 353 type(c_ptr), intent(in) :: spl 354 type(c_ptr), intent(in) :: acc 355 real(8), intent(out) :: y 356 end subroutine oct_spline_y 357 358 real(8) pure function oct_spline_eval_integ(spl, a, b, acc) 359 use iso_c_binding 360 implicit none 361 type(c_ptr), intent(in) :: spl 362 real(8), intent(in) :: a 363 real(8), intent(in) :: b 364 type(c_ptr), intent(in) :: acc 365 end function oct_spline_eval_integ 366 367 real(8) pure function oct_spline_eval_integ_full(spl, acc) 368 use iso_c_binding 369 implicit none 370 type(c_ptr), intent(in) :: spl 371 type(c_ptr), intent(in) :: acc 372 end function oct_spline_eval_integ_full 373 end interface 374 375contains 376 377 !------------------------------------------------------------ 378 subroutine spline_init_0(spl) 379 type(spline_t), intent(out) :: spl 380 381 ! No PUSH SUB, called too often 382 383 spl%spl = c_null_ptr 384 spl%acc = c_null_ptr 385 386 ! deliberately illegal values, for checking 387 spl%x_limit(1) = -1d0 388 spl%x_limit(2) = -2d0 389 390 end subroutine spline_init_0 391 392 393 !------------------------------------------------------------ 394 subroutine spline_init_1(spl) 395 type(spline_t), intent(out) :: spl(:) 396 397 integer :: i 398 399 PUSH_SUB(spline_init_1) 400 401 do i = 1, size(spl) 402 call spline_init_0(spl(i)) 403 end do 404 405 POP_SUB(spline_init_1) 406 end subroutine spline_init_1 407 408 409 !------------------------------------------------------------ 410 subroutine spline_init_2(spl) 411 type(spline_t), intent(out) :: spl(:, :) 412 413 integer :: i, j 414 415 PUSH_SUB(spline_init_2) 416 417 do i = 1, size(spl, 1) 418 do j = 1, size(spl, 2) 419 call spline_init_0(spl(i, j)) 420 end do 421 end do 422 423 POP_SUB(spline_init_2) 424 end subroutine spline_init_2 425 426 427 !------------------------------------------------------------ 428 subroutine spline_end_0(spl) 429 type(spline_t), intent(inout) :: spl 430 431 ! No PUSH SUB, called too often. 432 433 if(c_associated(spl%spl) .and. c_associated(spl%acc)) then 434 call oct_spline_end(spl%spl, spl%acc) 435 end if 436 spl%spl = c_null_ptr 437 spl%acc = c_null_ptr 438 439 end subroutine spline_end_0 440 441 442 !------------------------------------------------------------ 443 subroutine spline_end_1(spl) 444 type(spline_t), intent(inout) :: spl(:) 445 446 integer :: i 447 448 PUSH_SUB(spline_end_1) 449 450 do i = 1, size(spl) 451 call spline_end_0(spl(i)) 452 end do 453 454 POP_SUB(spline_end_1) 455 end subroutine spline_end_1 456 457 458 !------------------------------------------------------------ 459 subroutine spline_end_2(spl) 460 type(spline_t), intent(inout) :: spl(:, :) 461 462 integer :: i, j 463 464 PUSH_SUB(spline_end_2) 465 466 do i = 1, size(spl, 1) 467 do j = 1, size(spl, 2) 468 call spline_end_0(spl(i, j)) 469 end do 470 end do 471 472 POP_SUB(spline_end_2) 473 end subroutine spline_end_2 474 475 476 !------------------------------------------------------------ 477 subroutine spline_fit(nrc, rofi, ffit, spl) 478 integer, intent(in) :: nrc 479 real(8), intent(in) :: rofi(:) 480 real(8), intent(in) :: ffit(:) 481 type(spline_t), intent(inout) :: spl 482 483 !No PUSH SUB, called too often 484 485 spl%x_limit(1) = rofi(1) 486 spl%x_limit(2) = rofi(nrc) 487 call oct_spline_fit(nrc, rofi(1), ffit(1), spl%spl, spl%acc) 488 489 end subroutine spline_fit 490 491 !------------------------------------------------------------ 492 real(8) pure function spline_eval(spl, x) 493 type(spline_t), intent(in) :: spl 494 real(8), intent(in) :: x 495 496 spline_eval = oct_spline_eval(x, spl%spl, spl%acc) 497 end function spline_eval 498 499 500 !------------------------------------------------------------ 501 pure subroutine spline_eval8_array(spl, nn, xf) 502 type(spline_t), intent(in) :: spl 503 integer, intent(in) :: nn 504 real(8), intent(inout) :: xf(:) 505 506 call oct_spline_eval_array(nn, xf(1), spl%spl, spl%acc) 507 end subroutine spline_eval8_array 508 509 510 !------------------------------------------------------------ 511 pure subroutine spline_evalz_array(spl, nn, xf) 512 type(spline_t), intent(in) :: spl 513 integer, intent(in) :: nn 514 complex(8), intent(inout) :: xf(:) 515 516 call oct_spline_eval_arrayz(nn, xf(1), spl%spl, spl%acc) 517 end subroutine spline_evalz_array 518 519 !------------------------------------------------------------ 520 subroutine spline_sum(spl1, spl2, splsum) 521 type(spline_t), intent(in) :: spl1 522 type(spline_t), intent(in) :: spl2 523 type(spline_t), intent(out) :: splsum 524 525 integer :: npoints, i 526 real(8), allocatable :: x(:), y(:), y2(:) 527 528 PUSH_SUB(spline_sum) 529 530 npoints = oct_spline_npoints(spl1%spl, spl1%acc) 531 532 SAFE_ALLOCATE( x(1:npoints)) 533 SAFE_ALLOCATE( y(1:npoints)) 534 SAFE_ALLOCATE(y2(1:npoints)) 535 536 call oct_spline_x(spl1%spl, spl1%acc, x(1)) 537 call oct_spline_y(spl1%spl, spl1%acc, y(1)) 538 539 do i = 1, npoints 540 y2(i) = spline_eval(spl2, x(i)) 541 end do 542 543 y2 = y2 + y 544 call spline_fit(npoints, x, y2, splsum) 545 546 SAFE_DEALLOCATE_A(x) 547 SAFE_DEALLOCATE_A(y) 548 SAFE_DEALLOCATE_A(y2) 549 550 POP_SUB(spline_sum) 551 end subroutine spline_sum 552 553 554 !------------------------------------------------------------ 555 subroutine spline_times(a, spl) 556 FLOAT, intent(in) :: a 557 type(spline_t), intent(inout) :: spl 558 559 integer :: npoints, i 560 real(8), allocatable :: x(:), y(:) 561 562 PUSH_SUB(spline_times) 563 564 npoints = oct_spline_npoints(spl%spl, spl%acc) 565 SAFE_ALLOCATE(x(1:npoints)) 566 SAFE_ALLOCATE(y(1:npoints)) 567 568 call oct_spline_x(spl%spl, spl%acc, x(1)) 569 call oct_spline_y(spl%spl, spl%acc, y(1)) 570 call oct_spline_end(spl%spl, spl%acc) 571 do i = 1, npoints 572 y(i) = a*y(i) 573 end do 574 call spline_fit(npoints, x, y, spl) 575 576 SAFE_DEALLOCATE_A(x) 577 SAFE_DEALLOCATE_A(y) 578 579 POP_SUB(spline_times) 580 end subroutine spline_times 581 582 583 !------------------------------------------------------------ 584 real(8) function spline_integral_full(spl) result(res) 585 type(spline_t), intent(in) :: spl 586 587 PUSH_SUB(spline_integral_full) 588 589 res = oct_spline_eval_integ_full(spl%spl, spl%acc) 590 591 POP_SUB(spline_integral_full) 592 end function spline_integral_full 593 594 595 !------------------------------------------------------------ 596 real(8) pure function spline_integral_limits(spl, a, b) result(res) 597 type(spline_t), intent(in) :: spl 598 real(8), intent(in) :: a 599 real(8), intent(in) :: b 600 601 res = oct_spline_eval_integ(spl%spl, a, b, spl%acc) 602 end function spline_integral_limits 603 604 605 !------------------------------------------------------------ 606 real(8) function spline_dotp(spl1, spl2) result (res) 607 type(spline_t), intent(in) :: spl1 608 type(spline_t), intent(in) :: spl2 609 610 type(spline_t) :: aux 611 integer :: npoints, i 612 real(8), allocatable :: x(:), y(:) 613 614 PUSH_SUB(spline_dotp) 615 616 npoints = oct_spline_npoints(spl1%spl, spl1%acc) 617 SAFE_ALLOCATE(x(1:npoints)) 618 SAFE_ALLOCATE(y(1:npoints)) 619 620 call oct_spline_x(spl1%spl, spl1%acc, x(1)) 621 call oct_spline_y(spl1%spl, spl1%acc, y(1)) 622 do i = 1, npoints 623 y(i) = y(i)*oct_spline_eval(x(i), spl2%spl, spl2%acc) 624 end do 625 call spline_init(aux) 626 call spline_fit(npoints, x, y, aux) 627 res = oct_spline_eval_integ(aux%spl, x(1), x(npoints), aux%acc) 628 629 SAFE_DEALLOCATE_A(x) 630 SAFE_DEALLOCATE_A(y) 631 632 POP_SUB(spline_dotp) 633 end function spline_dotp 634 635 636 !------------------------------------------------------------ 637 subroutine spline_3dft(spl, splw, gmax) 638 type(spline_t), intent(in) :: spl 639 type(spline_t), intent(inout) :: splw 640 FLOAT, optional, intent(in) :: gmax 641 642 type(spline_t) :: aux 643 real(8) :: g, dg 644 integer :: np 645 integer :: npoints, i, j 646 real(8), allocatable :: x(:), y(:), y2(:), xw(:), yw(:) 647 648 PUSH_SUB(spline_3dft) 649 650 npoints = oct_spline_npoints(spl%spl, spl%acc) 651 SAFE_ALLOCATE( x(1:npoints)) 652 SAFE_ALLOCATE( y(1:npoints)) 653 SAFE_ALLOCATE(y2(1:npoints)) 654 655 call oct_spline_x(spl%spl, spl%acc, x(1)) 656 call oct_spline_y(spl%spl, spl%acc, y(1)) 657 658 ! Check whether splw comes with a defined grid, or else build it. 659 if(c_associated(splw%spl)) then 660 np = oct_spline_npoints(splw%spl, splw%acc) 661 SAFE_ALLOCATE(xw(1:np)) 662 SAFE_ALLOCATE(yw(1:np)) 663 call oct_spline_x(splw%spl, splw%acc, xw(1)) 664 ! But now we need to kill the input: 665 call spline_end(splw) 666 else 667 np = 200 ! hard coded value 668 dg = gmax/(np-1) 669 SAFE_ALLOCATE(xw(1:np)) 670 SAFE_ALLOCATE(yw(1:np)) 671 do i = 1, np 672 g = (i-1)*dg 673 xw(i) = g 674 end do 675 end if 676 677 ! The first point, xw(1) = 0.0 and it has to be treated separately. 678 do j = 1, npoints 679 y2(j) = CNST(4.0)*M_PI*y(j)*x(j)**2 680 end do 681 call spline_init(aux) 682 call spline_fit(npoints, x, y2, aux) 683 yw(1) = oct_spline_eval_integ_full(aux%spl, aux%acc) 684 call spline_end(aux) 685 686 do i = 2, np 687 do j = 1, npoints 688 y2(j) = (CNST(4.0)*M_PI/xw(i))*y(j)*x(j)*sin(xw(i)*x(j)) 689 end do 690 call spline_init(aux) 691 call spline_fit(npoints, x, y2, aux) 692 yw(i) = oct_spline_eval_integ_full(aux%spl, aux%acc) 693 call spline_end(aux) 694 end do 695 696 call spline_init(splw) 697 call spline_fit(np, xw, yw, splw) 698 699 SAFE_DEALLOCATE_A(x) 700 SAFE_DEALLOCATE_A(y) 701 SAFE_DEALLOCATE_A(y2) 702 SAFE_DEALLOCATE_A(xw) 703 SAFE_DEALLOCATE_A(yw) 704 705 POP_SUB(spline_3dft) 706 end subroutine spline_3dft 707 708 709 !------------------------------------------------------------ 710 subroutine spline_besselft(spl, splw, l, gmax) 711 type(spline_t), intent(in) :: spl 712 type(spline_t), intent(inout) :: splw 713 integer, intent(in) :: l 714 FLOAT, optional, intent(in) :: gmax 715 716 type(spline_t) :: aux 717 real(8) :: g, dg 718 integer :: np 719 integer :: npoints, i, j 720 real(8), allocatable :: x(:), y(:), y2(:), xw(:), yw(:) 721 722 PUSH_SUB(spline_besselft) 723 724 npoints = oct_spline_npoints(spl%spl, spl%acc) 725 SAFE_ALLOCATE( x(1:npoints)) 726 SAFE_ALLOCATE( y(1:npoints)) 727 SAFE_ALLOCATE(y2(1:npoints)) 728 729 call oct_spline_x(spl%spl, spl%acc, x(1)) 730 call oct_spline_y(spl%spl, spl%acc, y(1)) 731 732 ! Check whether splw comes with a defined grid, or else build it. 733 if(c_associated(splw%spl)) then 734 np = oct_spline_npoints(splw%spl, splw%acc) 735 SAFE_ALLOCATE(xw(1:np)) 736 SAFE_ALLOCATE(yw(1:np)) 737 call oct_spline_x(splw%spl, splw%acc, xw(1)) 738 ! But now we need to kill the input: 739 call spline_end(splw) 740 else 741 ASSERT(present(gmax)) 742 np = 1000 ! hard coded value 743 dg = gmax/(np-1) 744 SAFE_ALLOCATE(xw(1:np)) 745 SAFE_ALLOCATE(yw(1:np)) 746 do i = 1, np 747 g = real(i-1, 8)*dg 748 xw(i) = g 749 end do 750 end if 751 752 do i = 1, np 753 !$omp parallel do 754 do j = 1, npoints 755 y2(j) = y(j) * x(j)**2 * loct_sph_bessel(l, x(j)*xw(i)) 756 end do 757 call spline_init(aux) 758 call spline_fit(npoints, x, y2, aux) 759 yw(i) = sqrt(CNST(2.0)/M_PI)*oct_spline_eval_integ_full(aux%spl, aux%acc) 760 761 call spline_end(aux) 762 end do 763 764 call spline_init(splw) 765 call spline_fit(np, xw, yw, splw) 766 767 SAFE_DEALLOCATE_A(x) 768 SAFE_DEALLOCATE_A(y) 769 SAFE_DEALLOCATE_A(y2) 770 SAFE_DEALLOCATE_A(xw) 771 SAFE_DEALLOCATE_A(yw) 772 773 POP_SUB(spline_besselft) 774 end subroutine spline_besselft 775 776 777 !------------------------------------------------------------ 778 subroutine spline_cut(spl, cutoff, beta) 779 type(spline_t), intent(inout) :: spl 780 FLOAT, intent(in) :: cutoff 781 FLOAT, intent(in) :: beta 782 783 integer :: npoints, i 784 real(8), allocatable :: x(:), y(:) 785 FLOAT :: exp_arg 786 787 PUSH_SUB(spline_cut) 788 789 npoints = oct_spline_npoints(spl%spl, spl%acc) 790 SAFE_ALLOCATE(x(1:npoints)) 791 SAFE_ALLOCATE(y(1:npoints)) 792 793 call oct_spline_x(spl%spl, spl%acc, x(1)) 794 call oct_spline_y(spl%spl, spl%acc, y(1)) 795 call oct_spline_end(spl%spl, spl%acc) 796 do i = npoints, 1, -1 797 if(x(i)<cutoff) then 798 exit 799 end if 800 801 !To avoid underflows 802 exp_arg = -beta*(x(i)/cutoff - CNST(1.0))**2 803 if( exp_arg > CNST(-100)) then 804 y(i) = y(i) * exp(exp_arg) 805 else 806 y(i) = M_ZERO 807 end if 808 end do 809 call spline_fit(npoints, x, y, spl) 810 811 SAFE_DEALLOCATE_A(x) 812 SAFE_DEALLOCATE_A(y) 813 814 POP_SUB(spline_cut) 815 end subroutine spline_cut 816 817 818 !------------------------------------------------------------ 819 subroutine spline_div(spla, splb) 820 type(spline_t), intent(inout) :: spla 821 type(spline_t), intent(in) :: splb 822 823 integer :: npoints, i 824 real(8), allocatable :: x(:), y(:) 825 real(8) :: aa 826 827 PUSH_SUB(spline_div) 828 829 npoints = oct_spline_npoints(spla%spl, spla%acc) 830 831 SAFE_ALLOCATE(x(1:npoints)) 832 SAFE_ALLOCATE(y(1:npoints)) 833 834 call oct_spline_x(spla%spl, spla%acc, x(1)) 835 call oct_spline_y(spla%spl, spla%acc, y(1)) 836 call oct_spline_end(spla%spl, spla%acc) 837 838 ASSERT(splb%x_limit(2) >= splb%x_limit(1)) 839 840 do i = npoints, 1, -1 841 if(x(i) > splb%x_limit(2)) cycle 842 aa = spline_eval(splb, x(i)) 843 y(i) = y(i)/aa 844 end do 845 846 call spline_fit(npoints, x, y, spla) 847 848 SAFE_DEALLOCATE_A(x) 849 SAFE_DEALLOCATE_A(y) 850 851 POP_SUB(spline_div) 852 end subroutine spline_div 853 854 !------------------------------------------------------------ 855 !Force all the values of the spline to be positive 856 !------------------------------------------------------------ 857 subroutine spline_force_pos(spl) 858 type(spline_t), intent(inout) :: spl 859 860 integer :: npoints, i 861 real(8), allocatable :: x(:), y(:) 862 863 PUSH_SUB(spline_force_pos) 864 865 npoints = oct_spline_npoints(spl%spl, spl%acc) 866 867 SAFE_ALLOCATE(x(1:npoints)) 868 SAFE_ALLOCATE(y(1:npoints)) 869 870 call oct_spline_x(spl%spl, spl%acc, x(1)) 871 call oct_spline_y(spl%spl, spl%acc, y(1)) 872 call oct_spline_end(spl%spl, spl%acc) 873 874 do i = npoints, 1, -1 875 y(i) = max(y(i),M_ZERO) 876 end do 877 878 call spline_fit(npoints, x, y, spl) 879 880 SAFE_DEALLOCATE_A(x) 881 SAFE_DEALLOCATE_A(y) 882 883 POP_SUB(spline_force_pos) 884 end subroutine spline_force_pos 885 886 887 !------------------------------------------------------------ 888 subroutine spline_mult(spla, splb) 889 type(spline_t), intent(inout) :: spla 890 type(spline_t), intent(in) :: splb 891 892 integer :: npoints, i 893 real(8), allocatable :: x(:), y(:) 894 real(8) :: aa 895 896 PUSH_SUB(spline_mult) 897 898 npoints = oct_spline_npoints(spla%spl, spla%acc) 899 900 SAFE_ALLOCATE(x(1:npoints)) 901 SAFE_ALLOCATE(y(1:npoints)) 902 903 call oct_spline_x(spla%spl, spla%acc, x(1)) 904 call oct_spline_y(spla%spl, spla%acc, y(1)) 905 call oct_spline_end(spla%spl, spla%acc) 906 907 ASSERT(splb%x_limit(2) >= splb%x_limit(1)) 908 909 do i = npoints, 1, -1 910 if(x(i) > splb%x_limit(2)) then 911 aa = M_ZERO 912 else 913 aa = spline_eval(splb, x(i)) 914 end if 915 y(i) = y(i)*aa 916 end do 917 918 call spline_fit(npoints, x, y, spla) 919 920 SAFE_DEALLOCATE_A(x) 921 SAFE_DEALLOCATE_A(y) 922 923 POP_SUB(spline_mult) 924 end subroutine spline_mult 925 926 927 !------------------------------------------------------------ 928 subroutine spline_der(spl, dspl) 929 type(spline_t), intent(in) :: spl 930 type(spline_t), intent(inout) :: dspl 931 932 integer :: npoints, i 933 real(8), allocatable :: x(:), y(:) 934 935 PUSH_SUB(spline_der) 936 937 ! Use the grid of dspl if it is present, otherwise use the same one of spl. 938 if(.not. c_associated(dspl%spl)) then ! use the grid of spl 939 npoints = oct_spline_npoints(spl%spl, spl%acc) 940 SAFE_ALLOCATE(x(1:npoints)) 941 SAFE_ALLOCATE(y(1:npoints)) 942 call oct_spline_x(spl%spl, spl%acc, x(1)) 943 else ! use the grid of dspl 944 npoints = oct_spline_npoints(dspl%spl, dspl%acc) 945 SAFE_ALLOCATE(x(1:npoints)) 946 SAFE_ALLOCATE(y(1:npoints)) 947 call oct_spline_x(dspl%spl, dspl%acc, x(1)) 948 end if 949 do i = 1, npoints 950 y(i) = oct_spline_eval_der(x(i), spl%spl, spl%acc) 951 end do 952 call spline_fit(npoints, x, y, dspl) 953 954 SAFE_DEALLOCATE_A(x) 955 SAFE_DEALLOCATE_A(y) 956 957 POP_SUB(spline_der) 958 end subroutine spline_der 959 960 961 !------------------------------------------------------------ 962 subroutine spline_der2(spl, dspl) 963 type(spline_t), intent(in) :: spl 964 type(spline_t), intent(inout) :: dspl 965 966 integer :: npoints, i 967 real(8), allocatable :: x(:), y(:) 968 969 PUSH_SUB(spline_der2) 970 971 ! Use the grid of dspl if it is present, otherwise use the same one of spl. 972 if(.not. c_associated(dspl%spl)) then ! use the grid of spl 973 npoints = oct_spline_npoints(spl%spl, spl%acc) 974 SAFE_ALLOCATE(x(1:npoints)) 975 SAFE_ALLOCATE(y(1:npoints)) 976 call oct_spline_x(spl%spl, spl%acc, x(1)) 977 else ! use the grid of dspl 978 npoints = oct_spline_npoints(dspl%spl, dspl%acc) 979 SAFE_ALLOCATE(x(1:npoints)) 980 SAFE_ALLOCATE(y(1:npoints)) 981 call oct_spline_x(dspl%spl, dspl%acc, x(1)) 982 end if 983 do i = 1, npoints 984 y(i) = oct_spline_eval_der2(x(i), spl%spl, spl%acc) 985 end do 986 call spline_fit(npoints, x, y, dspl) 987 988 SAFE_DEALLOCATE_A(x) 989 SAFE_DEALLOCATE_A(y) 990 991 POP_SUB(spline_der2) 992 end subroutine spline_der2 993 994 995 !------------------------------------------------------------ 996 subroutine spline_print_0(spl, iunit) 997 type(spline_t), intent(in) :: spl 998 integer, intent(in) :: iunit 999 1000 integer :: np, i 1001 real(8), allocatable :: x(:), y(:) 1002 1003 PUSH_SUB(spline_print_0) 1004 1005 np = oct_spline_npoints(spl%spl, spl%acc) 1006 SAFE_ALLOCATE(x(1:np)) 1007 SAFE_ALLOCATE(y(1:np)) 1008 1009 call oct_spline_x(spl%spl, spl%acc, x(1)) 1010 call oct_spline_y(spl%spl, spl%acc, y(1)) 1011 do i = 1, np 1012 write(iunit, '(2es16.8)') x(i), y(i) 1013 end do 1014 1015 SAFE_DEALLOCATE_A(x) 1016 SAFE_DEALLOCATE_A(y) 1017 1018 POP_SUB(spline_print_0) 1019 end subroutine spline_print_0 1020 1021 1022 !------------------------------------------------------------ 1023 subroutine spline_print_1(spl, iunit) 1024 type(spline_t), intent(in) :: spl(:) 1025 integer, intent(in) :: iunit 1026 1027 character(len=4) :: fm 1028 integer :: np, i, n, j 1029 real(8), allocatable :: x(:), y(:) 1030 1031 PUSH_SUB(spline_print_1) 1032 1033 n = size(spl) 1034 if(n<=0) then 1035 POP_SUB(spline_print_1) 1036 return 1037 end if 1038 1039 write(fm,'(i4)') n + 1; fm = adjustl(fm) 1040 np = oct_spline_npoints(spl(1)%spl, spl(1)%acc) 1041 SAFE_ALLOCATE(x(1:np)) 1042 SAFE_ALLOCATE(y(1:np)) 1043 1044 call oct_spline_x(spl(1)%spl, spl(1)%acc, x(1)) 1045 call oct_spline_y(spl(1)%spl, spl(1)%acc, y(1)) 1046 do i = 1, np 1047 write(iunit, '('//trim(fm)//'es16.8)') x(i), (spline_eval(spl(j), x(i)), j = 1, size(spl)) 1048 end do 1049 1050 SAFE_DEALLOCATE_A(x) 1051 SAFE_DEALLOCATE_A(y) 1052 1053 POP_SUB(spline_print_1) 1054 end subroutine spline_print_1 1055 1056 1057 !------------------------------------------------------------ 1058 subroutine spline_print_2(spl, iunit) 1059 type(spline_t), intent(in) :: spl(:, :) 1060 integer, intent(in) :: iunit 1061 1062 character(len=4) :: fm 1063 integer :: np, i, n1, n2, j, k 1064 real(8), allocatable :: x(:), y(:) 1065 1066 PUSH_SUB(spline_print_2) 1067 1068 n1 = size(spl, 1); n2 = size(spl, 2) 1069 if(n1*n2<=0) then 1070 POP_SUB(spline_print_2) 1071 return 1072 end if 1073 1074 write(fm,'(i4)') n1*n2 + 1; fm = adjustl(fm) 1075 np = oct_spline_npoints(spl(1, 1)%spl, spl(1, 1)%acc) 1076 1077 SAFE_ALLOCATE(x(1:np)) 1078 SAFE_ALLOCATE(y(1:np)) 1079 1080 call oct_spline_x(spl(1, 1)%spl, spl(1, 1)%acc, x(1)) 1081 call oct_spline_y(spl(1, 1)%spl, spl(1, 1)%acc, y(1)) 1082 do i = 1, np 1083 write(iunit, '('//trim(fm)//'es16.8)') x(i), & 1084 ((spline_eval(spl(j, k), x(i)), j = 1, n1), k = 1, n2) 1085 end do 1086 1087 SAFE_DEALLOCATE_A(x) 1088 SAFE_DEALLOCATE_A(y) 1089 1090 POP_SUB(spline_print_2) 1091 end subroutine spline_print_2 1092 1093 1094 !------------------------------------------------------------ 1095 FLOAT function spline_cutoff_radius(spl, threshold) result(r) 1096 type(spline_t), intent(in) :: spl 1097 FLOAT, intent(in) :: threshold 1098 1099 integer :: ii, jj 1100 FLOAT, parameter :: dx = CNST(0.01) 1101 1102 ! No PUSH SUB, called too often. 1103 1104 ASSERT(spl%x_limit(2) >= spl%x_limit(1)) 1105 1106 jj = int(spl%x_limit(2)/dx) + 1 1107 1108 do ii = jj, 1, -1 1109 1110 r = dx*(ii-1) 1111 1112 ! The first point might not be inside range, so skip it, this 1113 ! should be done in a better way, but doing it introduces small 1114 ! numerical differences in many tests, so it is a lot of work. 1115 if(r > spl%x_limit(2)) cycle 1116 1117 if(abs(spline_eval(spl, r)) > threshold) exit 1118 end do 1119 1120 end function spline_cutoff_radius 1121 1122 ! ------------------------------------------------------- 1123 1124 FLOAT pure function spline_range_min(this) result(range) 1125 type(spline_t), intent(in) :: this 1126 1127 range = this%x_limit(1) 1128 1129 end function spline_range_min 1130 1131 ! ------------------------------------------------------- 1132 1133 FLOAT pure function spline_range_max(this) result(range) 1134 type(spline_t), intent(in) :: this 1135 1136 range = this%x_limit(2) 1137 1138 end function spline_range_max 1139 1140end module splines_oct_m 1141 1142!! Local Variables: 1143!! mode: f90 1144!! coding: utf-8 1145!! End: 1146