Galactic rotation curve Newton's calc

gistfile1.txt

using HCubature, LinearAlgebra, StatsPlots, QuadGK, StaticArrays



# Define a function that describes the density 
# as a function of the radial distance from the center of the galaxy. This is a
# dimensionless quantity (mass/area) / (mass/area at zero distance)
# r is a dimensionless distance, radius as a fraction of the galactic scale radius 
function rho(x) 
    exp(-sqrt(x[1]^2+x[2]^2)-5.0*abs(x[3]))
end


## given a test particle at position r which is a vector [a,b] 
## (I use [r,0.0] when I call the function below)
## return a function f(x) which when given another point x as a vector 
## in the galaxy computes the vector
## force that the point x exerts on the test particle at r. 
## Note that we ignore all mass within radius around .01 of the test particle, otherwise we have
## a singularity of "self gravitation" at zero distance.

function integrand(r,eps)
   function f(x) 
        n = norm(x-r)
        #(1.0-exp(-(n/eps)^2))* # no self gravity 
        (n > eps)*
            rho(x)*(x-r)/n^3 # cubed because (r-x)/norm(r-x) is the direction 1/norm(r-x)^2 is magnitude
   end
   f
end

## hcubature is a function which when given a function and an N dimensional "box" 
## (by the coordinates of two opposite corners of the box) computes the integral over the entire
## box of the integrand times differential volume in N dimensions. This is a test call

hcubature(integrand(SVector(1.0,0.0,0.0),.005),[-100.0,-100.0,-1.0],[100.0,100.0,1.0]; maxevals=100_000,rtol=.00125)



function rho2(x) 
    r = norm(x)
    return 1.0/(1.0+r) #*(r < 10.0)
end

function integrand2(r,eps)
    function f(x) 
         n = norm(x-r)
         (n > eps)* # no self gravity 
             rho2(x)*(x-r)/n^3 # cubed because (r-x)/norm(r-x) is the direction 1/norm(r-x)^2 is magnitude
    end
    f
 end
 

 hcubature(integrand2(SVector(1.2,0.0,0.0),.1),[-2.0,-2.0,-.10],[2.0,2.0,.10]; maxevals=100_000,rtol=.00125)




# try out some radii, from 0.1 to 10 galactic scale lengths

rvals = collect(0.1:0.1:13.0)

# calculate the accelerations experienced by test particle at location r for each r in rvals

accels = map(r -> hcubature(integrand([r,0.0,0.0],.02),[-100.0,-100.0,-5.0],[100.0,100.0,5.0]; maxevals=1_000_000,rtol=0.000125),rvals)

## calculate the mass of the disk inside the radius R for each R in rvals using quadgk
## Since this is a calculation over a disk, it's easier to perform along a radius rather than hcubature
## over a circular region

accels2 = map(r -> hcubature(integrand2(SVector(r,0.0,0.0),.002),[-20.0,-20.0,-.20],[20.0,20.0,.20],maxevals=1_000_000,rtol=.000125),rvals)


masses = map(R->quadgk(r -> rho([0.0,r,0.0])*2*pi*r,0,R)[1],rvals)


## plots of the various quantities

plot(rvals,[sqrt(rvals[i]*norm(accels[i][1])) for i in eachindex(rvals)]; 
    title= "Velocity vs Radius", xlab="Radius (fraction of avg galactic radius)", ylab="Velocity (dimensionless)") |> display

plot(rvals,[sqrt(rvals[i]*norm(accels2[i][1])) for i in eachindex(rvals)]; 
    title= "Velocity vs Radius", xlab="Radius (fraction of avg galactic radius)", ylab="Velocity (dimensionless)") |> display

plot(rvals,masses; title="Masses at Radius R",xlab="R",ylab="Mass") |>display