marker_set_AddSphere Subroutine

    pure subroutine marker_set_AddSphere(this,center,radius,vel,dl,tag)
      !> Adds an IB sphere.
      class(marker_set),   intent(inout) :: this                               !! A collection of tessellation elements
      real(wp),            intent(in)    :: center(3)                          !! Sphere center
      real(wp),            intent(in)    :: radius                             !! Sphere radius
      real(wp),            intent(in)    :: vel(3)                             !! Sphere velocity
      real(wp),            intent(in)    :: dl                                 !! Element size
      integer(kind=8),     intent(in),    &
                                optional :: tag                                !! Tag
      ! Work variables
      real(wp)              :: dSA                                             !! Infinitesimal Surface area
      real(wp)              :: d_phi                                           !! Infinitesimal phi angle (0 -> 2pi)
      real(wp)              :: d_theta                                         !! Infinitesimal theta angle (0 -> pi)
      real(wp), allocatable :: theta(:)                                        !! theta angle (0 -> pi)
      integer               :: n_theta
      integer               :: n_phi
      real(wp)              :: phi_mid
      real(wp)              :: theta_mid
      integer               :: i,j
      integer               :: element_count
      integer(kind=8)       :: facetID
      real(wp), parameter   :: Pi=4.0_wp*atan(1.0_wp)
      integer               :: tag_

      if (present(tag)) then
        tag_ = int(tag,kind=4)
      else
        tag_ = 0
      end if

      ! Setup target surface area and phi angle
      dSA=dl**2

      ! Split the lattitude into segments of size=dl
      n_theta = ceiling(Pi*radius/dl)
      d_theta = Pi/real(n_theta,wp)

      ! Create the nodal values for the lattitude angle
      allocate(theta(n_theta+1))
      do i=1,n_theta+1
        theta(i)= real(i-1,wp)*d_theta
      end do

      ! Count the total number of elements to seed
      element_count = 0
      do i=1,n_theta
        theta_mid = 0.5_wp*(theta(i)+theta(i+1))

        ! Split the longitude such that the area of each
        ! element is approximately equal to dl^2
        d_phi = dSA/(radius**2*(cos(theta(i))-cos(theta(i+1)) ))

        ! Add to total number of elements
        element_count = element_count + ceiling(2.0_wp*Pi/d_phi)
      end do
      ! Create n_phi*n_theta markers
      call this%Resize(this%count_ + element_count)

      ! Place markers at center of facets
      facetID = this%count_ - element_count
      select type(markers =>this%p)
      type is (marker_obj)
        do i=1,n_theta

          ! Lattitude angle at mid point
          theta_mid = 0.5_wp*(theta(i)+theta(i+1))

          ! Discretize the longitudinal direction
          d_phi = dSA/(radius**2*(cos(theta(i))-cos(theta(i+1)) ))
          n_phi = ceiling(2.0_wp*Pi/d_phi)
          d_phi = 2.0_wp*Pi/real(n_phi,wp)

          do j=1,n_phi

            ! Longitude angle at mid point
            phi_mid = 0.5_wp*(real(j-1,wp)+real(j,wp))*d_phi

            ! Facet ID
            facetID = facetID + int(1,kind=8)
            markers(facetID)%id = facetID

            ! Position
            markers(facetID)%p(1) = center(1) + radius*sin(theta_mid)*cos(phi_mid)
            markers(facetID)%p(2) = center(2) + radius*sin(theta_mid)*sin(phi_mid)
            markers(facetID)%p(3) = center(3) + radius*cos(theta_mid)

            ! Surface area
            markers(facetID)%SA = radius**2*d_phi *(cos(theta(i))-cos(theta(i+1)))

            ! Normal vector
            markers(facetID)%n = [sin(theta_mid)*cos(phi_mid),sin(theta_mid)*sin(phi_mid),cos(theta_mid)]

            ! Velocity
            markers(facetID)%v = vel

            ! Tag
            markers(facetID)%s = tag_

            ! Force
            markers(facetID)%f = 0.0_wp

          end do
        end do
      end select

      deallocate(theta)

      return
    end subroutine marker_set_AddSphere