This notebook should be run AFTER the intermediate files have been generated. This file only contains “part 1” (the first case) of the calculation. It should be edited to work on different parts of the calculation.

Do not forget to set your directory appropriately.

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

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


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));ταphys[Te_,n_]:= 1.17 (10^18)((Te^(3/2))/n);


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

Do[

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

rhs1[n_,nα_,Ti_,Te_,c2_,τEi_,Shi_,She_]:=Shi-1.5 ((n Ti)/(τEi))+1.5 ((n(Te-Ti))/(τei[Te,Ti,n]));
rhs2[n_,nα_,Ti_,Te_,c2_,τEi_,Shi_,She_]:=She-1.5 ((n Te)/(c2 τEi))+1.5 ((n(Ti-Te))/(τei[Te,Ti,n]))+((nα)/(cτα τEi ))Eα-cB (n^2) (Te^(1/2));
rhs3[n_,nα_,Ti_,Te_,c2_,c4_,τEi_,Shi_,She_]:=((n^2)/4)σvfun[Ti]-((nα)/(cτα τEi ))-((nα)/(c4 τEi ));

SetDirectory[];

SetDirectory["../intermediate files"];

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

reftab=Import[filename];

nT=Length[reftab];

nTdo=nT;
Sheloc=0;

ταref=N[ταphys[20,nloc]];
τeiref=N[τei[20,20,nloc]];


nmax=100;


vectemp={0,0,0,0,0,0,0,0,0,0,0,0,0};


summarytab=Table[{0,0,0,0,0,0,0,0,0,0,0,0,0},{iii,1,nTdo}];

SetDirectory[];

SetDirectory["../intermediate files/part 1"]

Do[

Tiref=reftab[[iT,1]];
nloc=reftab[[iT,3]];
τEiloc=reftab[[iT,5]];
c2loc=reftab[[iT,6]];
c4loc=reftab[[iT,7]];

τhat = reftab[[iT,4]];

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

completed=0;

Shiup=2. (10^22);
Shidown=1. (10^16);

dShi=(10.^(-3));

valup = 2.;
valdown=-0.4;
(* "random" values, they don't matter too much *)

Print["@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"];
Print["iT = ", iT, "   τhat = ",τhat, "   Tiref =", Tiref];
Print["@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"];


Do[

If[completed>0,Goto[skip]];

If[EvenQ[i+1],
Shiloc=((Shiup+Shidown)/2),

fp=((valup-valdown)/(Shiup-Shidown));
Shiloc=Shidown-(valdown/fp)
];

(*Print["i = ", i];
Print["Shiup = ", Shiup, "   Shidown = ", Shidown];
Print["Shiloc = ",Shiloc];*)

eqdiff1NL:=-1.5nloc  Ti'[t]+rhs1[nloc,nα[t],Ti[t],Te[t],c2loc,τEiloc,Shiloc,Sheloc];

eqdiff2NL:=-1.5 nloc  Te'[t]+rhs2[nloc,nα[t],Ti[t],Te[t],c2loc,τEiloc,Shiloc,Sheloc];

eqdiff3NL:=-nα'[t]+rhs3[nloc,nα[t],Ti[t],Te[t],c2loc,c4loc,τEiloc,Shiloc,Sheloc];

tmax=100 τEiloc;

solNL=NDSolve[{eqdiff1NL==0,eqdiff2NL==0,eqdiff3NL==0,Ti[0]==1.001,Te[0]==1.0011,nα[0]==0},{Ti[t],Te[t],nα[t]},{t,0,tmax}];
Ti1NL[t_]=First[Chop[Ti[t]/.solNL]];
Te1NL[t_]=First[Chop[Te[t]/.solNL]];
nα1NL[t_]=First[Chop[nα[t]/.solNL]];

Tifin=Ti1NL[tmax];
Tefin=Te1NL[tmax];

funval=Tifin/Tiref-1.;

heattab[[i]]={τEiloc,Shiloc,c2loc,c4loc,Tifin,ταref,τeiref,τhat,N[nloc Tifin],N[nloc Tefin],N[nloc Tifin]+N[nloc Tefin]};

If[funval>0,vectemp={τEiloc,Shiloc,c2loc,c4loc,Tifin,ταref,τeiref,τhat,N[nloc Tifin],N[nloc Tefin],N[nloc Tifin]+N[nloc Tefin],Tiref,τEiloc*Shiloc}];

(*Print["Tifin = ", Tifin, "    funval = ",funval];*)
If[funval*valdown<0,
Shiup=Shiloc;
valup=funval,
Shidown=Shiloc;
valdown=funval];


If[(Abs[((Shiup-Shidown)/Shiup)]<dShi),
completed=1;
Print["tau_Ei = ",τEiloc,"   ","Shiloc = ", Shiloc];
Print["Shiup = ", Shiup, "   ", "Shidown = ", Shidown];
Print["Tifin = ", Tifin]
];

Label[skip],

{i,1,nmax}];

filename="tauEi"<>ToString[τEiloc]<>"c_tau_alpha_"<>ToString[cτα]<>".dat";
Export[filename,heattab,"Table"];

MatrixForm[heattab];

summarytab[[iT]]=vectemp,

{iT,1,nTdo}
];

summarytabcase1=summarytab;

SetDirectory[];

SetDirectory["../output tables"]

Export["summary_case1.dat",summarytab,"Table"];

SetDirectory[];

SetDirectory["../output tables"];

fileoutname = "summary_part1_case_"<>ToString[cτα]<>".dat";
Export[fileoutname,summarytab,"Table"],
{iτα,1,nτα}

]

(*Plot[Ti1NL[t],{t,0,tmax}]
Plot[Ti1NL'[t],{t,0,tmax}]
Plot[nα1NL[t],{t,0,tmax}]*)