Resolution of the TOV equations
Morgane Fortin & Jérôme Margueron, Orsay 2008
CONSTANTS :
kboltz = 8.617 * 10 ^ - 11 H* Boltzmann constant in MeV.K^-1 *L;
hbc = 197.3 H* MeV.fm *L;
mneut = 938.5 H* neutron mass in MeV *L;
H*cm,g,s units*L
msol = 1.989 * 10 ^ 33; H*mass of the Sun in g*L
G = 6.6720 * 10 ^ - 8;
c = 2.9979 * 10 ^ 10;
rsol = 2 * G * msol  Hc ^ 2L; H*Schwarzschild radius of the Sun, in cm*L
rhosol = msol * 3  H4 * Pi * rsol ^ 3L;
rhodrip = 5.409 * 10 ^ 11;
rhosat = 2.66 * 10 ^ 14  2.; H* expected saturation nuclear matter density*L
IMPORT THE EOS (SLY) FROM A TABLE :
obtained from http://www.ioffe.ru/astro/NSG/NSEOS/
tabs = Import@"sly.eos", 8"Table"<D;
prs@d_D :=
Interpolation@Transpose@8tabs@@All, 3DD, tabs@@All, 4DD<D, InterpolationOrder ® 2D@dD;
rhs@p_D := Interpolation@Transpose@8tabs@@All, 4DD, tabs@@All, 3DD<D,
InterpolationOrder ® 2D@pD;
LogLogPlotAprs@xD, 9x, 100, 1015 =, AxesLabel ® 8"density Hgcm^3L", "Pressure Hdyncm^2L"<E
Pressure Hdyncm^2L
1035
1031
1027
1023
1019
density Hgcm^3L
105
108
1011
1014
2
tov.nb
CALCULATION OF THE DENSITY PROFILE
H* defines the boudnary condition *L
rho0 = 7.6 * rhosat; p0 = prs@rho0D; rmin = 0.001; rmax = 5;
H* solves the TOV equations *L
solutions = NDSolve@
8m¢ @rD Š If@p@rD > 0.0, 3 * rhs@p@rD * p0D  rhosol * r ^ 2, 0.0D,
p¢ @rD * r * Hm@rD - rL Š
If@p@rD > 0.0, 0.5 * rhs@p@rD * p0D * c ^ 2  p0 * m@rD * H1.0 + p@rD * p0  Hrhs@p@rD * p0D * c ^ 2LL *
H1.0 + 3 * p0  Hrhosol * c ^ 2L * p@rD * r ^ 3  m@rDL, 0.0D,
p@rminD Š 1.0, m@rminD Š rho0 * 4.0  3.0 * Pi * rmin ^ 3  msol<,
8m, p<,
8r, rmin, rmax<, Method -> 8"ExplicitRungeKutta", "DifferenceOrder" -> 8<D;
H* store the results *L
mass@r_D := m@rD . solutions@@1DD;
pression@r_D := p@rD . solutions@@1DD;
density@r_D := rhs@p@rD * p0D . solutions@@1DD;
H* calculation of the total radius: *L
rts = 1; rstep = 0.001; While@pression@rtsD > 0, rts = rts + rstepD; rts = rts - rstep;
H* calculation of the core radius *L
rcs = 1; rstep = 0.001;
While@rhs@pression@rcsD * p0D > rhosat  2., rcs = rcs + rstepD; rcs = rcs - rstep;
H* convert in ordinary units *L
oupress@r_D := p0 * pression@r  rsol * 10 ^ 5D;
oumass@r_D := mass@r  rsol * 10 ^ 5D;
oudens@r_D := density@r  rsol * 10 ^ 5D;
H* print the results *L
Print@"Mass = ", mass@rtsD, " sol. mass; Radius = ", rts * rsol  10 ^ 5, " km; "D
Print@"Core radius = ", rcs * rsol  10 ^ 5, "km; ",
"Crust thickness = ", Hrts - rcsL * rsol  10 ^ 5, "km;"D
LogPlot@oupress@rD, 8r, 0.1, rts * rsol  10 ^ 5<,
AxesLabel ® 8"radius HkmL", "Pressure Hdyncm^2L"<, PlotStyle ® 8Blue<D
LogPlot@oudens@rD, 8r, 0.1, rts * rsol  10 ^ 5<,
AxesLabel ® 8"radius HkmL", "density Hgcm^3L"<, PlotStyle ® 8Blue<D
Plot@oumass@rD, 8r, 0.1, rts * rsol  10 ^ 5<,
AxesLabel ® 8"radius HkmL", "Mass Hsolar massL"<, PlotStyle ® 8Blue<D
Mass = 1.43146 sol. mass; Radius = 11.603 km;
Core radius = 10.909km;
Crust thickness = 0.693992km;
Pressure Hdyncm^2L
1 ´ 1035
5 ´ 1034
1 ´ 1034
5 ´ 1033
1 ´ 1033
5 ´ 1032
2
4
6
8
10
radius HkmL
tov.nb
density Hgcm^3L
1.0 ´ 1015
7.0 ´ 1014
5.0 ´ 1014
3.0 ´ 1014
2.0 ´ 1014
1.5 ´ 1014
2
4
6
8
10
6
8
10
radius HkmL
Mass Hsolar massL
1.4
1.2
1.0
0.8
0.6
0.4
0.2
2
4
radius HkmL
3
4
tov.nb
Calculate the M HRL relation
imax = 22;
Print@"central density,
mass,
radius"D;
Print@"
gcm3,
sol. mass,
km"D;
For@i = 1, i < imax, i ++,
rho0 = rhosat * H1.0 + i  5.01L ^ 2; H* define the central density *L
p0 = prs@rho0D; rmin = 0.001; rmax = 5;
H* solves the TOV equations *L
solutions = NDSolve@
8m¢ @rD Š If@p@rD > 0.0, 3 * rhs@p@rD * p0D  rhosol * r ^ 2, 0.0D,
p¢ @rD * r * Hm@rD - rL Š
If@p@rD > 0.0, 0.5 * rhs@p@rD * p0D * c ^ 2  p0 * m@rD * H1.0 + p@rD * p0  Hrhs@p@rD * p0D * c ^ 2LL *
H1.0 + 3 * p0  Hrhosol * c ^ 2L * p@rD * r ^ 3  m@rDL, 0.0D,
p@rminD Š 1.0, m@rminD Š rho0 * 4.0  3.0 * Pi * rmin ^ 3  msol<,
8m, p<,
8r, rmin, rmax<, Method -> 8"ExplicitRungeKutta", "DifferenceOrder" -> 8<D;
H* store the results *L
mass@r_D := m@rD . solutions@@1DD;
pression@r_D := p@rD . solutions@@1DD;
density@r_D := rhs@p@rD * p0D . solutions@@1DD;
H* calculation of the total radius: *L
rts = 1; rstep = 0.001; While@pression@rtsD > 0, rts = rts + rstepD; rts = rts - rstep;
Print@rho0, "
", mass@rtsD, "
", rts * rsol  10 ^ 5DD
tov.nb
central density,
gcm3,
1.91393 ´ 1014
2.60383 ´ 10
14
3.3997 ´ 1014
mass,
radius
sol. mass,
0.0971632
0.146068
0.244564
km
18.2062
16.2867
13.7794
4.30156 ´ 1014
0.397662
12.2881
5.30939 ´ 1014
0.596195
11.8303
6.42319 ´ 1014
0.822207
11.7565
14
1.05433
11.727
8.96873 ´ 1014
1.27373
11.6591
1.04005 ´ 1015
1.46726
11.5941
1.19382 ´ 1015
1.62831
11.4405
1.35818 ´ 1015
1.75634
11.2752
15
1.85452
11.095
1.71871 ´ 1015
1.92705
10.909
1.91487 ´ 1015
1.97833
10.7141
2.12163 ´ 1015
2.01286
10.5103
2.33899 ´ 1015
2.03434
10.3213
15
2.04595
10.153
2.80549 ´ 1015
2.05006
9.97576
3.05464 ´ 1015
2.0486
9.82811
3.31439 ´ 1015
2.04302
9.6834
3.58473 ´ 1015
2.03447
9.5446
7.64297 ´ 10
1.53315 ´ 10
2.56694 ´ 10
5

rooster CSPE MC

solidNotes

A Mapping Method for VISSR Date

Program : UTM Conversions (excel)