> 
doCalc:= proc( xi , u) #u is the \bar(H): normalize magnetic field magnitude,
# where H = bar(H)*H__a
# Import Packages
uses ArrayTools, Student:Calculus1, LinearAlgebra,
ListTools, RootFinding, plots, ListTools:
local gamma__1:= .1093,
alpha__3:= 0.1104e2,
k__1:= 6*10^(12),
d:= 0.2e3,
theta0:= 0.1e3,
eta__1:= 0.240e1, chi:= 1.219*10^(6),
alpha:= 1alpha__3^2/(gamma__1*eta__1),
theta_init:= theta0*sin(Pi*z/d),
H__a:= Pi*sqrt(k__1/chi)/d,
H := u*H__a,
c:=alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2alpha__3*xi/eta__1  chi*H^2)),
w := chi*H^2*alpha/(4*k__1*q^2/d^2alpha__3*xi/eta__1  chi*H^2),
n:= 50,
g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp, nstar,soln3, Imagroot1, Imagroot2, AreTherePurelyImaginaryRoots;
# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes
g:= q(1alpha)*tan(q)+ c*tan(q) + w*q:
f:= subs(q = x+I*y, g):
b1:= evalc(Re(f)) = 0:
b2:= y(1alpha)*tanh(y) (alpha__3*xi*alpha/(eta__1*(4*k__1*y^2/d^2+alpha__3*xi/eta__1)))*tanh(y) = 0;
qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);
# Compute Odd asymptote
ModifiedOddAsym:= abs(``~(OddAsymptotes, qstar));
qstarTemporary:= min(ModifiedOddAsym);
indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
qstar2:= OddAsymptotes(indexOfqstar2);
# Compute complex roots
AreThereComplexRoots:= type(true, 'truefalse');
try
soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = infinity .. 0});
qcomplex1:= subs(soln1, x+I*y);
qcomplex2:= subs(soln2, x+I*y);
catch:
AreThereComplexRoots:= type(FAIL, 'truefalse');
end try;
AreTherePurelyImaginaryRoots:= type(true, 'truefalse');
try
# Compute the rest of the Roots
soln3:= simplify~(fsolve({ b2}, { y = 0 .. infinity}),zero);
Imagroot1:=subs(soln3, I*y);
Imagroot2:= Imagroot1;
catch:
AreTherePurelyImaginaryRoots:= type(FAIL, 'truefalse');
end try;
## odd and all asymptotes
OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
if (xi > 0) then
if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
then gg:= [qcomplex1, qcomplex2,op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])),
seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif AreThereComplexRoots and AreTherePurelyImaginaryRoots
then gg:= [qcomplex1, qcomplex2, Imagroot2, op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])), seq(op (Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif not AreThereComplexRoots and AreTherePurelyImaginaryRoots
then gg:= [Imagroot2, op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif not AreThereComplexRoots and not AreTherePurelyImaginaryRoots
then gg:= [op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
end if:
else
if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
then gg:= [qcomplex1, qcomplex2,op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])),
seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif AreThereComplexRoots and AreTherePurelyImaginaryRoots
then gg:= [qcomplex1, qcomplex2, Imagroot2, op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])), seq(op (Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif not AreThereComplexRoots and AreTherePurelyImaginaryRoots
then gg:= [op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
elif not AreThereComplexRoots and not AreTherePurelyImaginaryRoots
then gg:= [op(Roots(g, q = 0.1e3 .. AllAsymptotes[1])), seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)];
end if:
end if:
# Remove the repeated roots if any & Redefine n
qq:= MakeUnique(gg):
m:= numelems(qq):
## Compute the time constants
for i to m do
p[i] := simplify(evalf(gamma__1*alpha/(4*k__1*qq[i]^2/d^2  alpha__3*xi/eta__1 chi*H^2)),':zero'):
end do:
## Compute θ_n from initial conditions
for i to m do
Efun[i] := cos(qq[i]2*qq[i]*z/d)cos(qq[i]):
end do:
## Compute integral coefficients for offdiagonal elements θ_n matrix
printlevel := 2:
for i to m do
for j from i+1 to m do
b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):
b[j, i] := b[i, j]:
aa[i, j] := b[i, j]:
aa[j, i] := b[j, i]:
end do :
end do:
## Calculate integral coefficients for diagonal elements in theta_n matrix
for i to m do
aa[i, i] := int(Efun[i]^2, z = 0 .. d):
end do:
## Calculate integrals for RHS vector
for i to m do
F[i] := int(theta_init*Efun[i], z = 0 .. d):
end do:
## Create matrix A and RHS vector B
A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):
B := Vector([seq(F[i], i = 1 .. m)]):
## Calculate inverse of A and solve A*x=B
#Ainv := 1/A:
#r := MatrixVectorMultiply(Ainv, B):
r := LinearSolve(A,B);
## Define Theta(z,t)
theta_sol := add(r[i]*Efun[i]*exp(t/p[i]), i = 1 .. m):
## Compute v_n for n times constant
for i to m do
v[i] := (2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]) + alpha__3*chi*H^2/(2*eta__1*qq[i]*gamma__1*alpha):
end do:
## Compute v(z,t) from initial conditions
for i to m do
Vfun[i] := d*sin(qq[i]2*qq[i]*z/d)+(2*zd)*sin(qq[i]):
end do:
## Define v(z,t)
v_sol := add(v[i]*r[i]*Vfun[i]*exp(t/p[i]), i = 1 .. m):
## compute minimum value of tau
minp:=min(seq( Re(p[i]), i = 1 .. m) );
member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):
## Return all the plots
return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, H, H__a,qq;
end proc:


Run the calculation for supplied value of 'xi



> 
ans:=[doCalc(0.1, 2)]:
evalf(ans[9]);
evalf(ans[10]);


Contour plot for pmin (minimum of tau) for different xi and H values


> 
with(plots):
d:= 0.2e3:

> 
testproc := proc(j, u, k) option remember;
local calcvals;
if not [j,u,k]::[numeric,numeric,posint] then
return 'procname'(args);
end if;
calcvals:=doCalc(j,u);
evalf(calcvals[k]);
end proc:

> 
(gridji,rngj,rngi):=[5,5],0..3,2.0..0.0;
PP[gridji,4]:=plot3d(max(1e2,testproc(i,j,4)), j=rngj, i=rngi, grid=gridji):
PP[gridji,7]:= plot3d(Im(testproc(i,j,7)), j=rngj, i=rngi, grid=gridji):

> 
## Surface plot using 3d
PP[gridji,4]:

> 
## This is the pmin values
data[gridji,4]:=op([1,3],PP[gridji,4]);

> 
## Convert Array to Matrix to allow the use listcontplot
ConMatrix := Matrix(data[gridji,4]);

> 
## Assign the min and max value of pmin to minv and maxv
(minv,maxv):=[min[defined],max[defined]](data[gridji,4])[];

> 
## GRB color and the contours
(color1,color2):="Red","Yellow":
conts:=[minv, seq(minv+(i1)*(maxvminv)/(NN+21),i=2..NN+1), maxv]/2.0;

> 
colorlist:=ColorTools:Color~(((c1,c2,N)>[seq(c1+(i1)*(c2c1)/(N1),
i=1..N)])([ColorTools:Color(color2)[]],
[ColorTools:Color(color1)[]],
nops(conts))):

> 
## Use the Conmatrix and the conts to plot contours and transform the axes
subsindets(plots:listcontplot(ConMatrix, contours=conts,filledregions,thickness=0),
specfunc(anything,THICKNESS),u>THICKNESS(0.2)):
LCP[gridji,4]:=plottools:transform((x,y)>[lhs(rngj)+(rhs(rngj)lhs(rngj))*(x1)/(gridji[1]1),
lhs(rngi)+(rhs(rngi)lhs(rngi))*(y1)/(gridji[2]1)])(%):
display(LCP[gridji,4],
seq(plot(2.0,1..1,legend=evalf[4](conts[i]),thickness=15,
color=colorlist[nops(conts)i+1]),
i=1..nops(conts)),
size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

> 
## This deals with the Imaginary part

> 
## Surface plot using plot3d for Imginary part
PP[gridji,7]:

> 
data[gridji,7]:=op([1,3],PP[gridji,7]):

> 
## Convert Array to matrix
AMatrix := Matrix(data[gridji,7]);

> 
(minv,maxv):=[min,max](data[gridji,7])[];

> 
(color1,color2):="Red","Yellow":
conts1:=[minv, seq(minv+(i1)*(maxvminv)/(NN+21),i=2..NN+1), maxv]/~1.5;

> 
colorlist:=ColorTools:Color~(((c1,c2,N)>[seq(c1+(i1)*(c2c1)/(N1),
i=1..N)])([ColorTools:Color(color2)[]],
[ColorTools:Color(color1)[]],
nops(conts1))):

> 
subsindets(plots:listcontplot(AMatrix,
contours=conts1~1e9,filledregions,thickness=0),
specfunc(anything,THICKNESS),u>THICKNESS(0.2)):
LCP[gridji,7]:=plottools:transform((x,y)>[lhs(rngj)+(rhs(rngj)lhs(rngj))*(x1)/(gridji[1]1),
lhs(rngi)+(rhs(rngi)lhs(rngi))*(y1)/(gridji[2]1)])(%):
display(LCP[gridji,7],
seq(plot(2.0,1..1,legend=evalf[4](conts1[i]),thickness=15,
color=colorlist[nops(conts)i+1]),
i=1..nops(conts1)),
size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);


