This notebook calculates the Lawson product for ignition as a function of temperature for the single-fluid and two-fluid models.

Do not forget to set the output directory.

For the TF model, c2=τEe/τEi, c4=τEα/τEi.

The SF limit is recovered setting c2=1, c4>>1, τα<<τEi, τei<<τEi.

τα[Te_,n_]:=1.17 (10^18)((Te^(3/2))/n);

τei[Te_,Ti_,n_]:=(((78 ((n (10^(-20)))/(((((Te+Ti)/2)))^(3/2)))))^(-1));

σvJF[T_]=Exp[k0σ+k1σ Log[T]+k2σ ((Log[T])^2)+k3σ ((Log[T])^3)+k4σ ((Log[T])^4)];
k0σ=-60.4593;
k1σ=6.1371;
k2σ=-0.8609;
k3σ=0.0356;
k4σ=-0.0045;

σvfun[x_]=σvJF[x];

Eα=3.5 (10^3);

kB=0.052;
eV=1.6022 (10^(-19));
cB = ((10 )^(-1))eV kB;

Single - Fluid Part

eqSFB[n_,T_,τE_]=n T τE-((12(T^2))/(Eα σvfun[T]-4cB (T^(1/2))))

nmax=1000;

soltabB=Table[{0,0},{i,1,nmax}];soltabbigB=Table[{0,0,0,0},{i,1,nmax}];

nloc=(10^20);
Tmin = 6.;
Tmax=100;
dT=((Tmax-Tmin)/(nmax-1.));

Do[

Tloc = Tmin+(i-1.)dT;

solloc = FindRoot[eqSFB[nloc,Tloc,τEloc],{τEloc,2,5}];
τEsol=τEloc/.solloc;
soltabB[[i]]={Tloc,τEsol};
soltabbigB[[i]]={Tloc,τEsol, nloc Tloc τEsol/(10.^20),2nloc Tloc τEsol/(10.^20)},

{i,1,nmax}]

MatrixForm[soltabB]; (* remove the ";" to see a table with the output *)
ListPlot[soltabB]

Two - Fluid Part

eq0D1=-(Ti/(τEi))+((Te-Ti)/(τei))
eq0D2=-(3/2)((n Te)/(τEe))+((nα Eα)/(τα))-CB (n^2)√Te -(3/2)((n(Te-Ti))/(τei))
eq0D3=((n^2)/4)σv -((nα)/(τα))-((nα)/(τEα))

s1=Solve[eq0D3==0,nα]
s2=Solve[eq0D1==0,τEi]

eq0D2b=ReplaceAll[eq0D2,{nα->First[nα/.s1],τEi->First[τEi/.s2]}]

eq0D2c=ReplaceAll[eq0D2b,{τEe ->c2 τEi,τEα ->c4 τEi}]
eq0D2d=FullSimplify[ReplaceAll[eq0D2c,{τEi->First[τEi/.s2]}]]

eq0D2e=ReplaceAll[eq0D2d,{CB->cB,τei->τei[Te,Ti,n],τα->τα[Te,n],σv->σvfun[Ti]}]

eq0Dlast[n_,Ti_,Te_,c2_,c4_]=eq0D2e
eqτEi[n_,Ti_,Te_]=ReplaceAll[First[τEi/.s2],{τei->τei[Te,Ti,n]}]

nmax=1000;

soltabTFB=Table[{0,0},{i,1,nmax}];soltabTFBtrip=Table[{0,0},{i,1,nmax}];soltabTFbigB=Table[{0,0,0,0},{i,1,nmax}];

nloc=(10^20);
Timin = 6.;
Timax=100;
dTi=((Timax-Timin)/(nmax-1.));

(* These set the values of c2 and c4 *)
c2loc=1.;
c4loc=1.;

Do[

Tiloc = Timin+(i-1.)dTi;

sollocTFB = FindRoot[eq0Dlast[nloc,Tiloc,Te,c2loc,c4loc]==0,{Te,Tiloc+.1}];
Tesol=Te/.sollocTFB;

τEsol=eqτEi[nloc,Tiloc,Tesol];

soltabTFB[[i]]={Tiloc,τEsol};
soltabTFBtrip[[i]]={Tiloc,nloc (Tiloc +Tesol)τEsol/(10.^20)};
soltabTFbigB[[i]]={Tiloc,τEsol,nloc Tiloc τEsol/(10.^20),nloc (Tiloc +Tesol)τEsol/(10.^20)},

{i,1,nmax}]

ListPlot[{soltabB,soltabTFB},GridLines->Automatic]
ListPlot[soltabTFBtrip]

MatrixForm[soltabTFBtrip];

MatrixForm[soltabTFbigB];

Output Part

(* Set directory here *)

SetDirectory[];

SetDirectory["/Users/Luca/Documents/lavoro/TF ignition/minimum Sh search/extra notebooks for paper/notebooks for webpage"]

Export["nTtau_SF_with_Brem.dat",soltabbigB,"Table"]
Export["nTtau_TF_with_Brem.dat",soltabTFbigB,"Table"]

FileNames[]