correct_for_higher_modes Function

  1. module_numerics.f90
  2. numerics
  3. correct_for_higher_modes

public function correct_for_higher_modes(ydat, nmodes, startfit, endfit, order, npar)

Uses

  • proc~~correct_for_higher_modes~~UsesGraph proc~correct_for_higher_modes correct_for_higher_modes module~gsl_interface gsl_interface proc~correct_for_higher_modes->module~gsl_interface module~kinds kinds module~gsl_interface->module~kinds iso_c_binding iso_c_binding module~gsl_interface->iso_c_binding

Function that fits higher order terms to data containing amplitude as function of l-values and corrects the sum over the evolved l-modes with the contribution from the not evolved l-modes.

Arguments

Type IntentOptional AttributesName
real(kind=wp), intent(in), dimension(:):: ydat

A 1d-array containing the amplitudes of the -modes. It is assumed, for now, that the data starts at .
integer(kind=ip), intent(in) :: nmodes

The number of -modes.
integer(kind=ip), intent(in) :: startfit

At what -value should the fit begin.
integer(kind=ip), intent(in) :: endfit

At what -value should the fit end.
integer(kind=ip), intent(in) :: order

What is the lowest order term that should be fitted.
integer(kind=ip), intent(in) :: npar

How many terms should be included in the correction. One more term is included in the fit, but the coefficient for the last fit coefficient is only used as an error estimator.

Return Value real(kind=wp), dimension(3)

A 1d-array of size 3 that on return contains the sum of evolved modes (first entry), the corrected sum (second entry) and the error estimate (third entry).


Calls

proc~~correct_for_higher_modes~~CallsGraph proc~correct_for_higher_modes correct_for_higher_modes proc~ldep ldep proc~correct_for_higher_modes->proc~ldep proc~lsum lsum proc~correct_for_higher_modes->proc~lsum

Called by

proc~~correct_for_higher_modes~~CalledByGraph proc~correct_for_higher_modes correct_for_higher_modes proc~sf_extract sf_extract proc~sf_extract->proc~correct_for_higher_modes interface~sf_extract sf_extract interface~sf_extract->proc~sf_extract

Contents

Source Code

correct_for_higher_modes

Source Code

  function correct_for_higher_modes ( ydat, nmodes, &
                                      startfit, endfit, order, npar )
  !! Function that fits higher order terms to data containing amplitude as
  !! function of l-values and corrects the sum over the evolved l-modes with
  !! the contribution from the not evolved l-modes.

    use gsl_interface

    real(wp), dimension(:), intent(in) :: ydat
    !! A 1d-array containing the amplitudes of the \(\ell\)-modes. It is
    !! assumed, for now, that the data starts at \(\ell=0\).
    integer(ip), intent(in) :: nmodes
    !! The number of \(\ell\)-modes.
    integer(ip), intent(in) :: startfit
    !! At what \(\ell\)-value should the fit begin.
    integer(ip), intent(in) :: endfit
    !! At what \(\ell\)-value should the fit end.
    integer(ip), intent(in) :: order
    !! What is the lowest order term that should be fitted.
    integer(ip), intent(in) :: npar
    !! How many terms should be included in the correction. One more term
    !! is included in the fit, but the coefficient for the last fit
    !! coefficient is only used as an error estimator.
    real(wp), dimension(3) :: correct_for_higher_modes
    !! A 1d-array of size 3 that on return contains the sum of evolved modes
    !! (first entry), the corrected sum (second entry) and the error
    !! estimate (third entry).
    integer(ip) :: n
    real(wp), dimension(npar+1,endfit-startfit+1) :: xm
    real(wp), dimension(endfit-startfit+1) :: yv
    real(wp), dimension(npar+1) :: cv
    real(wp), dimension(npar+1,npar+1) :: covm
    real(wp) :: chisq
    integer(ip) :: ierr, i, j, l
    real(wp) :: c1, c2, c3

    n = endfit-startfit+1
    do i = 1, n
      l = startfit+i-1
      yv(i) = ydat(startfit+i)
      do j = 1, npar+1
        xm(j,i) = ldep(l,order+2*(j-1))
      end do
    end do

    call multifit_linear ( n, npar+1, reshape(xm,(/n*(npar+1)/)), yv, cv, &
                           reshape(covm,(/(npar+1)*(npar+1)/)), chisq, ierr )

    correct_for_higher_modes(1) = sum(ydat)
    correct_for_higher_modes(2) = correct_for_higher_modes(1)
    do j = 1, npar
      correct_for_higher_modes(2) = correct_for_higher_modes(2) + &
                                   cv(j)*lsum(nmodes,order+2*(j-1))
    end do
    correct_for_higher_modes(3) = cv(npar+1)*lsum(nmodes,order+2*npar)
  end function correct_for_higher_modes