trgev Subroutine

private subroutine trgev(cphi, sphi, p, vcos, vsin)

Compute arrays of \f$ \cos(m \phi) \f$ and \f$ \sin(m \phi) \f$

All values are computed recurrently from input \f$ \cos(\phi) \f$ and \f$ \sin(\phi) \f$ without accessing arccos or arcsin functions.

@param[in] cphi: \f$ \cos(\phi) \f$. -1 <= cphi <= 1 @param[in] sphi: \f$ \sin(\phi) \f$. -1 <= sphi <= 1 @param[in] p: Maximal value of \f$ m \f$, for which to compute \f$ \cos(m \phi) \f$ and \f$ \sin(m\phi) \f$. p >= 0 @param[out] vcos: Array of \f$ \cos(m\phi) \f$ for \f$ m=0..p \f$. Dimension is (p+1) @param[out] vsin: Array of \f$ \sin(m\phi) \f$ for \f$ m=0..p \f$. Dimension is (p+1)

Arguments

Type IntentOptional Attributes Name
real(kind=rp), intent(in) :: cphi
real(kind=rp), intent(in) :: sphi
integer, intent(in) :: p
real(kind=rp), intent(out) :: vcos(p+1)
real(kind=rp), intent(out) :: vsin(p+1)

Called by

proc~~trgev~~CalledByGraph proc~trgev trgev proc~fmm_m2m_rotation_work fmm_m2m_rotation_work proc~fmm_m2m_rotation_work->proc~trgev proc~fmm_l2l_rotation_work fmm_l2l_rotation_work proc~fmm_l2l_rotation_work->proc~trgev proc~fmm_m2l_rotation_work fmm_m2l_rotation_work proc~fmm_m2l_rotation_work->proc~trgev proc~fmm_m2m fmm_m2m proc~fmm_m2m->proc~fmm_m2m_rotation_work proc~fmm_l2l fmm_l2l proc~fmm_l2l->proc~fmm_l2l_rotation_work proc~fmm_m2l fmm_m2l proc~fmm_m2l->proc~fmm_m2l_rotation_work proc~tree_p2m tree_p2m proc~tree_p2m->proc~fmm_m2m proc~cart_propfar_at_ipart cart_propfar_at_ipart proc~cart_propfar_at_ipart->proc~fmm_l2l proc~tree_m2m tree_m2m proc~tree_m2m->proc~fmm_m2m proc~tree_m2l tree_m2l proc~tree_m2l->proc~fmm_m2l proc~tree_l2l tree_l2l proc~tree_l2l->proc~fmm_l2l proc~cart_propnear_at_ipart cart_propnear_at_ipart proc~cart_propnear_at_ipart->proc~fmm_m2l proc~fmm_solve_for_multipoles fmm_solve_for_multipoles proc~fmm_solve_for_multipoles->proc~tree_p2m proc~fmm_solve_for_multipoles->proc~tree_m2m proc~fmm_solve_for_multipoles->proc~tree_m2l proc~fmm_solve_for_multipoles->proc~tree_l2l proc~elec_prop_d2m elec_prop_D2M proc~elec_prop_d2m->proc~cart_propfar_at_ipart proc~prepare_fmm_ipd prepare_fmm_ipd proc~elec_prop_d2m->proc~prepare_fmm_ipd proc~field_extd2d field_extD2D proc~field_extd2d->proc~cart_propfar_at_ipart proc~prepare_fmm_ext_ipd prepare_fmm_ext_ipd proc~field_extd2d->proc~prepare_fmm_ext_ipd proc~elec_prop_d2d elec_prop_D2D proc~elec_prop_d2d->proc~cart_propfar_at_ipart proc~elec_prop_d2d->proc~prepare_fmm_ipd proc~elec_prop_m2m elec_prop_M2M proc~elec_prop_m2m->proc~cart_propfar_at_ipart proc~preapare_fmm_static preapare_fmm_static proc~elec_prop_m2m->proc~preapare_fmm_static proc~elec_prop_m2d elec_prop_M2D proc~elec_prop_m2d->proc~cart_propfar_at_ipart proc~elec_prop_m2d->proc~preapare_fmm_static proc~cart_prop_at_ipart cart_prop_at_ipart proc~cart_prop_at_ipart->proc~cart_propfar_at_ipart proc~cart_prop_at_ipart->proc~cart_propnear_at_ipart proc~fmm_solve fmm_solve proc~fmm_solve->proc~tree_m2m proc~fmm_solve->proc~tree_m2l proc~fmm_solve->proc~tree_l2l proc~prepare_fmm_ext_ipd->proc~fmm_solve_for_multipoles proc~preapare_fmm_static->proc~fmm_solve_for_multipoles proc~prepare_polelec prepare_polelec proc~prepare_polelec->proc~elec_prop_d2m proc~prepare_polelec->proc~elec_prop_d2d proc~prepare_polelec->proc~elec_prop_m2d proc~prepare_fixedelec prepare_fixedelec proc~prepare_fixedelec->proc~elec_prop_m2m proc~tmatvec_otf TMatVec_otf proc~tmatvec_otf->proc~field_extd2d proc~prepare_fmm_ipd->proc~prepare_fmm_ext_ipd proc~ommp_get_polelec_energy ommp_get_polelec_energy proc~ommp_get_polelec_energy->proc~prepare_polelec proc~energy_mm_pol energy_MM_pol proc~ommp_get_polelec_energy->proc~energy_mm_pol proc~ommp_set_external_field ommp_set_external_field proc~ommp_set_external_field->proc~prepare_polelec proc~polelec_geomgrad polelec_geomgrad proc~polelec_geomgrad->proc~prepare_polelec proc~energy_mm_pol->proc~prepare_polelec proc~fixedelec_geomgrad fixedelec_geomgrad proc~fixedelec_geomgrad->proc~prepare_fixedelec proc~energy_mm_mm energy_MM_MM proc~energy_mm_mm->proc~prepare_fixedelec proc~ommp_get_full_ele_energy ommp_get_full_ele_energy proc~ommp_get_full_ele_energy->proc~ommp_get_polelec_energy proc~ommp_get_fixedelec_energy ommp_get_fixedelec_energy proc~ommp_get_full_ele_energy->proc~ommp_get_fixedelec_energy proc~c_ommp_get_polelec_energy C_ommp_get_polelec_energy proc~c_ommp_get_polelec_energy->proc~ommp_get_polelec_energy proc~ommp_set_external_field_nomm ommp_set_external_field_nomm proc~ommp_set_external_field_nomm->proc~ommp_set_external_field proc~c_ommp_set_external_field C_ommp_set_external_field proc~c_ommp_set_external_field->proc~ommp_set_external_field proc~c_ommp_set_external_field_nomm C_ommp_set_external_field_nomm proc~c_ommp_set_external_field_nomm->proc~ommp_set_external_field proc~ommp_polelec_geomgrad ommp_polelec_geomgrad proc~ommp_polelec_geomgrad->proc~polelec_geomgrad proc~ommp_full_geomgrad ommp_full_geomgrad proc~ommp_full_geomgrad->proc~polelec_geomgrad proc~ommp_full_geomgrad->proc~fixedelec_geomgrad proc~ommp_fixedelec_geomgrad ommp_fixedelec_geomgrad proc~ommp_fixedelec_geomgrad->proc~fixedelec_geomgrad proc~ommp_get_fixedelec_energy->proc~energy_mm_mm proc~ommp_get_full_energy ommp_get_full_energy proc~ommp_get_full_energy->proc~ommp_get_full_ele_energy proc~c_ommp_polelec_geomgrad C_ommp_polelec_geomgrad proc~c_ommp_polelec_geomgrad->proc~ommp_polelec_geomgrad proc~c_ommp_full_geomgrad C_ommp_full_geomgrad proc~c_ommp_full_geomgrad->proc~ommp_full_geomgrad proc~c_ommp_fixedelec_geomgrad C_ommp_fixedelec_geomgrad proc~c_ommp_fixedelec_geomgrad->proc~ommp_fixedelec_geomgrad proc~c_ommp_get_full_ele_energy C_ommp_get_full_ele_energy proc~c_ommp_get_full_ele_energy->proc~ommp_get_full_ele_energy proc~c_ommp_get_fixedelec_energy C_ommp_get_fixedelec_energy proc~c_ommp_get_fixedelec_energy->proc~ommp_get_fixedelec_energy proc~c_ommp_get_full_energy C_ommp_get_full_energy proc~c_ommp_get_full_energy->proc~ommp_get_full_energy

Contents

Source Code

trgev

Source Code

    subroutine trgev(cphi, sphi, p, vcos, vsin)
        ! Inputs
        integer, intent(in) :: p
        real(rp), intent(in) :: cphi, sphi
        ! Output
        real(rp), intent(out) :: vcos(p+1), vsin(p+1)
        ! Local variables
        integer :: m
        real(rp) :: c4phi, s4phi
        !! Treat values of p from 0 to 3 differently
        select case(p)
            ! Do nothing if p < 0
            case (:-1)
                return
            ! p = 0
            case (0)
                vcos(1) = 1.0
                vsin(1) = 0.0
                return
            ! p = 1
            case (1)
                vcos(1) = 1.0
                vsin(1) = 0.0
                vcos(2) = cphi
                vsin(2) = sphi
                return
            ! p = 2
            case (2)
                vcos(1) = 1.0
                vsin(1) = 0.0
                vcos(2) = cphi
                vsin(2) = sphi
                vcos(3) = cphi**2 - sphi**2
                vsin(3) = 2 * cphi * sphi
                return
            ! p = 3
            case (3)
                vcos(1) = 1.0
                vsin(1) = 0.0
                vcos(2) = cphi
                vsin(2) = sphi
                vcos(3) = cphi**2 - sphi**2
                vsin(3) = 2 * cphi * sphi
                vcos(4) = vcos(3)*cphi - vsin(3)*sphi
                vsin(4) = vcos(3)*sphi + vsin(3)*cphi
                return
            ! p >= 4
            case default
                vcos(1) = 1.0
                vsin(1) = 0.0
                vcos(2) = cphi
                vsin(2) = sphi
                vcos(3) = cphi**2 - sphi**2
                vsin(3) = 2 * cphi * sphi
                vcos(4) = vcos(3)*cphi - vsin(3)*sphi
                vsin(4) = vcos(3)*sphi + vsin(3)*cphi
                vcos(5) = vcos(3)**2 - vsin(3)**2
                vsin(5) = 2 * vcos(3) * vsin(3)
                c4phi = vcos(5)
                s4phi = vsin(5)
        end select
        ! Define cos(m*phi) and sin(m*phi) recurrently 4 values at a time
        do m = 6, p-2, 4
            vcos(m:m+3) = vcos(m-4:m-1)*c4phi - vsin(m-4:m-1)*s4phi
            vsin(m:m+3) = vcos(m-4:m-1)*s4phi + vsin(m-4:m-1)*c4phi
        end do
        ! Work with leftover
        select case(m-p)
            case (-1)
                vcos(p-1) = vcos(p-2)*cphi - vsin(p-2)*sphi
                vsin(p-1) = vcos(p-2)*sphi + vsin(p-2)*cphi
                vcos(p) = vcos(p-2)*vcos(3) - vsin(p-2)*vsin(3)
                vsin(p) = vcos(p-2)*vsin(3) + vsin(p-2)*vcos(3)
                vcos(p+1) = vcos(p-2)*vcos(4) - vsin(p-2)*vsin(4)
                vsin(p+1) = vcos(p-2)*vsin(4) + vsin(p-2)*vcos(4)
            case (0)
                vcos(p) = vcos(p-1)*cphi - vsin(p-1)*sphi
                vsin(p) = vcos(p-1)*sphi + vsin(p-1)*cphi
                vcos(p+1) = vcos(p-1)*vcos(3) - vsin(p-1)*vsin(3)
                vsin(p+1) = vcos(p-1)*vsin(3) + vsin(p-1)*vcos(3)
            case (1)
                vcos(p+1) = vcos(p)*cphi - vsin(p)*sphi
                vsin(p+1) = vcos(p)*sphi + vsin(p)*cphi
        end select
    end subroutine trgev