feldc Subroutine

private subroutine feldc(me, v, b)

Alternate version of feldg to be used with cartesian coordinates

Type Bound

shellig_type

Arguments

Type IntentOptional Attributes Name
class(shellig_type), intent(inout) :: me
real(kind=wp), intent(in), dimension(3) :: v

cartesian coordinates in earth radii (6371.2 km) x-axis pointing to equator at 0 longitude y-axis pointing to equator at 90 long. z-axis pointing to north pole
real(kind=wp), intent(out) :: b(3)

field components

Called by

proc~~feldc~~CalledByGraph proc~feldc shellig_module::shellig_type%feldc proc~igrfc shellig_module::shellig_type%igrfc proc~igrfc->proc~feldc proc~get_flux_c_ radbelt_module::radbelt_type%get_flux_c_ proc~get_flux_c_->proc~igrfc none~get_flux radbelt_module::radbelt_type%get_flux none~get_flux->proc~get_flux_c_ proc~get_flux_c radbelt_module::get_flux_c proc~get_flux_c->none~get_flux proc~get_flux_g radbelt_module::get_flux_g proc~get_flux_g->none~get_flux proc~get_flux_g_c radbelt_c_module::get_flux_g_c proc~get_flux_g_c->none~get_flux interface~get_flux radbelt_module::get_flux interface~get_flux->proc~get_flux_c interface~get_flux->proc~get_flux_g

Source Code

    subroutine feldc(me, v, b)

        class(shellig_type), intent(inout) :: me
        real(wp), dimension(3), intent(in) :: v !! cartesian coordinates in earth radii (6371.2 km)
                                          !! x-axis pointing to equator at 0 longitude
                                          !! y-axis pointing to equator at 90 long.
                                          !! z-axis pointing to north pole
        real(wp), intent(out) :: b(3) !! field components

        real(wp) :: f, rq, s, t, x, xxx, y, yyy, z, zzz
        integer :: i, ih, ihmax, il, imax, k, last, m

        xxx = v(1)
        yyy = v(2)
        zzz = v(3)

        rq = 1.0_wp / (xxx * xxx + yyy * yyy + zzz * zzz)
        me%xi = [xxx, yyy, zzz] * rq

        ihmax = me%nmax * me%nmax + 1
        last = ihmax + me%nmax + me%nmax
        imax = me%nmax + me%nmax - 1
        do i = ihmax, last
            me%h(i) = me%g(i)
        end do
        do k = 1, 3, 2
            i = imax
            ih = ihmax
            do
                il = ih - i
                f = 2.0_wp / real(i - k + 2, wp)
                x = me%xi(1) * f
                y = me%xi(2) * f
                z = me%xi(3) * (f + f)
                i = i - 2
                if ((i - 1) >= 0) then
                    if ((i - 1) > 0) then
                        do m = 3, i, 2
                            me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - &
                                                                       me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2))
                            me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - &
                                                                       me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1))
                        end do
                    end if
                    me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih))
                    me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih))
                end if
                me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2))
                ih = il
                if (i < k) exit
            end do
        end do

        s = 0.5_wp * me%h(1) + 2.0_wp * (me%h(2) * me%xi(3) + me%h(3) * me%xi(1) + me%h(4) * me%xi(2))
        t = (rq + rq) * sqrt(rq)

        b(1) = t * (me%h(3) - s * xxx)
        b(2) = t * (me%h(4) - s * yyy)
        b(3) = t * (me%h(2) - s * zzz)

    end subroutine feldc