# Modeling a stellar atmosphere

I'm trying to model a stellar atmosphere using Mathematica. Basically I'm working in a physical and chemical model which combined should give me a profile of the chemical abundance in the atmosphere. The physical model is based in two equations:

t[r_,ϕ_]:= Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/r)^2])^(1/4);
dens[r_, ϕ_] = n*Exp[-((2*cg* mf* mh)/(kb*Te[ϕ]*Ra[ϕ]*AU))*(1-(Ra[ϕ]/r)^(1/2))];


which gives me both the temperature and density as a function of the radius. Given those statements I write the differential equations that will give the rate of formation and destruction of the chemical species: $\displaystyle \frac{d(AB)}{dt}=kn(A)n(B)$, where $k$ is the rate coefficient and it depends on the temperature $t$. Now, I have built the code (with your help) to solve the system of ODE and it is working pretty fine:

AU = 1.49597871*10^13;
To = 2200;
ΔT = 440;
Ro = 1.8;
ΔR = 0.1;
n = 0.5*10^15;
cg = 6.67384*10^-8;
kb = 1.3806488*10^-16;
mh = 1.0078250321*(1.660538921*10^-24);
mf = 1.98892*10^33;
Te[ϕ_]:= 2200 + 440*Cos[ϕ];
Ra[ϕ_]:= 2 + 0.25*Cos[ϕ + Pi];
t[r_,ϕ_]:= Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/r)^2])^(1/4);

k1[r_,ϕ_]:= (6.99*10^-14)*((t[r,ϕ]/300)^2.8)*Exp[-1950/t[r,ϕ]];
k2[r_,ϕ_]:= (1.59*10^-11)*((t[r,ϕ]/300)^1.2)*Exp[-9610/t[r,ϕ]];
k17[r_,ϕ_]:= (3.14*10^-13)*((t[r,ϕ]/300)^2.7)*Exp[-3150/t[r,ϕ]];
k18[r_,ϕ_]:= (2.05*10^-12)*((t[r,ϕ]/300)^1.52)*Exp[-1736/t[r,ϕ]];
k62[r_,ϕ_]:= (1.77*10^-11)*Exp[178/t[r,ϕ]];
k63[r_,ϕ_]:= (1.85*10^-11)*((t[r,ϕ]/300)^0.95)*Exp[-8571/t[r,ϕ]];
k94[r_,ϕ_] := (1.65*10^-12)*((t[r,ϕ]/300)^1.14)*Exp[-50/t[r,ϕ]];
k138[r_,ϕ_] := (5.94*10^-17)*((t[r,ϕ]/300)^0.17)*Exp[65.9/t[r,ϕ]];
k141[r_,ϕ_] := (1*10^-11)*Exp[-4800/t[r,ϕ]];
k143[r_,ϕ_] := (1*10^-11)*Exp[-4800/t[r,ϕ]];
k144[r_,ϕ_] := 1*10^-9;

r1 = k1[r,ϕ] nH[r] nOH[r];
r2 = k2[r,ϕ] nH[r] nH2O[r];
r17 = k17[r,ϕ] nH2[r] nO[r];
r18 = k18[r,ϕ] nH2[r] nOH[r];
r62 = k62[r,ϕ] nO[r] nOH[r];
r63 = k63[r,ϕ] nO[r] nH2O[r];
r94 = k94[r,ϕ] nOH[r] nOH[r];
r138 = k138[r,ϕ] nTi[r] nO[r];
r141 = k141[r,ϕ] nTi[r] nOH[r];
r143 = k143[r,ϕ] nTi[r] nH2O[r];
r144 = k144[r,ϕ] nTiO[r] nH2O[r];
eqns = {nH'[r] == (r17 + r18 + r62) - (r1 + r2),
nOH'[r] == (r2 + r17 + r63) - (r1 + r18 + r62 + r94 + r141),
nO'[r] == (r1 + r94) - (r62 + r63 + r138),
nH2'[r] == (r1 + r2) - (r17 + r18),
nH2O'[r] == (r18 + r94) - (r2 + r143 + r144),
nTi'[r] == -(r138 + r141 + r143 + r144),
nTiO'[r] == (r138 + r141 + r143) - (r144)};

eqEqn = {nH[r] + nOH[r] + nO[r] + nH2[r] + nH2O[r] + nTi[r] +
nTiO[r] + nConst[r] == 0.500029396};

ic = {nH[2] == 2*10^-7, nOH[2] == 9.8*10^-7,
nO[2] == 2*10^-5, nH2[2] == 5*10^-1,
nH2O[2] == 8*10^-6, nTi[2] == 1.96*10^-7,
nTiO[2] == 2*10^-8, nConst[2] == 0};

sol = Table[
NDSolve[{eqns, eqEqn, ic} /. ϕ -> ϕi, {nH, nOH, nO, nH2,
nH2O, nTi, nTiO}, {r, 2, 9}, MaxSteps -> Infinity,
AccuracyGoal -> 4, PrecisionGoal -> 4,
Method -> {"TimeIntegration" -> "StateSpace"}], {ϕi, 0,
2*Pi, Pi/5}];

solTable = DeleteCases[sol, {}];


But I'm missing one thing. I'm not seeing how could I add the density equation

dens[r_, ϕ_] = n*Exp[-((2*cg* mf* mh)/(kb*Te[ϕ]*Ra[ϕ]* AU))*(1-(Ra[ϕ]/r)^(1/2))];


to this code.

It should multiply the species abundances n(X) and then the result to be used to solve the ODE. So for example, when $r=2$, $nOH[2]=nOH[0]*dens[2]$. I tried to multiply the initial conditions to the density, but I get a error:

NDSolve::ndinnt: "Initial condition is not a number or a rectangular array of numbers."

Could anyone help me with that? What should I do to incorporate the density equation to solve this problem? Moreover, after the calculation, the final result should be given divided for the density. So, for example, $nOH[2]=nOH[0]*dens[2]$ is used to calculate the ODE, but the final result should be given by $nOH[2]=nOH[2]/dens$. Got it? Is there some way to do that? Thank you very much.

• Your code does not run as written, because equations such as r1 = k1 nH[r] nOH[r]; should be r1 = k1[r, \[Phi]] nH[r] nOH[r];. Even with this correction, introducing dens in the way you describe overspecifies your dependent variables. Commented Jun 22, 2015 at 5:09

As noted in the comment above, k1 etc must be called with its argument list. With this change the r1 etc equations become

r1 = k1[r, ϕ] dens[r, ϕ]^2 nH[r] nOH[r];
r2 = k2[r, ϕ] dens[r, ϕ]^2 nH[r] nH2O[r];
r17 = k17[r, ϕ] dens[r, ϕ]^2 nH2[r] nO[r];
r18 = k18[r, ϕ] dens[r, ϕ]^2  nH2[r] nOH[r];
r62 = k62[r, ϕ] dens[r, ϕ]^2  nO[r] nOH[r];
r63 = k63[r, ϕ] dens[r, ϕ]^2  nO[r] nH2O[r];
r94 = k94[r, ϕ] dens[r, ϕ]^2  nOH[r] nOH[r];
r138 = k138[r, ϕ] dens[r, ϕ]^2 nTi[r] nO[r];
r141 = k141[r, ϕ] dens[r, ϕ]^2  nTi[r] nOH[r];
r143 = k143[r, ϕ] dens[r, ϕ]^2 nTi[r] nH2O[r];
r144 = k144[r, ϕ] dens[r, ϕ]^2 nTiO[r] nH2O[r];


I interpret the density discussion in the question to mean that expressions of the form

D[species1[r],r] == c species2[r] species3[r]


should be replaced by expressions of the form

D[species1[r] dens[r], r] == c species2[r] dens[r] species3[r] dens[r]


which can be rewritten as

(D[species1[r], r] dens[r] + species1[r] D[dens[r], r])/dens[r]^2 ==
c species2[r] species3[r]


The operation which performs this transformation is

Table[eqns[[i, 1]] = eqns[[i, 1]] /. Derivative[1][z_][r] ->
((Derivative[1][z][r] + z[r] dldens[r, ϕ])/dens[r, ϕ], {i, 7}];


which yields

eqns = {(dldens[r, ϕ] nH[r] + nH'[r])/dens[r, ϕ] == -r1 + r17 + r18 - r2 + r62,
(dldens[r, ϕ] nOH[r] + nOH'[r])/dens[r, ϕ] == -r1 - r141 + r17 - r18
+ r2 - r62 + r63 - r94,
(dldens[r, ϕ] nO[r] + nO'[r])/dens[r, ϕ] == r1 - r138 - r62 - r63 + r94,
(dldens[r, ϕ] nH2[r] + nH2'[r])/dens[r, ϕ] == r1 - r17 - r18 + r2,
(dldens[r, ϕ] nH2O[r] + nH2O'[r])/dens[r, ϕ] == -r143 - r144 + r18 - r2 + r94,
(dldens[r, ϕ] nTi[r] + nTi'[r])/dens[r, ϕ] == -r138 - r141 - r143 - r144,
(dldens[r, ϕ] nTiO[r] + nTiO'[r])/dens[r, ϕ] == r138 + r141 + r143 - r144};


with

dens[r_, ϕ_] := n*Exp[-((2*cg*mf*mh)/(kb*Te[ϕ]*Ra[ϕ]*AU))*(1 - (Ra[ϕ]/r)^(1/2))];
dldens[r_, ϕ_] := D[dens[r, ϕ], r]/dens[r, ϕ];


Then,

sol = Table[NDSolveValue[{eqns, ic} /. ϕ -> ϕi, {nH, nOH, nO, nH2, nH2O, nTi, nTiO},
{r, 2, 9}, {ϕi, 0, 2*Pi, Pi/5}];


solves the equations without difficulty, and the first element of sol, for instance, can be plotted as

LogPlot[Evaluate[Table[sol[[1, i]][r], {i, 7}]], {r, 2, 9},
PlotLegends -> {nH, nOH, nO, nH2, nH2O, nTi, nTiO}]


Note that nConst[r] is omitted from the calculation for convenience. It can be computed at the end using Solve on eqEqn, if desired.

• Thank you very much @bbgodfrey! But oddly I get an error message when trying to run it: "NDSolveValue::ndsz: At r == 2., step size is effectively zero; singularity or stiff system suspected." I tried some adjustments in the code but it didn't work. Commented Jun 23, 2015 at 0:24
• My computation runs quickly and without error. Do the equations you feed to NDSolveValue look like (dldens[r, ϕ] nH[r] + nH]'[r])/ dens[r, ϕ] == -r1 + r17 + r18 - r2 + r62? I do know that replacing dens[r, ϕ] by dens[r, ϕ]^2 in these equations produces the error you describe. Commented Jun 23, 2015 at 0:33
• Yes, precisely like that. I may have misunderstood the operation Table[eqns[[i, 1]] = eqns[[i, 1]] /. Derivative[1][z_][r] -> ((Derivative[1][z][r] + z[r] dldens[r, ϕ])/dens[r, ϕ], {i, 7}];, to transform the expression because everything else looks exactly as you said. Commented Jun 23, 2015 at 0:53
• If your equations look like the one in my last comment after applying the Table` operation, I would expect that you would obtain the same results that I do. I shall add all seven equations to my answer. Commented Jun 23, 2015 at 1:13
• Thank you very much for your help @bbgodfrey! Commented Jun 23, 2015 at 1:24