shape_regions Subroutine

    subroutine shape_regions(self, r, rs, d)
        type(ShadowCones), intent(inout) :: self
        real(wp), intent(in) :: r, rs, d                    !! scene units
        real(wp) :: sa, ca, ta, sb, cb, tb, apex, len
        integer :: k
        self%region%drawn = .false.
        self%drawn_length_ = 0.0_wp
        if (r > 0.0_wp .and. rs > 0.0_wp .and. d > r + rs) then
            sa = (rs - r) / d                               ! the umbra's half-angle (sine; < 0: widening)
            ca = sqrt(1.0_wp - sa**2)
            ta = sa / ca
            sb = (rs + r) / d                               ! the penumbra's
            cb = sqrt(1.0_wp - sb**2)
            tb = sb / cb
            apex = huge(1.0_wp)
            if (sa > 0.0_wp) apex = r / sa
            len = self%length_
            if (len <= 0.0_wp) then
                len = 100.0_wp * r
                if (sa > 0.0_wp) len = 2.0_wp * apex
            end if
            self%drawn_length_ = len
            ! the umbra: from the outer tangents' circle on the body to the apex (or len)
            associate (g => self%region(umbra))
                g%z0 = r * sa
                g%r0 = r * ca
                g%z1 = min(apex, len)
                g%r1 = max(0.0_wp, g%r0 - (g%z1 - g%z0) * ta)
                if (g%z1 >= apex) g%r1 = 0.0_wp
                g%drawn = g%z1 > g%z0
            end associate
            ! the antumbra: from the apex to len
            associate (g => self%region(antumbra))
                g%z0 = apex
                g%r0 = 0.0_wp
                g%z1 = len
                g%r1 = (len - apex) * ta
                g%drawn = sa > 0.0_wp .and. len > apex
            end associate
            ! the penumbra: from the inner tangents' circle (on the light's side) to len
            associate (g => self%region(penumbra))
                g%z0 = -r * sb
                g%r0 = r * cb
                g%z1 = len
                g%r1 = g%r0 + (len - g%z0) * tb
                g%drawn = len > g%z0
            end associate
        end if
        self%region(umbra)%drawn = self%region(umbra)%drawn .and. iand(self%regions_, SHADOW_UMBRA) /= 0
        self%region(penumbra)%drawn = self%region(penumbra)%drawn .and. iand(self%regions_, SHADOW_PENUMBRA) /= 0
        self%region(antumbra)%drawn = self%region(antumbra)%drawn .and. iand(self%regions_, SHADOW_ANTUMBRA) /= 0
        do k = 1, n_regions
            call write_region(self%region(k))
        end do
    end subroutine shape_regions