spline_filter.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module spline_filter_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use splines_oct_m
  use messages_oct_m
  use namespace_oct_m
 
  implicit none
 
  private
 
  public ::               &
    spline_filter_ft,     &
    spline_filter_bessel, &
    spline_filter_mask,   &
    spline_filter_mask_init
 
  integer, parameter :: mask_n = 201
  real(real64) :: mask_x(mask_n), mask_y(mask_n)
 
contains
 
  !
  !----------------------------------------------------------------------------
  subroutine spline_filter_ft(spl, threshold, fs, rs)
    type(spline_t),      intent(inout) :: spl
    real(real64),        intent(in)    :: threshold
    real(real64), optional,     intent(in)    :: fs(2)
    real(real64), optional,     intent(in)    :: rs(2)
 
    type(spline_t) :: splw
 
    push_sub(spline_filter_ft)
 
    if (present(fs)) then
      call spline_init(splw)
      call spline_3dft(spl, splw, threshold, m_two*fs(1))
      call spline_cut(splw, fs(1), fs(2))
      call spline_3dft(splw, spl, threshold)
      call spline_times(real(m_one/(m_two*m_pi)**3, real64), spl, threshold)
      call spline_end(splw)
    end if
    if (present(rs)) then
      call spline_cut(spl, rs(1), rs(2), threshold)
    end if
 
    pop_sub(spline_filter_ft)
  end subroutine spline_filter_ft
 
 
  !----------------------------------------------------------------------------
  subroutine spline_filter_bessel(spl, l, qmax, alpha, beta_fs, rcut, beta_rs, threshold)
    type(spline_t), intent(inout) :: spl
    integer,        intent(in)    :: l
    real(real64),   intent(in)    :: qmax, alpha, beta_fs, rcut, beta_rs
    real(real64),   intent(in)    :: threshold
 
    type(spline_t) :: splw
 
    push_sub(spline_filter_bessel)
 
    call spline_init(splw)
    call spline_besselft(spl, splw, l, gmax=m_four*qmax)
    call spline_cut(splw, alpha*qmax, beta_fs)
    call spline_besselft(splw, spl, l)
    call spline_end(splw)
    call spline_cut(spl, rcut, beta_rs, threshold)
 
    pop_sub(spline_filter_bessel)
  end subroutine spline_filter_bessel
 
 
  !----------------------------------------------------------------------------
  subroutine spline_filter_mask_init()
 
    integer :: iunit, i
 
    push_sub(spline_filter_mask_init)
 
    iunit = io_open(trim(conf%share)//"/filter_mask.data", action='read', status='old', die=.true.)
 
    do i = 1, mask_n
      read(iunit, *) mask_x(i), mask_y(i)
    end do
 
    call io_close(iunit)
 
    pop_sub(spline_filter_mask_init)
  end subroutine spline_filter_mask_init
 
 
  !----------------------------------------------------------------------------
  subroutine spline_filter_mask(spl, l, rmax, qmax, alpha, gamma, threshold)
    type(spline_t), intent(inout) :: spl
    integer,        intent(in)    :: l
    real(real64),   intent(in)    :: rmax
    real(real64),   intent(in)    :: qmax
    real(real64),   intent(in)    :: alpha
    real(real64),   intent(in)    :: gamma
    real(real64),   intent(in)    :: threshold
 
    type(spline_t) :: mask, splw
    real(real64) :: local_mask_x(mask_n), rcut, beta
 
    push_sub(spline_filter_mask)
 
    rcut = gamma*rmax
 
    ! If rut is zero, the code below fails
    assert(rcut > m_epsilon)
 
    ! we define the mask function as f(r/rcut)
    local_mask_x = mask_x*rcut
    call spline_fit(mask_n, local_mask_x, mask_y, mask)
 
    ! apply the mask function
    call spline_div(spl, mask)
 
    ! transform to bessel space
    call spline_init(splw)
    call spline_besselft(spl, splw, l, gmax=m_four*qmax)
 
    ! apply a cutoff
    beta = log(1.e5_real64)/(alpha - m_one)**2
    call spline_cut(splw, qmax/alpha, beta)
 
    ! transform back to real space
    call spline_besselft(splw, spl, l)
    call spline_end(splw)
 
    ! remove the mask function
    call spline_mult(spl, mask, threshold)
 
    call spline_end(mask)
 
    pop_sub(spline_filter_mask)
 
  end subroutine spline_filter_mask
 
end module spline_filter_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: