4914.sagews

# Calculating the disc radial integral
Rcore = var('Rcore')
RM = var('RM')
L = var('L')
r = var('r')
assume(RM>0)
f(Rcore,RM,L) = integral(r*exp(-(r-Rcore)/L),r,Rcore,RM)+integral(r,r,0,Rcore)
f

(Rcore, RM, L) |--> -(L^2*e^(Rcore/L) + L*RM*e^(Rcore/L))*e^(-RM/L) + L^2 + L*Rcore + 1/2*Rcore^2

N(f(16000,5200*10,5200))

2.37947082535011e8

# Calculating the disc z integral
ZM = var('ZM')
ZS = var('ZS')
z = var('z')
assume(ZM>0)
g(ZM,ZS) = integral(exp(-abs(z)/ZS),z,-ZM,ZM)
g

(ZM, ZS) |--> 2*(ZS*e^(ZM/ZS) - ZS)*e^(-ZM/ZS)

N(g(950*25,950))

1899.99999997361

# Total disc integral
Ndisc = N(2*pi()*f(16000,130000,5200)*g(20000,950))
Ndisc

2.84412152383230e12

# Calculating the bulge radial integral
fb = var('fb')
Lb = var('Lb')
Rbm = var('Rbm')
assume(Rbm>0)
h(Lb,Rbm) = integral(r^2*exp(-r/Lb),r,0,Rbm)
h

(Lb, Rbm) |--> 2*Lb^3 - (2*Lb^3 + 2*Lb^2*Rbm + Lb*Rbm^2)*e^(-Rbm/Lb)

# Total bulge integral
Nbulge = N(4*pi()*2*h(2600,70000))
Nbulge

8.83466118971905e11

# Calculate normalization factor for the integral
rho0 = N(3e11/(Ndisc+Nbulge))
rho0

0.0804810050755278

# Check normalization factor
N(rho0*(Ndisc+Nbulge))

3.00000000000000e11

# Calculate the cut-off factor