100 Reputation

2 years, 145 days

Appreciation...

@acer Thank you very much for your assistance.

Explanations...

Please can you explain what these codes do?

(minv,maxv):=[min,max](op([1,3],%))[]:
conts:=[seq(minv+(u-1)*(maxv-minv)/(nconts+2-1),u=1..nconts+2)][2..nconts+1]:
colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(u-1)*(c2-c1)/(N-1),
u=1..N)])([ColorTools:-Color(color2)[]],
[ColorTools:-Color(color1)[]],
nops(conts))):

Thank you very much.

Setcolors...

Please how do I set colors in listcontplot? I wanna know the values that correspond to each color in the contour plot.

 > restart:   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.1104e-2,                        k__1:= 6*10^(-12),                        d:= 0.2e-3,                        theta0:= 0.1e-3,                        eta__1:= 0.240e-1, chi:= 1.219*10^(-6),                        alpha:= 1-alpha__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^2-alpha__3*xi/eta__1 - chi*H^2)),                        w := chi*H^2*eta__1*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),                        n:= 20,                        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; # Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes                  g:= q-(1-alpha)*tan(q)- (w*q + c)*tan(q):                  f:= subs(q = x+I*y, g):                  b1:= evalc(Re(f)) = 0:                  b2:= evalc(Im(f)) = 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 Odd asymptote                  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; # Compute the rest of the Roots                  OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);                  AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));                  if   AreThereComplexRoots                  then gg:= [ qcomplex1,                              qcomplex2,                              op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),                              seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)                           ];                  elif not AreThereComplexRoots                  then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),                              seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)                            ];                  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] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2):                  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 off-diagonal 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): ## 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*z-d)*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): ##              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;                  end proc:
 > # Run the calculation for supplied value of 'xi'
 > # Put the returned quantities in a simple list
 > ans:=[doCalc(-0.06, 5)]: ans1:=[doCalc(0.06, 5)]: evalf(ans[9]); evalf(ans[10]);
 (1)
 > ##
 > with(plots): d:= 0.2e-3: ## director Plot negative activity display( plot(ans[1], z=0..d, color = green),          plot([seq(eval(ans[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legendstyle = [font                 = ["ROMAN", 12],location = right], labeldirections = ["horizontal", "vertical"]),                     axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020=                         "200"]], axis[2]=[tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006=                     "6e-05", 0.00008= "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," ()                     "), typeset(theta(z, t)," (degrees)")], labelfont = ["TimesNewRoman", 14],                             axesfont = [14, 14]); ## director angle Plot postive activity display( plot(ans1[1], z=0..d, color = green, legend = theta[0]),          plot([seq(eval(ans1[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legend = [t[1] =                  .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]), axis[1]=[tickmarks=                [0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]], axis[2]=                        [tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006="6e-05", 0.00008=                       "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," ()"), typeset(theta(z, t)," (degrees)")], labelfont =["TimesNewRoman", 14], axesfont = [14, 14]); ## v Plot negative activity plt2:= plot([seq(eval(ans[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]): display( plt2, axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," ()"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);   ## v Plot postive activity plt21:= plot([seq(eval(ans1[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legend = [t[1] =                .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]): display( plt21,axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," ()"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);
 > # q plot
 > #
 > plot(ans[8], q=0..5*Pi,view=[DEFAULT,-10..10], labels =[q, r(q)], labelfont = ["TimesNewRoman", 14], labeldirections = ["horizontal", "vertical"] );
 > # p values
 > evalf(ans[5]):
 > # Director angle plot for at d/2 as a funtion of time
 > theta_f := subs(z=d/2, Re(ans[2])): plot(theta_f, t=0..10, view=[DEFAULT, DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed); theta_f := subs(z=d/2, ans1[2]): plot(theta_f, t=0..1, view=[DEFAULT,DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed);
 > #Convert and store xi and doCalc into an array # ans1:= CodeTools:-Usage(Array([seq(seq( [j, u, doCalc(j,u)], j=-2..0, 1), u =0..2,1)])):
 > #
 > with(ListTools):
 > testproc := proc(j, u)    global calcvals:=doCalc(j,u);    evalf(calcvals[4]), evalf(Im(calcvals[7])); end proc;
 (2)
 > with(plots): setcolors(["BlueViolet","Coral"]): whatjs:=[op([seq(1..4)]),op([seq(6..10)])]: CodeTools:- Usage(listcontplot([seq([seq(testproc(i,j)[1], i=-2..0,0.1)],j=1..3)],filledregions=true));
 > CodeTools:- Usage(listcontplot([seq([seq(testproc(i,j)[2], i=-2..0, 0.1)],j=1..3)],filledregions=true));
 memory used=20.42GiB, alloc change=32.00MiB, cpu time=2.50m, real time=4.84m, gc time=13.31s
 > CodeTools:- Usage(listcontplot3d([seq([seq(testproc(i,j)[1], i=-2..0, 0.1)],j=1..3)],filledregions=true));
 memory used=20.42GiB, alloc change=0 bytes, cpu time=3.92m, real time=6.55m, gc time=19.97s
 > CodeTools:- Usage(listcontplot3d([seq([seq(testproc(i,j)[2], i=-2..0, 0.1)],j=1..3)],filledregions=true));
 memory used=20.41GiB, alloc change=0 bytes, cpu time=2.51m, real time=4.99m, gc time=13.38s
 > CodeTools:- Usage(matrixplot([seq([seq(testproc(i,j)[2], i=-2..0, 0.1)],j=1..3)],filledregions=true));
 > #for i from -2 to 0 by 0.1 do #   for j from 1 to 10 do #       testproc(i,j); #       print(i, j); # end do; # end do;
 >

Thank you very much.

@tomleslie

I really understand what is happening. In case, I have to use contourplot instead of listconplot.  Is there any way I can change.

CodeTools:- Usage(listconplot([seq([seq(testproc(i,j)[2], i=-2..0,0.1)],j=1..3)],filledregions=true));

to something say

Usage(contourplot([seq([seq(testproc(i,j)[2], i=-2..0,0.1)],j=1..3)],filledregions=true));

Because, I want the x-axis to be -2..0 and y 1..3.

I tried it but it returned an error.

Okay, thank you.

@Carl Love

Okay, thank you.

removing repeated entry without using Ma...

@tomleslie

Hi Tomleslie,

## please why do I have repeated values for j and u? And how can solve this problem without using Makeunique function.
## I also want to do contourplot for u vs j for complex values of p_min.

In the attached file you will see what I'm talking about at the end of the file.

Case_3_IjuptilK_300522.mw

@tomleslie Okay, thank you very muc...

@tomleslie

Okay, thank you very much.

@tomleslie

Okay, thanks

Legend to the centre...

Hi Tom,

Please, how can I place the legend in the center of the curve?

Okay, thanks

SCALING THE X AND Y AXIS IN PLOT...

@tomleslie

Hi Tom,

Please how can I change the scale of a plot: For instance, on the theta_plot, I wanna change 0.00005 to 50 micrometer on the x-axis and convert the 0.0001 to 1*10^-4(standard form) on the y-axis.

Scaling.mw