submodule_self_force_observer_implementation.f90 Source File

submodule(self_force_observer) self_force_observer_implementation
!! The implementation of the interface in [[self_force_observer]].

  implicit none

contains

  module procedure sf_init

    use parameters, only : lmin, lmax, fit_high_l, first_l, last_l
    implicit none

    integer(ip) :: n_unique
    select type(object)
    type is (scal_schw)
      this%p => object
    class default
      print*,'Self-force observer initialized with invalid type'
      stop
    end select

    allocate ( this%vname, source = this%p%ename )

    if (size(rad)/=1) then
      print*,'Self-force observer should be initialized with 1 extraction radius'
      stop
    end if
    this%nradii = size(rad)
    this%radii = rad
    this%ioo_id = -1

    call find_indices ( rad, coord, this%elem_index, this%node_index )

    n_unique = n_unique_values ( this%p%ll )
    if (.not. fit_high_l) then
      allocate ( this%fl(lmin:lmax), this%ftl(lmin:lmax), &
                 this%fphil(lmin:lmax), this%frl(lmin:lmax) )
    else
      allocate ( this%fl(lmin:lmax+3), this%ftl(lmin:lmax+3), &
                 this%fphil(lmin:lmax+3), this%frl(lmin:lmax+3) )
    end if
    if ( fit_high_l ) then
      if ( lmin>0 ) then
        print*,'fit_high_l should only be .true. when lmin=0'
        stop
      end if
      if ( first_l < lmin ) then
        print*,'Error: first_l < lmin'
        stop
      end if
      if ( last_l > lmax ) then
        print*,'Error: last_l < lmax'
        stop
      end if
    end if

  end procedure sf_init


  module procedure sf_extract

    use orbit_base
    use gsl_interface, only : legendre_sphPlm
    use parameters, only : mass, use_generic_orbit, lmax, fit_high_l, &
                           first_l, last_l, fit_order, nfit
    use numerics, only : correct_for_higher_modes
    use time_info, only : get_current_time

    implicit none

    real(wp) :: r, phi, ur, En, Lz, m_fold_factor, y_lm, ftl_new, frl_new
    complex(wp) :: phifactor, phase, psip, psitp, psirp
    integer(ip) :: i, j, k

    i = this%elem_index(1)
    j = this%node_index(1)
    this%fl = 0.0_wp
    this%ftl = 0.0_wp
    this%fphil = 0.0_wp
    this%frl = 0.0_wp
!    print*,'SF extract at time = ', get_current_time()

    call orbit_info%get_orbit ( r, phi, ur, En, Lz )    
!    print*,'phi = ', phi

    associate ( nmodes => this%p%nmodes, &
                mm => this%p%mm, &
                ll => this%p%ll, &
                var => this%p%eq_data, &
                fl => this%fl, &
                ftl => this%ftl, &
                fphil => this%fphil, &
                frl => this%frl )
      do k = 1, nmodes
        if ( mm(k) == 0 ) then
          m_fold_factor = 1.0_wp
        else
          m_fold_factor = 2.0_wp
        end if

        phifactor = zi*mm(k)
        phase = cmplx ( cos(mm(k)*phi), sin(mm(k)*phi), wp )
        y_lm = legendre_sphPlm ( ll(k), mm(k), 0.0_wp )
        psip = 0.5_wp * ( var(1,k)%elems(i)%var(j) + &
                          var(1,k)%elems(i+1)%var(1) )
        psitp = 0.5_wp * ( var(2,k)%elems(i)%var(j) + &
                           var(2,k)%elems(i+1)%var(1) )
        psirp = 0.5_wp * ( var(3,k)%elems(i)%var(j) + &
                           var(3,k)%elems(i+1)%var(1) )
!        if (k<=16) then
!          print*,'k = ', k,': ', psitp, ftl(ll(k))
!        end if
        fl(ll(k)) = fl(ll(k)) + m_fold_factor*y_lm*real(phase*psip,wp)
        ftl(ll(k)) = ftl(ll(k)) + m_fold_factor*y_lm*real(phase*psitp,wp)
!        if (k<=16) then
!          print*,m_fold_factor*y_lm,m_fold_factor*y_lm*real(phase*psitp,wp),ftl(ll(k))
!        end if
        fphil(ll(k)) = fphil(ll(k)) &
                      + m_fold_factor*y_lm*real(phifactor*phase*psip,wp)
        frl(ll(k)) = frl(ll(k)) + m_fold_factor*y_lm*real(phase*psirp,wp)

      end do 
!      print*,'ftl = ', ftl(0:20)
      do k = 0, lmax
        if (use_generic_orbit) then
          call this%p%time_dep_coord%tdc_to_tortoise ( i, 2, ftl(k), &
                                                       frl(k), &
                                                       ftl_new, frl_new )
          ftl(k) = ftl_new
          frl(k) = frl_new
        end if
      end do
!      print*,'ftl = ', ftl(0:20)

      frl = frl/(r-2.0_wp*mass)-fl/r**2
      fl = fl/r
      ftl = ftl/r
      fphil = fphil/r

      if (fit_high_l) then
        fl(lmax+1:lmax+3) = correct_for_higher_modes ( fl(0:lmax), &
                                    lmax+1, first_l, last_l, fit_order, nfit )
        ftl(lmax+1:lmax+3) = correct_for_higher_modes ( ftl(0:lmax), &
                                    lmax+1, first_l, last_l, fit_order, nfit )
        frl(lmax+1:lmax+3) = correct_for_higher_modes ( frl(0:lmax), &
                                    lmax+1, first_l, last_l, fit_order, nfit )
        fphil(lmax+1:lmax+3) = correct_for_higher_modes ( fphil(0:lmax), &
                                    lmax+1, first_l, last_l, fit_order, nfit )
      end if
    end associate

  end procedure sf_extract


  module procedure sf_output

    use output_base
    use time_info

    implicit none

    integer(ip) :: ioo_id, tmp_id
    character(len=len_trim(this%vname)+10) :: filename

    ioo_id = this%ioo_id
    if (ioo_id<0) then
      filename = trim(this%vname) // '.fl.asc'
      ioo_id = next_available_io_id ()
      this%ioo_id = ioo_id
      print*,'Opening ', filename, ' with id ', ioo_id
      open(ioo_id, file=filename, status='replace', action='write')
      filename = trim(this%vname) // '.ftl.asc'
      tmp_id = next_available_io_id ()
      print*,'Opening ', filename, ' with id ', tmp_id
      open(tmp_id, file=filename, status='replace', action='write')
      filename = trim(this%vname) // '.fphil.asc'
      tmp_id = next_available_io_id ()
      print*,'Opening ', filename, ' with id ', tmp_id
      open(tmp_id, file=filename, status='replace', action='write')
      filename = trim(this%vname) // '.frl.asc'
      tmp_id = next_available_io_id ()
      print*,'Opening ', filename, ' with id ', tmp_id
      open(tmp_id, file=filename, status='replace', action='write')
    end if
    write(ioo_id,'(*(es23.15e3,1x))') get_current_time ( ), this%fl(:) 
    write(ioo_id+1,'(*(es23.15e3,1x))') get_current_time ( ), this%ftl(:) 
    write(ioo_id+2,'(*(es23.15e3,1x))') get_current_time ( ), this%fphil(:) 
    write(ioo_id+3,'(*(es23.15e3,1x))') get_current_time ( ), this%frl(:) 
    
  end procedure sf_output


  module procedure close_sf_observer
  end procedure close_sf_observer


  function n_unique_values ( var ) result(n)
  !! Helper function that finds the number of unique values in a 1d integer
  !! array.
  !!
  !! This function does not seem to be used, so can probably safely be
  !! removed.

    use iso_c_binding

    integer(c_int), dimension(:), intent(in) :: var
    !! An 1d integer arrays.
    integer(ip) :: n
    !! The return value is the number of unique values in the input array.
    integer(c_int), dimension(size(var)) :: tmp
    integer(ip) :: i

    tmp(1) = var(1)
    n = 1

    do i=2, size(var)
      if ( .not. any ( tmp(1:n) == var(i) ) ) then
        n = n+1
        tmp(n) = var(i)
      end if
    end do
   
  end function n_unique_values

end submodule self_force_observer_implementation