reading the numerical data and calculating everything in Weyl-Papapetrou coordinates

In[23]:=

pathread = "D:\mathematica notebook\NS-model" ;

SetDirectory[pathread] ;

Dimensions[eosFilesTab]

Count[%, 0]

eosFilesTabNew = Table[eosFilesTab[[i]], {i, 1, %}]

Out[26]=

{3}

Out[28]=

1

Out[29]=

{L-e_c0.3790-r-6.41528e-01.txt}

In[112]:=

<<NumericalMath`ListIntegrate`

In[1042]:=

overallstarttime = SessionTime[] ;

<<Statistics`LinearRegression`

rhoSeries = Series[2 Log[(1 - M/(2 r))/(1 + M/(2 r))] - Log[1 - M^2/(4 r^2)], {r, ∞, 6}] ;

rhoSeriesTab = Table[SeriesCoefficient[rhoSeries, i], {i, 0, 6}] ;

TableForm[Table[GegeP_n[z_] = FullSimplify[π/2^(1/2) ((-1)^n Gamma[n + 2])/(2^(n + 1/2) n ! Gamma[n + 3/2]) (1 - z^2)^(-1/2) D[(1 - z^2)^(n + 1/2), {z, n}]], {n, 0, 5}]] ;

assympt = 300 ;

DRiscoTabP = Table[{Null, Null, Null, Null, Null, Null, Null, Null}, {i, 1, Dimensions[eosFilesTabNew][[1]]}] ;

DRiscoTabM = Table[{Null, Null, Null, Null, Null, Null, Null, Null}, {i, 1, Dimensions[eosFilesTabNew][[1]]}] ;

DEtildaPlot = Table[Null, {i, 1, Dimensions[eosFilesTabNew][[1]]}] ;

modIn = 1 ;

starttime = SessionTime[] ;

SetDirectory[pathread] ;

dataFileName = StringDrop[eosFilesTabNew[[modIn]], -4] ;

Print[dataFileName] ;

ModelTypeTab[[modIn + 1, 1]] = dataFileName ;

fileToRead1 = dataFileName<>"-p9.txt" ;

fileToRead2 = dataFileName<>".txt" ;

tableMoments = Import[fileToRead1, "Table"] ;

tabM1 = Flatten[tableMoments[[Dimensions[tableMoments]]]] ;

tabM2 = Flatten[tableMoments[[4]]] ;

Print[Table[{tabM1[[i]], tabM2[[i]]}, {i, 1, Dimensions[tabM1][[1]]}]//MatrixForm] ; tableData = Import[fileToRead2, "Table"] ;

tableData[[Dimensions[tableMoments][[1]] + 21]] ;

tableData[[Dimensions[tableMoments][[1]] + 23]] ;

numtable = Table[tableData[[Dimensions[tableMoments][[1]] + 22 + i]], {i, 1, Dimensions[tableData][[1]] - Dimensions[tableMoments][[1]] - 22}] ;

numericaldata = numtable ;

Clear[a, M, a2, c] ;

Clear[numericaldata] ;

rhoOFm = Table[{coslist[[i]], rholist[[i]]}, {i, kN - 1 jpol + 1, kN - 0 jpol}] ; omegaOFm = Table[{coslist[[i]], omegalist[[i]]/29.979}, {i, kN - 1 jpol + 1, kN - 0 jpol}] ;

Nωm = Interpolation[omegaOFm] ; Nrhom = Interpolation[rhoOFm] ; (*=== === === === === === === === === === === === === === === === === === === === === === === === == *)

Print["=============== Correcting the moments =================="] ;

Clear[b] ;

BTab = Table[{1/rlist[[(i - 1) * jpol + j]], coslist[[(i - 1) * jpol + j]], Exp[gammalist[[(i - 1) * jpol + j]]]}, {i, 601 - assympt, 601}, {j, 1, 301}] ;

kmax = 2 ; vFit = Underoverscript[∑, k = 0, arg3] b_k x^k ; BCoeffG0 = FindFit[BTabIntG0, vFit, Table[b_k, {k, 0, kmax}], x, MaxIterations→100000] ;

Print[Coefficients of G (μ) for the function B: , BCoeffG0] ; 0

2 M Print[--- = , -M^2/4/.{M->Mgeo}, -> b ] ; 4 2

Print["b = ", b_2/M^2/.{M->Mgeo}/.BCoeffG0, " , b based on Laarakkers & Poisson should be -1/4"] ;

Print["============= Introducing the numerical solution on the equatorial plane ===================="] ;

j = 1 ;

Print["================ introducing the numerical solution on axis z ================================"] ;

Clear[j] ; (*=== === === === === === === === === === === === === === === === === === === === === === === === == *)

Nωρ = Interpolation[Table[{Nωlist[[i, 1]], {Nωlist[[i, 2]], NDωlist[[i - 302, 2]]}}, {i, 303, SDIV - 1}]] ;

Nfρ = Interpolation[Table[{Nflist[[i, 1]], {Nflist[[i, 2]], NDflist[[i - 302, 2]]}}, {i, 303, SDIV - 1}]] ;

Nωρ2 = Interpolation[Table[{Nωlist[[i, 1]], Nωlist[[i, 2]]}, {i, 301, SDIV}]] ;

Nfρ2 = Interpolation[Table[{Nflist[[i, 1]], Nflist[[i, 2]]}, {i, 301, SDIV}]] ;

Nωρ3 = Interpolation[Table[{Nωlist[[i, 1]], Nωlist[[i, 2]]}, {i, 1, SDIV}]] ;

Nfρ3 = Interpolation[Table[{Nflist[[i, 1]], Nflist[[i, 2]]}, {i, 1, SDIV}]] ;

Nωr = Interpolation[omegaOFr] ;

Nrhor = Interpolation[rhoOFr] ;

Print["======== Weyl coordinates of the numerical ISCOs and the eq. surface ========="] ;

rhoISCO2 = FindRoot[(Nfρ2[ρ]^(-1) ρ^2 - Nfρ2[ρ] Nωρ2[ρ]^2)^(1/2) == (Rs + hISCOminus), {ρ, NRs}]

rhoISCO1 = FindRoot[(Nfρ2[ρ]^(-1) ρ^2 - Nfρ2[ρ] Nωρ2[ρ]^2)^(1/2) == (Rs + hISCOplus), {ρ, NRs}]

Print["============== computing the frequencies ================"] ;

Out[1053]=

L-e_c0.3790-r-6.41528e-01

=============== Correcting the moments ==================

Coefficients of G (μ) for the function B: {b_0→0.999928, b_1→2.38862*10^-7, b_2→ -0.935529} 0

2 M --- = -1.10036 -> b 4 2

b = -0.212551 , b based on Laarakkers & Poisson should be -1/4

S corrected= -100.765 3

M corrected= -32.9295 2

a= 1.46615

M= 2.09796

M two-soliton= 631.657 4

============= Introducing the numerical solution on the equatorial plane ====================

================ introducing the numerical solution on axis z ================================

======== Weyl coordinates of the numerical ISCOs and the eq. surface =========

Out[1105]=

{ρ→18.8577}

Out[1106]=

{ρ→16.798}

Out[1107]=

{ρ→16.798}

============== computing the frequencies ================

total time = 0 h 1 min 10 sec

Created by Mathematica  (September 26, 2012) Valid XHTML 1.1!