subroutine sphatmraddens
    use defvar
    use util
    use functions
    implicit real*8 (a-h,o-z)
    
    ! Declare variables
    real*8, allocatable :: radx(:), rady(:), radz(:), radpos(:), density(:)
    integer :: ntheta, nphi
    real*8 :: theta, phi
    integer :: irad, itheta, iphi, idx
    real*8 :: rad, x, y, z
    integer :: ifinish, iprogstp, iprogcrit, nradpt
    integer :: itmp
    parameter (truncrho = 1D-8, rlow = 0D0, rhigh = 12D0, nradpt = 200, ntheta = 50, nphi = 100)
    
    ! Allocate arrays
    allocate(radx(nradpt), rady(nradpt), radz(nradpt), radpos(nradpt), density(nradpt*ntheta*nphi))
    
    ifinish = 0
    iprogstp = 20
    iprogcrit = iprogstp
    write(*,*) "Calculating..."

    ! Parallel loop
    !$OMP PARALLEL DO SHARED(density, radpos, ifinish, iprogcrit) PRIVATE(irad, itheta, iphi, rad, x, y, z, theta, phi, idx) schedule(dynamic) NUM_THREADS(nthreads)
    do irad = 1, nradpt
        rad = rlow + (rhigh - rlow) * (irad - 1) / (nradpt - 1)
        radpos(irad) = rad

        do itheta = 0, ntheta-1
            theta = pi * itheta / (ntheta - 1) ! theta ranges from 0 to pi

            do iphi = 0, nphi-1
                phi = 2 * pi * iphi / nphi ! phi ranges from 0 to 2*pi

                x = rad * sin(theta) * cos(phi)
                y = rad * sin(theta) * sin(phi)
                z = rad * cos(theta)
                idx = (irad - 1) * ntheta * nphi + itheta * nphi + iphi + 1
                
                ! Compute density
                density(idx) = fdens(x, y, z)
                
            end do
        end do

        ifinish = ifinish + 1
        if (ifinish >= iprogcrit) then
            call showprog(ifinish, nradpt)
            iprogcrit = iprogcrit + iprogstp
        end if
    end do
    !$OMP END PARALLEL DO

    ! Output results
    open(10, file="density.txt", status="replace")
    itmp = 0
    do irad = 1, nradpt
        do itheta = 0, ntheta-1
            theta = pi * itheta / (ntheta - 1)
            do iphi = 0, nphi-1
                phi = 2 * pi * iphi / nphi
                idx = (irad - 1) * ntheta * nphi + itheta * nphi + iphi + 1
                if (density(idx) > truncrho) then
                    itmp = itmp + 1
                    write(10, "('    r=', f12.6, ' theta=', f12.6, ' phi=', f12.6, ' density=', f25.10, 'D0')") radpos(irad), theta, phi, density(idx)
                end if
            end do
        end do
    end do
    close(10)

    write(*,*) "The result has been output to density.txt in the current folder"
end subroutine sphatmraddens