This notebook calculates the minimum Lawson product for ignition for a given profile and set of constants (e.g., τEα/τEi). The search algorithm for the minimum is not very advanced and may fail or occasionally converge slowly.
The notebook is set up to have η=2 (see paper), but it will work for different values if so desired.

NotebookAutoSave->True;

(* First define the variables proper of this notebook, then proceed with the usual cycle *)

SetDirectory[]

SetDirectory[""]

(* enter your directory here *)

nvalues=30;
datain=Table[{0,0,0,0,0,0,0,0,0,0,0},{i,1,nvalues}];

nave0=8.72273`*^19;
(* for consistency *)

ηmin=0.1;
ηmax=2;
dη=((ηmax-ηmin)/(nvalues-1.));

c2loc=1;
τref=6.;
Tref=60.;

Do[
ηloc=ηmin+dη(iloc-1.);
datain[[iloc]]={0.5, 1.5, ηloc,2.,τref,c2loc τref,Tref,0,0,0,0},
{iloc,1,nvalues}]

(* The first two values set the (fixed) temperature profile. use 0.5 and 1.5 for H-mode, 1.5 and 2.5 for L-mode. *)

Now proceed with the usual stuff

use Jeff Freidberg's function for <σv>
also use physical expression for (Subscript[τ, α]), taken from Betti ad Freidberg, PoF B 4, 1465 (1992)

note that we have defined:
(Subscript[τ, eq])=c1 τEi;
	(Subscript[τ, Ee])= c2 τEi;
	(Subscript[τ, α loss])=c4 τEi.

Eα=3.5 (10^3); (* [keV] *)

(* from Freidberg *)
(* n in (m^(-3)), T in keV *)
kB=0.052;
eV=1.6022 (10^(-19));
cB = ((10 )^(-1))eV kB;

σvJF[T_?NumericQ]=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_?NumericQ]=σvJF[x];


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

(* We assign input profiles for density and temperatures, so that we can study the effect of the profile on ignition. Profiles are kept fixed in time, i.e. they are picked once and for all at the beginning of each simulation. Ion and electron temperature profiles do not need to be the same.
To avoid numerical issues, functions do not go to 0 at the edge, but to a finite value tedge or nedge. *)

tedge = (10.^(-5));
nedge = (10.^(-5));

tempμν[r_?NumericQ,μ_?NumericQ,ν_?NumericQ]:=(1-tedge)(((1-(r^ν)))^μ)+tedge (* this is Th *)
nηθ[r_?NumericQ,η_?NumericQ,θ_?NumericQ]:=(1-nedge)(((1-(r^θ)))^η)+nedge
(* this is nh *)

(* Something needs to be made up for nα. We can use an arbitrary profile, or try and make one up that's derived from the system equations. We first try getting the nα profile from the α evolution equation. Observe that we solve for nαhat, we still take time-derivatives of nα[t] *)

(* Define an intermediate function and use funny variable names just in case; lval stands for local value. *)


fnα[r_?NumericQ,ηη_?NumericQ,θθ_?NumericQ,μμi_?NumericQ,ννi_?NumericQ,μμe_?NumericQ,ννe_?NumericQ,nlval_?NumericQ,Tilval_?NumericQ,Telval_?NumericQ,τEαlval_?NumericQ]:=(1/4)(((Re[nlval] nηθ[r,ηη,θθ]))^2)σvfun[Re[Tilval] tempμν[r,μμi,ννi]]/((1/(τα[Re[Telval] tempμν[r,μμe,ννe],Re[nlval] nηθ[r,ηη,θθ]]))+(1/(τEαlval)));


nαh[r_?NumericQ,ηη_?NumericQ,θθ_?NumericQ,μμi_?NumericQ,ννi_?NumericQ,μμe_?NumericQ,ννe_?NumericQ,nlval_?NumericQ,Tilval_?NumericQ,Telval_?NumericQ,τEαlval_?NumericQ]:=fnα[r,ηη,θθ,μμi,ννi,μμe,ννe,nlval,Tilval,Telval,τEαlval]/fnα[0,ηη,θθ,μμi,ννi,μμe,ννe,nlval,Tilval,Telval,τEαlval]

Notice that integrals are over the area, so in r dr dθ (here θ is the angular variable). For convenience and since nothing depends on the angle, we’ll do an integral in r dr only.
This also means that we are neglecting the effect of curvature (can put those back in later, if so desired).

Now we need to define the RHS of each differential equation (after integration in r dr). It seems reasonable to divide the nα equation by nαhat before integrating. This leaves only one messy integral.
Contrary to the original SF file, we leave the profile parameters in the definitions, so that if necessary one can iterate over parameters.

f1b[μ_?NumericQ,ν_?NumericQ,η_?NumericQ,θ_?NumericQ]:= NIntegrate[nηθ[r,η,θ]tempμν[r,μ,ν]r,{r,0,1}];
f2b[μ_?NumericQ,ν_?NumericQ,η_?NumericQ,θ_?NumericQ]:= NIntegrate[((nηθ[r,η,θ])^2)((tempμν[r,μ,ν])^(1/2))r,{r,0,1}]
f3b[μ_?NumericQ,ν_?NumericQ,η_?NumericQ,θ_?NumericQ]:= NIntegrate[(1/(τα[tempμν[r,μ,ν],nηθ[r,η,θ]]))r,{r,0,1}]



(* ion energy equation *)
rhs1b[n_?NumericQ,nα_?NumericQ,Ti_?NumericQ,Te_?NumericQ,c2_?NumericQ,τEi_?NumericQ,Shi_?NumericQ,She_?NumericQ,η_?NumericQ,θ_?NumericQ,μi_?NumericQ,νi_?NumericQ,μe_?NumericQ,νe_?NumericQ]:= Shi-1.5 ((n Ti)/(τEi))f1b[μi,νi,η,θ]+1.5 NIntegrate[ FullSimplify[((n nηθ[r,η,θ](Te tempμν[r,μe,νe]-Ti tempμν[r,μi,νi]))/(τei[Te tempμν[r,μe,νe],Ti  tempμν[r,μi,νi],n nηθ[r,η,θ]]))r],{r,0,1}]

(* electron energy equation *)
rhs2b[n_?NumericQ,nα_?NumericQ,Ti_?NumericQ,Te_?NumericQ,c2_?NumericQ,c4_?NumericQ,τEi_?NumericQ,Shi_?NumericQ,She_?NumericQ,η_?NumericQ,θ_?NumericQ,μi_?NumericQ,νi_?NumericQ,μe_?NumericQ,νe_?NumericQ]:=She-1.5 ((n Te)/(c2 τEi))f1b[μe,νe,η,θ]-1.5 NIntegrate[ FullSimplify[((n nηθ[r,η,θ](Te tempμν[r,μe,νe]-Ti tempμν[r,μi,νi]))/(τei[Te tempμν[r,μe,νe],Ti  tempμν[r,μi,νi],n nηθ[r,η,θ]]))r],{r,0,1}]+nα NIntegrate[(1/(τα[Te  tempμν[r,μe,νe],n nηθ[r,η,θ]] ))nαh[r,η,θ,μi,νi,μe,νe,n,Ti,Te,c4 τEi]r,{r,0,1}]Eα-cB (n^2) (Te^(1/2)) f2b[μe,νe,η,θ]

(* α density equation - remember: this is divided by nαh *)
rhs3b[n_?NumericQ,nα_?NumericQ,Ti_?NumericQ,Te_?NumericQ,c4_?NumericQ,τEi_?NumericQ,η_?NumericQ,θ_?NumericQ,μi_?NumericQ,νi_?NumericQ,μe_?NumericQ,νe_?NumericQ]:=((n^2)/4) NIntegrate[((((nηθ[r,η,θ])^2)σvfun[Ti tempμν[r,μi,νi]])/(nαh[r,η,θ,μi,νi,μe,νe,n,Ti,Te,c4 τEi]))r,{r,0,1}]-nα  (n/(Te^(3/2))) f3b[μe,νe,η,θ]-(1/2)((nα)/(c4 τEi ))

valuelist0={3.,180.,17,16.,18.,15.,19.,14.,20.,13.,21.,12.,22.,11.,23.,10.,24.,9.,25.,26.,27.};

nmax0=Length[valuelist0];

(* use a list of values for bracketing;
if no brackets are found within the list, either there are no solutions or they are closer than 0.25 keV *)

μloc=datain[[1,1]];
νloc=datain[[1,2]];
(* these values are constant from input, so they can be assigned once and for all at the beginning *)
μetest=μloc;
μitest=μloc;

νetest=νloc;
νitest=νloc;


(* Notice that these values are not saved, so be sure you use the right ones! *)

c2loc=1;
c4loc=1;
(*@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@;
BE CAREFUL AND ENTER THE CORRECT VALUES FOR c2 and c4;
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@*)
nloc= (10^20);
Sheloc=0;
Shiloc=0;

nmax=100;

dtriplesearch=(10.^(-2)); (* RELATIVE distance between roots *)
dτsearch=1. (10.^(-2));
dT=.0125; (* distance between guesses in root finder *)

(*
Here we try to find the low-temperature root of the energy-balance that gives the minimum triple product for ignition.
We know that the system ignites for the tau_E and parameters in the input.
We assign a value to Ti, then evolve the system.
We know that if the system is at the ignition point Ti' is going to be >~0, with Te>Ti.
We also want to find the low root close to the minimum (where there is only one root), so we adjust τE until the two roots are closer than the prescribed distance.
*)



Do[


dτEisearch = 0.4;
sτ=1.;

ηloc=datain[[ivalues,3]];
θloc=datain[[ivalues,4]];
τEiloc=datain[[ivalues,5]]-dτEisearch;
(*τEiloc=4.-dτEisearch;*)

Tiup=datain[[ivalues,7]];
Tidown=5.;
(* Tiup and Tidown are bracketing the "low" ignition root *)

ηtest=ηloc;
θtest=θloc;

nloc = nave0 / (2  NIntegrate[nηθ[r,ηtest,θtest]r,{r,0,1}]);

Print["@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"];
Print["ivalues = ", ivalues];
Print[μitest,"   ",νitest,"   ",ηtest,"   ",θtest];
Print["@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"];

nave = 2 nloc NIntegrate[nηθ[r,ηtest,θtest]r,{r,0,1}];




ipass = 0;

(* since we are looking for a minimum, we need three points*)
τtab={0,0,0};
tripletab={0,0,0};

Label[repeatlabel];

If[ipass==0,
valuelist=valuelist0,
valuelist=valuelist1
];


nmax1=Length[valuelist];

TiL={0,0};
TiR={0,0};
TiprimL={0,0};
TiprimR={0,0};
Tiprim1=0;
Tiprim2 = 0;
(* saving Ti' twice, it looks convenient, but may change later *)

If[ipass<3,
τEiloc=τEiloc+sτ dτEisearch;
];


Print["_______________________________________________________"];
Print["τEiloc = ", τEiloc];
Print["_______________________________________________________"];

completed=0;

(* Define and fill matrices here, better in case of repeat *)
heattab=Table[{0,0,0,0},{iii,1,nmax+2}];
heattabbig=Table[{0,0,0,0,0,0,0,0,0,0,0,0,0, 0},{iii,1,nmax+2}];
heattabhuge=Table[{0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0},{iii,1,nmax+2}];
(* VERY IMPORTANT: THE FIRST ROW OF HEATTAB CONTAINS THE PROFILE VARIABLE VALUES *)
(* heattabhuge has an extra variable, Ti' *)

heattab[[1]]={μitest,νitest,ηtest,θtest};
heattabbig[[1]]={μitest,νitest,ηtest,θtest,0,0,0,0,0,0,0,0,0, 0};
heattabhuge[[1]]={μitest,νitest,ηtest,θtest,τEiloc,0,0,0,0,0,0,0,0, 0, 0};



(* First Do loop for bracketing *)

Do[

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

Tinow=valuelist[[i]];

Tifix[t_]=Tinow;

Print["i = ", i];
Print["Tinow = ",Tinow];

(* replace Ti with Tifix everywhere *)

(* Throw out the ion energy equation. Just calculate Ti' a posteriori. *)
(*
eqdiff1NL:=-1.5f1b[μitest,νitest,ηtest,θtest]nloc  Ti'[t]+rhs1b[nloc,nα[t],Ti[t],Te[t],c2loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest];
*)

eqdiff2NL:=-1.5f1b[μetest,νetest,ηtest,θtest] nloc  Te'[t]+rhs2b[nloc,nα[t],Tifix[t],Te[t],c2loc,c4loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest];

eqdiff3NL:=-(1/2)nα'[t]+rhs3b[nloc,nα[t],Tifix[t],Te[t],c4loc,τEiloc,ηtest,θtest,μitest,νitest,μetest,νetest];


tmax=30 τEiloc;

solNL=NDSolve[{(*eqdiff1NL==0,*)eqdiff2NL==0,eqdiff3NL==0,(*Ti[0]==100.002,*)Te[0]==100.001,nα[0]==0,WhenEvent[Abs[Te[t]]<Tinow,"StopIntegration"]},{(*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]];

tt=Te1NL["Domain"];
tdomain=Last[Last[tt]];

(*Tifin=Re[Ti1NL[tdomain]];*)
Tefin=Re[Te1NL[tdomain]];
nαfin=Re[nα1NL[tdomain]];
Tifin = Tinow;

If[Tefin<0,Tefin=0];

Tiprim=rhs1b[nloc,nαfin,Tifin,Tefin,c2loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest]/(1.5f1b[μitest,νitest,ηtest,θtest]nloc );

If[i==1, 
TiL[[1]]=Tinow;
TiprimL[[1]]=Tiprim;
Tiprim1=Tiprim];
If[i==2,
TiR[[2]]=Tinow;
TiprimR[[2]]=Tiprim;
 Tiprim2=Tiprim];

If[Tiprim1 Tiprim2<0,
Print["***********************************************"];
Print["***********************************************"];
Print["***********************************************"];
Print["***********************************************"];
Print["***********************************************"];
Print["SOMETHING WRONG: FIRST TWO VALUES BRACKET A ROOT!!"];
Print["***********************************************"];
Print["***********************************************"];
Print["***********************************************"];
Print["***********************************************"];
Print["***********************************************"];
Print["***********************************************"];
];

(* exit condition: root is bracketed *)
If[i>2,
If[Tiprim Tiprim1<0,
completed=1;
TiL[[2]]=Tinow;
TiprimL[[2]]=Tiprim;
TiR[[1]]=Tinow;
TiprimR[[1]]=Tiprim;
j=i;
]
];


piave= 2Tifin nloc NIntegrate[nηθ[r,ηtest,θtest]tempμν[r,μitest,νitest]r,{r,0,1}];
peave=2 Tefin nloc NIntegrate[nηθ[r,ηtest,θtest]tempμν[r,μetest,νetest]r,{r,0,1}];
ptot=piave+peave;
(* To keep things clean, do 
NOT multiply by eV. Also include the factor of 2. *)

power = Eα eV((nloc^2)/4) NIntegrate[((nηθ[r,ηtest,θtest])^2)σvfun[Tifin tempμν[r,μitest,νitest]]r,{r,0,1}];


heattabbig[[i+1]]={μitest,νitest,ηtest,θtest,N[τEiloc],c2loc*N[τEiloc],Tifin,Tefin,nαfin, nloc/nave,power, piave,peave,ptot*N[τEiloc]};


heattabhuge[[i+1]]={μitest,νitest,ηtest,θtest,N[τEiloc],c2loc*N[τEiloc],Tifin,Tefin,nαfin, nloc/nave,power, piave,peave,ptot*N[τEiloc],Tiprim};

heattab[[i+1]]={nloc, nave,N[τEiloc],Tifin};




Print["Tiprim = ", Tiprim];


Label[skipbrack],


{i,1,nmax1}];

If[completed==0,
If[sτ<0,
dτEisearch=dτEisearch/2
];
sτ=1;
Goto[repeatlabel]
];




(* Second Do loop for minimum finding *)
(* Now this needs to be repeated once per root (low and high). *)


Tiroot={0,0};
tripleroot={0,0};

Do[

Tiup=TiR[[k]];
Tidown=TiL[[k]];
Tiprimup=TiprimR[[k]];
Tiprimdown=TiprimL[[k]];
completed=0;

Do[

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

If[EvenQ[i],
Tinow=((Tiup+Tidown)/2),

fp=((Tiprimup-Tiprimdown)/(Tiup-Tidown));
Tinow=Tidown-(Tiprimdown/fp)
];

Print["i = ", i];
Print["Tiup = ", Tiup, "   Tidown = ", Tidown];
Print["Tinow = ",Tinow];

Tifix[t_]=Tinow;

(* replace Ti with Tifix everywhere *)

(* Throw out the ion energy equation. Just calculate Ti' a posteriori. *)
(*
eqdiff1NL:=-1.5f1b[μitest,νitest,ηtest,θtest]nloc  Ti'[t]+rhs1b[nloc,nα[t],Ti[t],Te[t],c2loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest];
*)

eqdiff2NL:=-1.5f1b[μetest,νetest,ηtest,θtest] nloc  Te'[t]+rhs2b[nloc,nα[t],Tifix[t],Te[t],c2loc,c4loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest];

eqdiff3NL:=-nα'[t]+rhs3b[nloc,nα[t],Tifix[t],Te[t],c4loc,τEiloc,ηtest,θtest,μitest,νitest,μetest,νetest];


tmax=30 τEiloc;

solNL=NDSolve[{(*eqdiff1NL==0,*)eqdiff2NL==0,eqdiff3NL==0,(*Ti[0]==100.002,*)Te[0]==100.001,nα[0]==0,WhenEvent[Abs[Te[t]]<Tinow,"StopIntegration"]},{(*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]];

tt=Te1NL["Domain"];
tdomain=Last[Last[tt]];

(*Tifin=Re[Ti1NL[tdomain]];*)
Tefin=Re[Te1NL[tdomain]];
nαfin=Re[nα1NL[tdomain]];
Tifin = Tinow;

If[Tefin<0,Tefin=0];

Tiprim=rhs1b[nloc,nαfin,Tifin,Tefin,c2loc,τEiloc,Shiloc,Sheloc,ηtest,θtest,μitest,νitest,μetest,νetest]/(1.5f1b[μitest,νitest,ηtest,θtest]nloc );

If[i==1,Tiprimref=Tiprim; Tiprimup=Tiprim];
If[i==2, Tiprimdown=Tiprim];


If[(Abs[Tiup-Tidown]<dT)∨(Abs[Tiprim]<(10.^(-4))),
completed=1;
Print["tau_Ei = ",τEiloc,"   ","Tifin = ", Tifin]
];


piave= 2Tifin nloc NIntegrate[nηθ[r,ηtest,θtest]tempμν[r,μitest,νitest]r,{r,0,1}];
peave=2 Tefin nloc NIntegrate[nηθ[r,ηtest,θtest]tempμν[r,μetest,νetest]r,{r,0,1}];
ptot=piave+peave;
(* To keep things clean, do 
NOT multiply by eV. Also include the factor of 2. *)

power = Eα eV((nloc^2)/4) NIntegrate[((nηθ[r,ηtest,θtest])^2)σvfun[Tifin tempμν[r,μitest,νitest]]r,{r,0,1}];


heattabbig[[i+1]]={μitest,νitest,ηtest,θtest,N[τEiloc],c2loc*N[τEiloc],Tifin,Tefin,nαfin, nloc/nave,power, piave,peave,ptot*N[τEiloc]};


heattabhuge[[i+1]]={μitest,νitest,ηtest,θtest,N[τEiloc],c2loc*N[τEiloc],Tifin,Tefin,nαfin, nloc/nave,power, piave,peave,ptot*N[τEiloc],Tiprim};

heattab[[i+1]]={nloc, nave,N[τEiloc],Tifin};


If[Tiprim*Tiprimdown<0,
Tiup=Tinow;
Tiprimup=Tiprim,
Tidown=Tinow;
Tiprimdown=Tiprim];

Print["Tiprim = ", Tiprim];

Label[skip],

{i,j+1,nmax}];

Tiroot[[k]]=Tinow;
tripleroot[[k]]=ptot*N[τEiloc];

heattab[[nmax+2]]={nloc, nave,datain[[ivalues,5]],datain[[ivalues,6]]};
heattabbig[[nmax+2]]={nloc, nave,datain[[ivalues,5]],datain[[ivalues,6]],0,0,0,0,0, 0,0, 0,0,0};


nustr=ToString[νitest];
mustr=ToString[μitest];
filename="igntab_mu_"<>mustr<>"_nu_"<>nustr<>"_eta_"<>ToString[ηtest]<>"_theta_"<>ToString[θtest]<>".dat";
bigfilename = "big"<>filename;
hugefilename = "huge"<>filename;

If[k==2,
filename = "high_"<>filename;
bigfilename = "high_"<>bigfilename;
hugefilename = "high_"<>hugefilename
];

Export[filename,heattab,"Table"];

Export[bigfilename,heattabbig,"Table"];

Export[hugefilename,heattabhuge,"Table"],
{k,1,1}];

If[ipass>2,
If[ipass<10,
chroot=(((((9.-ipass)/9.)))^2),
chroot=0],

chroot=1];

tripleroot[[2]]=(10.^36);
Tiroot[[2]]=100.(chroot)+50.(1.- chroot);

Print["Ti high = ", Tiroot[[2]], "   Ti low = ", Tiroot[[1]]];
Print["triple high = ", tripleroot[[2]], "   triple low = ", tripleroot[[1]]];

If[ipass<3,
τtab[[ipass+1]]=τEiloc;
tripletab[[ipass+1]]=Min[tripleroot],
τtab[[iMax]]=τEiloc;
tripletab[[iMax]]=Min[tripleroot]
];
(* It may be worth checking that the search is converging *)

iMax=First[Ordering[tripletab,-1]];
convdist=Max[Abs[tripletab[[1]]-tripletab[[2]]],Abs[tripletab[[1]]-tripletab[[3]]],Abs[tripletab[[3]]-tripletab[[2]]]];
convdistτ=Max[Abs[τtab[[1]]-τtab[[2]]],Abs[τtab[[1]]-τtab[[3]]],Abs[τtab[[3]]-τtab[[2]]]];

If[ipass<3,
convdist=1.,
convdist=(convdist/(Max[tripletab]));
convdistτ=((convdistτ)/(Max[τtab]))
];

If[(*Abs[Tiroot[[2]]-Tiroot[[1]]]>dTsearch,*)
convdist>dtriplesearch || convdistτ>dτsearch,

xx=τtab;
yy=tripletab;

If[ipass>1,
parsol=Solve[{ap ((xx[[1]])^2)+bp xx[[1 ]]+cp==yy[[1]],ap ((xx[[2]])^2)+bp xx[[2]] +cp==yy[[2]],ap ((xx[[3]])^2)+bp xx[[3]] +cp==yy[[3]]},{ap,bp,cp}];
aloc=First[ap/.parsol];
bloc=First[bp/.parsol];
δτloc=Min[Abs[τtab[[1]]+(bloc/(2 aloc))],Abs[τtab[[2]]+(bloc/(2 aloc))],Abs[τtab[[3]]+(bloc/(2 aloc))]];
τEiloc=-(bloc/(2 aloc));
If[aloc<0.||τEiloc<1.5||δτloc>0.5,
Print[" XXXXXXXXXXXXXX Failed parabola iteration XXXXXXXXXXXXXX"];
Print["Discarded τEiloc = ", τEiloc];
Print["aloc = ", aloc, "  δτloc = ", δτloc];
(* just do linear extrapolation using the smallest two values *)
iii1=Ordering[tripletab,1][[1]];
iii2=Ordering[tripletab,2][[2]];
τEiloc=Max[2τtab[[iii1]]-τtab[[iii2]],1.85+RandomReal[{0.74,1.5}]Abs[τtab[[iii1]]-τtab[[iii2]]]]

];

];

Print["τEis = ",τtab];
Print["triples = ",tripletab];
Print["XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"];
Print["XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"];
Print["XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"];
Print["continue root search: roots are too far"];
Print["XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"];
Print["XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"];
Print["XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"];

(* Define a new list of values to be tested: this should save A LOT of time*)

ipass=ipass+1;

Print["ipass = ", ipass];
Print["XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"];

Tcenter = ((Tiroot[[2]]+Tiroot[[1]])/2);
droots=Tiroot[[2]]-Tiroot[[1]];
nlist=Max[10,Ceiling[(droots/4)]];
dTlist = (droots/nlist);

valuelist1=Table[0,{ilist,1,nlist+2}];

valuelist1[[1]]=Max[Tiroot[[1]]-3,6.];
valuelist1[[2]]=Min[Tiroot[[2]]+3,180.];

valuelist1[[1]]=3.;
valuelist1[[2]]=120.;

Do[

valueloc=Tcenter+(((-1))^ilist)Floor[(ilist/2)]dTlist;
valuelist1[[ilist+2]]=valueloc,
{ilist,1,nlist}
];


Goto[repeatlabel];
]
,

{ivalues,16,30}]
(*{ivalues,1,nvalueslin}]*)