22 #include "./base/base_uses.f90"
27 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'splines_methods'
49 TYPE(spline_data_type),
INTENT(INOUT) :: spl
50 INTEGER,
INTENT(IN) :: nn
54 IF (
ASSOCIATED(spl%y))
THEN
59 IF (
ASSOCIATED(spl%y2))
THEN
64 ALLOCATE (spl%y(1:nn))
65 ALLOCATE (spl%y2(1:nn))
101 TYPE(spline_data_type),
INTENT(INOUT) :: spl
102 REAL(kind=
dp),
INTENT(IN) :: dx
103 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: y1a, y1b
106 REAL(kind=
dp) :: p, s
107 REAL(kind=
dp),
POINTER :: ww(:)
110 spl%xn = spl%x1 + (n - 1)*dx
113 spl%h26 = dx**2/6.0_dp
115 IF (
PRESENT(y1a))
THEN
117 ww(1) = 3.0_dp*((spl%y(2) - spl%y(1))/dx - y1a)/dx
124 p = 0.5_dp*spl%y2(i - 1) + 2.0_dp
125 spl%y2(i) = -0.5_dp/p
126 ww(i) = (3.0_dp*(spl%y(i + 1) - 2.0_dp*spl%y(i) + spl%y(i - 1))/(dx*dx) &
127 - 0.5_dp*ww(i - 1))/p
129 IF (
PRESENT(y1b))
THEN
130 spl%y2(n) = (3.0_dp*(y1b - (spl%y(n) - spl%y(n - 1))/dx)/dx - &
131 0.5_dp*ww(n - 1))/(0.5_dp*spl%y2(n - 1) + 1.0_dp)
136 spl%y2(i) = spl%y2(i)*spl%y2(i + 1) + ww(i)
162 TYPE(spline_data_p_type),
DIMENSION(:),
POINTER :: spl_p
163 REAL(kind=
dp),
INTENT(IN) :: xxi
164 REAL(kind=
dp),
INTENT(OUT) :: y1
165 TYPE(spline_factor_type),
POINTER :: spl_f
166 TYPE(cp_logger_type),
POINTER :: logger
169 REAL(kind=
dp),
PARAMETER :: f13 = 1.0_dp/3.0_dp
171 INTEGER :: i, output_unit
172 REAL(kind=
dp) :: a, b, h26, invh, x4, xx, xx0, y2hi, &
176 xx = spl_f%rscale(1)*xx0
178 h26 = spl_p(1)%spline_data%h26
179 invh = spl_p(1)%spline_data%invh
180 IF (xx >= spl_p(1)%spline_data%xn)
THEN
185 yy = spl_p(1)%spline_data%xn - spl_p(1)%spline_data%h
186 WRITE (output_unit, fmt=
'(/,80("*"),/,"*",1X,"Value of r in Input =",F11.6,'// &
187 '" not in the spline range. Using =",F11.6,T80,"*",/,80("*"))') sqrt(1.0_dp/xx), sqrt(1.0_dp/yy)
190 i = int((xx - spl_p(1)%spline_data%x1)*invh + 1)
191 a = (spl_p(1)%spline_data%x1 - xx)*invh + real(i, kind=
dp)
194 ylo = spl_p(1)%spline_data%y(i)
195 yhi = spl_p(1)%spline_data%y(i + 1)
196 y2lo = spl_p(1)%spline_data%y2(i)
197 y2hi = spl_p(1)%spline_data%y2(i + 1)
198 potential_s = (a*ylo + b*yhi - ((a + 1.0_dp)*y2lo + (b + 1.0_dp)*y2hi)*a*b*h26)*spl_f%fscale(1)
199 y1 = invh*((yhi - ylo) + ((f13 - a*a)*y2lo - (f13 - b*b)*y2hi)*3.0_dp*h26)
200 y1 = 2.0_dp*y1*x4*spl_f%dscale(1)
219 TYPE(spline_data_type),
INTENT(IN) :: spl
220 REAL(kind=
dp),
INTENT(IN) :: xx
221 REAL(kind=
dp),
INTENT(OUT),
OPTIONAL :: y1
224 REAL(kind=
dp),
PARAMETER :: f13 = 1.0_dp/3.0_dp
227 REAL(kind=
dp) :: a, b, h26, invh, y2hi, y2lo, yhi, ylo
231 i = int((xx - spl%x1)*invh + 1)
233 a = (spl%x1 - xx)*invh + real(i, kind=
dp)
239 spline_value = a*ylo + b*yhi - ((a + 1.0_dp)*y2lo + (b + 1.0_dp)*y2hi)*a*b*h26
240 IF (
PRESENT(y1)) y1 = invh*((yhi - ylo) + &
241 ((f13 - a*a)*y2lo - (f13 - b*b)*y2hi)*3.0_dp*h26)
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Defines the basic variable types.
integer, parameter, public dp
routines for handling splines
pure subroutine, public init_splinexy(spl, nn)
allocates storage for function table to be interpolated both x and y are allocated
pure subroutine, public init_spline(spl, dx, y1a, y1b)
allocates storage for y2 table calculates y2 table and other spline parameters
real(kind=dp) function, public potential_s(spl_p, xxi, y1, spl_f, logger)
calculates the potential interpolated with splines value at a given point and the first derivative....
real(kind=dp) function, public spline_value(spl, xx, y1)
calculates the spline value at a given point (and possibly the first derivative) WITHOUT checks and w...
routines for handling splines_types