Run this file first to generate intermediate files needed when calculating the heating power needed to reach ignition.

This notebook contains several different cases for convenience.

Do not forget to set your output directory.

σ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=100;

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

nloc=(10^20);
Tmin = 6.;
Tmax=40;
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};
soltabBnew[[i]]={Tloc,τEsol,nloc};
soltabbigB[[i]]={Tloc,τEsol, nloc Tloc τEsol/(10.^20),2nloc Tloc τEsol/(10.^20)},

{i,1,nmax}]

MatrixForm[soltabB];
ListPlot[soltabB]

Output part

SetDirectory[];

SetDirectory["../intermediate files"]

Export["reftab_SF.dat",soltabBnew,"Table"]

Two Fluid part

Part 1: Changing (Subscript[τ, α])

valarray={0.1,0.4,0.7,1,10};
nτα=Length[valarray];

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

c2loc=1;
c4loc=(1./small);

Do[

cτα=valarray[[iτα]];
τα[Te_,n_]:=cτα  τEi;


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,τα->cτα  τEi}];
eq0D2d=FullSimplify[ReplaceAll[eq0D2c,{τEi->First[τEi/.s2]}]];


eq0D2e=ReplaceAll[eq0D2d,{CB->cB,τei->τei[Te,Ti,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]}];


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

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


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}];

soltabTFBnewcase1=soltabTFBnew;

Print[ListPlot[soltabTFB]];

Print[ListPlot[{soltabB,soltabTFB},GridLines->Automatic,PlotRange->All]];

SetDirectory[];

SetDirectory["../intermediate files"];

filename="reftab_TF_part1_case_"<>ToString[cτα]<>".dat";

Export[filename,soltabTFBnew,"Table"],

{iτα,1,nτα}

]

Case 2 : 2 Values for τα, Physical (Subscript[τ, ei]), c2=0.1, large c4

small = (10^(-6));

valarray={small,1};
nτα=Length[valarray];

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


c2loc=.1;
c4loc=(1./small);

Do[

cτα=valarray[[iτα]];

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

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]}];


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

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


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}];

soltabTFBnewcase1=soltabTFBnew;

Print[ListPlot[soltabTFB]];

Print[ListPlot[{soltabB,soltabTFB},GridLines->Automatic,PlotRange->All]];

SetDirectory[];

SetDirectory["../intermediate files"];

filename="reftab_TF_part2_case_"<>ToString[cτα]<>".dat";

Export[filename,soltabTFBnew,"Table"],

{iτα,1,nτα}

]

Case 3 : Vary c2, large c4, All Rest Physical

small = (10^(-6));

valarray={0.1,0.4,0.7,1};
nc2=Length[valarray];

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


c4loc=(1./small);

Do[

c2loc=valarray[[ic2]];



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]}];


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

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


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}];

soltabTFBnewcase1=soltabTFBnew;

Print[ListPlot[soltabTFB]];

Print[ListPlot[{soltabB,soltabTFB},GridLines->Automatic,PlotRange->All]];

SetDirectory[];

SetDirectory["../intermediate files"];

filename="reftab_TF_part3_case_"<>ToString[c2loc]<>".dat";

Export[filename,soltabTFBnew,"Table"],

{ic2,1,nc2}

]

Case 4 : c2=1, Vary c4, All Rest Physical

small = (10^(-6));

valarray={0.1,0.4,0.7,1,10};
nc4=Length[valarray];

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

c2loc=1;


Do[

c4loc=valarray[[ic4]];


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]}];


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

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


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}];

soltabTFBnewcase1=soltabTFBnew;

Print[ListPlot[soltabTFB]];

Print[ListPlot[{soltabB,soltabTFB},GridLines->Automatic,PlotRange->All]];

SetDirectory[];

SetDirectory["../intermediate files"];

filename="reftab_TF_part4_case_"<>ToString[c4loc]<>".dat";

Export[filename,soltabTFBnew,"Table"],

{ic4,1,nc4}

]

Case 2

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

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]}]

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

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

c2loc=1;
c4loc=1;


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}]

soltabTFBnewcase2=soltabTFBnew;

ListPlot[soltabTFB]

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

Output part

SetDirectory[];

SetDirectory["../intermediate files"]

Export["reftab_TF_case2.dat",soltabTFBnew,"Table"]

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

Case 3

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

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]}]

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

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

c2loc=1;
c4loc=1;


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}]

soltabTFBnewcase3=soltabTFBnew;

ListPlot[soltabTFB]

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

Output part

SetDirectory[];

SetDirectory["../intermediate files"]

Export["reftab_TF_case3.dat",soltabTFBnew,"Table"]

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

Case 4

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

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]}]

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

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

c2loc=1;
c4loc=1;


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}]

soltabTFBnewcase4=soltabTFBnew;

ListPlot[soltabTFB]

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

Output part

SetDirectory[];

SetDirectory["../intermediate files"]

Export["reftab_TF_case4.dat",soltabTFBnew,"Table"]

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

Case 5

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

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]}]

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

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

c2loc=1;
c4loc=(10^6);


Do[

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

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

τEsol=eqτEi[nloc,Tiloc,Tesol];
(* Differently from the other version, this is τEi, so we need to calculate τhat separately *)

τEhat=2((c2loc ((τEsol)^2))/(c2loc τEsol+τEsol));

soltabTFBnew[[i]]={Tiloc,Tesol,nloc,τEhat,τEsol,c2loc,c4loc ,τei[Tesol,Tiloc,nloc]};

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

{i,1,nmax}]

soltabTFBnewcase5=soltabTFBnew;

ListPlot[soltabTFB]

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

Output part

SetDirectory[];

SetDirectory["../intermediate files"]

Export["reftab_TF_case5.dat",soltabTFBnew,"Table"]

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