Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- rscalc[rm2m_] := rm2m + 2 m*(1 + Log[(rm2m)/(2 m)])
- xmax = 700; xmin = -700; delx = 0.03; imax = Floor[(xmax - xmin)/delx];
- rstar = xmax; m = 1.0; l = 2.0; rcour = 0.99; ntot = 3;
- If[rstar > 4 m, r = rstar, r = 2 m*E^((rstar/(2 m)) - 1)];
- rm2m = r - 2 m;
- Do[rstest = rscalc[rm2m];
- drsdr = (rm2m + 2 m)/rm2m;
- rtempm2m = rm2m + (rstar - rstest)/drsdr;
- rm2m = rtempm2m;, {i, 7}];
- mesh = NDSolve[{rm2mc'[rs] == rm2mc[rs]/(rm2mc[rs] + 2 m),
- rm2mc[rstar] == rm2m}, rm2mc, {rs, xmax, xmin}, AccuracyGoal -> 25,
- PrecisionGoal -> 25, WorkingPrecision -> 50]
- Plot[Evaluate[rm2mc[lex] /. mesh], {lex, -20, 20}, PlotRange -> All]
- xtest = -250;
- b = Evaluate[rm2mc[xtest] /. mesh]
- xchck = rscalc[b];
- err = xchck - xtest
- V[rm2m_] := (rm2m/(rm2m +
- 2 m))*((l (l + 1))/(rm2m + 2 m)^2 + (2 m/(rm2m + 2 m)^3));
- Plot[V[nax - 2 m], {nax, 0, 10}]
- Vsx[x_] := V[Evaluate[rm2mc[x] /. mesh]];
- Plot[Vsx[x], {x, -20, 20}]
- XTab = Table[xmin + i*delx, {i, imax + 1}];
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement