(*This notebook shows the use of the Mathematica library MatrixMethod2D.m to characterize the unstable of an unstable Mestel disc, with following Zans's notations : q=6, N=6. This allows the reproduction of the figure C.1 from Fouvry&Pichon&Magorrian&Chavanis (2015).
If ever you were to use this library, please cite Fouvry&PIchon&Magorrian&Chavanis (2015), where the details of this implementation are presented.*)

(*Place the library file MatrixMethod2D.m in the same directory as this test notebook MestelMode.m.*)
(*Loading the package.*)
<<(NotebookDirectory[]<>"MatrixMethod2D.m");
(*adressfiles is a variable from MatrixMethod2D controls the location where the output files will be saved.*)
adressfiles = NotebookDirectory[];

(*Definition of the parameters from MatrixMethod2D.m.*)
(*Definition of the basis used.*)
RKA=5;
kKA=7;
UbasisMES = UbasisKA;
(*Definition of the list of the basis elements (ψ^((p))) considered.*)
lmode = 2;
nmode = 7;
(*Creation of the list of the basis elements, via the variables : listbasis2D and nbelementsbasis2D. Since we are looking for mϕ=2 mode, listbasis won't be used afterwards.*)
Generatelistbasis2D[];
(*Bound to compute the sum on m1 in the response matrix expression.*)
m1sumMatbound = 7;
(*Construction of the (rp,ra)-grid.
One must have rGridminMat-ΔrGridMat > 0, so that listGridlarge can effectively be constructed.*)
rGridminMat = 0.15;
rGridmaxMat = 4.8;
ΔrGridMat = 0.1;
(*Creation of the list (rp,ra) grid-list in the variables : listGrid, nbelementsGrid, listGridlarge, nbelementsGridlarge.
The "large" grid is used to compute the derivatives via finite differences.*)
GeneratelistGrid[];
(*Getting the number of evaluation of WMat which will be needed.
This is the most expensive part of the computation. One should account for 1s to compute one element. However, this is made in parallel, and can therefore be significantly fastened when run on multiple cores.*)
GetnbevalsWMat2D[lmode,nmode,m1sumMatbound];
(*Definition of the physical parameters of the disc.*)
V0=1;
G = 1;
Ri=1;
νt= 6.;
μt = 5.;
R0=11.5;
Rmax = 20;
(*Unstable value of the velocity dispersion.*)
q = 6.;
σr0 = Sqrt[(V0^(2))/(q+1)];
(*Active fraction of the disc. The disc is fully massive.*)
ξ = 1.;

(*Definition of the total Mestel potential.*)
ψMES[r_?NumericQ] := V0^(2) Log[r/Rmax];

(*Definition of the system DF in the (E,L)-space.
One should be VERY careful of having the correct normalization constant.*)
(*Tapering functions.*)
Tinner[Lf_?NumericQ] := (Lf^(νt))/((Ri V0)^(νt) + Lf^(νt))
Touter[Lf_?NumericQ] := (1 + ((Lf)/(R0 V0))^(μt))^(-1)
CnormDF = (V0^(2))/(2^(1+q/2) π^(3/2) G σr0^(q+2) Gamma[(q/2)+(1/2)] Rmax^(q+1));
DFELMES[Ef_?NumericQ,Lf_?NumericQ] := ξ CnormDF (Lf)^(q) Exp[- (Ef)/(σr0^(2))] Tinner[Lf] Touter[Lf];

(*Computation and dump of the WMat elements, i.e. the Fourier transform w.r.t. to the angles of the basis elements. This is the most expensive part of the computation.
This part of the code is commented, and should not be executed, if one uses the data file DumpWMat2D_....txt.*)
(*PrepareWMat2D[lmode,nmode,m1sumMatbound,ψMES,UbasisMES]//AbsoluteTiming
DumpWMat2D[lmode,nmode,m1sumMatbound,rGridminMat,rGridmaxMat,ΔrGridMat,ψMES,UbasisMES];*)

(*Computing : E, J2=L, J1, Ω1, Ω2, the Jacobian (E,L)-> (rp,ra) and ∂F/∂J.
If one uses the prepared data files : DumpE_...txt, DumpJ2_...txt, DumpJ1_...txt, DumpOmega1_...txt, DumpOmega2_...txt, DumpJac_...txt, DumpDFGrad_...txt, this cell does not have to be executed.*)
(*PrepareΩ[ψMES]//AbsoluteTiming
DumpΩ[rGridminMat,rGridmaxMat,ΔrGridMat,ψMES];
PrepareJac[ψMES]//AbsoluteTiming
DumpJac[rGridminMat,rGridmaxMat,ΔrGridMat,ψMES];
PrepareDFGrad[lmode,m1sumMatbound,ψMES,DFELMES]//AbsoluteTiming
DumpDFGrad[lmode,m1sumMatbound,rGridminMat,rGridmaxMat,ΔrGridMat,ψMES,DFELMES]*)

(*Loading all the precomputed quantities : WMat, Frequencies, Jacobian and DF Gradients.
Be sure to have enough RAM on your system of execution.*)
LoadWMat2D[lmode,nmode,m1sumMatbound,rGridminMat,rGridmaxMat,ΔrGridMat,ψMES,UbasisMES]//AbsoluteTiming
LoadΩ[rGridminMat,rGridmaxMat,ΔrGridMat,ψMES]//AbsoluteTiming
LoadJac[rGridminMat,rGridmaxMat,ΔrGridMat,ψMES]//AbsoluteTiming
LoadDFGrad[lmode,m1sumMatbound,rGridminMat,rGridmaxMat,ΔrGridMat,ψMES,DFELMES]//AbsoluteTiming

(*One can then compute one element of the response matrix. The computation of the first matrix element is much longer than the subsequent one, because of the preparation of memoized functions.*)
(*For comparison, check that you obtain the same result (except of course for the timing of the execution).*)
ZMat2D[lmode,nmode,6.5,0.03,ψMES,UbasisMES,DFELMES]// AbsoluteTiming
ZMat2D[lmode,nmode,6.6,0.04,ψMES,UbasisMES,DFELMES]// AbsoluteTiming

(*Computation of numerous values of the response matrix for various (ω,η).
The final instruction of this block does not have to be executed, since the values are already given in the data file DumpZMatCoeff2D_...txt*)
(*Array of ω for which the computation of ZMat was performed.*)
tabω = {0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1.,1.01,1.02,1.03,1.04,1.05,1.06,1.07,1.08,1.09,1.1,1.11,1.12,1.13,1.14,1.15,1.16,1.17,1.18,1.19,1.2,1.21,1.22,1.23,1.24,1.25,1.26,1.27,1.28,1.29,1.3};
nbω = tabω//Length;
(*Array of η for which the computation of ZMat was performed.*)
tabη = {0.16,0.18,0.20,0.22,0.24};
nbη = tabη//Length;
(*Function to compute Det[I-qM(ω,η)].*)
TabDetq[η_,q_] :=Map[{Re[#],Im[#]}&,Table[Det[IdentityMatrix[nmode+1] - q ZMat2D[lmode,nmode,tabω[[ind]],η,ψMES,UbasisMES,DFELMES]],{ind,1,nbω,1}]];

(*Computation of ZMat for the all the (ω,η) elements in the arrays tabω and tabη. This does not need to be executed, since one can use the data file DataZMatCoeff2D_...txt.*)
(*Do[Print[TabDetq[tabη[[ind]],1.];//AbsoluteTiming],{ind,1,nbη,1}]*)

(*To dump the already computed values of ZMat, for immediate later use.*)
(*DumpZMatCoeff2D[lmode,nmode,m1sumMatbound,rGridminMat,rGridmaxMat,ΔrGridMat,ψMES,UbasisMES,DFELMES]*)

(*Using the given data file DumpZMatCoeff2D_...txt, one can load the already computed values of ZMatCoeff2D.*)
LoadZMatCoeff2D[lmode,nmode,m1sumMatbound,rGridminMat,rGridmaxMat,ΔrGridMat,ψMES,UbasisMES,DFELMES];//AbsoluteTiming

(*Reproducing the figures presented in the appendix C from Fouvry&Pichon&Magorrian&Chavanis. Most of the instructions are only cosmetic.*)
(*Values of η from tabη which will be represented.*)
ηp1 = 0.24;
ηp2 = 0.22;
ηp3 = 0.20;
ηp4 = 0.18;
ηp5 = 0.16;
(*Edges of the first figure.*)
xmincontour = -0.02;
xmaxcontour = 0.02;
ymincontour = -0.02;
ymaxcontour = 0.02;
(*Edges of the second figure.*)
xminnorm = 0.5;
xmaxnorm = 1.3;
yminnorm = -4;
ymaxnorm = -1;
(*Choice of the colors.*)
color1= ColorData["SolarColors"][0.1];
color2=ColorData["SolarColors"][0.3];
color3=ColorData["SolarColors"][0.5];
color4=ColorData["SolarColors"][0.7];
color5=ColorData["SolarColors"][0.9];
basestyleplot = {FontFamily -> "Avenir Next Condensed",FontSize -> 13};
(*Creating the axes labels for the plots.*)
s1 = Style["Re[det[","Avenir Next Condensed",13];
s2 = Style["I","Avenir Next Condensed",13,Bold];
s3 = Style["-","Avenir Next Condensed",13];
s4 = OverHat[Style["M","Avenir Next Condensed",13,Bold]];
s5 = Style["(ω)]]","Avenir Next Condensed",13];
sRefull = Apply[StringJoin,ToString[#,StandardForm]&/@{s1,s2,s3,s4,s5}];
s1 = Style["Im[det[","Avenir Next Condensed",13];
s2 = Style["I","Avenir Next Condensed",13,Bold];
s3 = Style["-","Avenir Next Condensed",13];
s4 = OverHat[Style["M","Avenir Next Condensed",13,Bold]];
s5 = Style["(ω)]]","Avenir Next Condensed",13];
sImfull = Apply[StringJoin,ToString[#,StandardForm]&/@{s1,s2,s3,s4,s5}];
s1 = Style["Log |det[","Avenir Next Condensed",13];
s2 = Style["I","Avenir Next Condensed",13,Bold];
s3 = Style["-","Avenir Next Condensed",13];
s4 = OverHat[Style["M","Avenir Next Condensed",13,Bold]];
s5 = Style["(ω)]|","Avenir Next Condensed",13];
sLogfull = Apply[StringJoin,ToString[#,StandardForm]&/@{s1,s2,s3,s4,s5}];
(*Plots of the Nyquist contours.*)
plcontour =  ListLinePlot[TabDetq[#,1.]&/@{ηp1,ηp2,ηp3,ηp4,ηp5},PlotRange->{{xmincontour,xmaxcontour},{ymincontour,ymaxcontour}},PlotStyle->{{Dashed,Thickness[0.007],color1},{Dashed,Thickness[0.007],color2},{Dashed,Thickness[0.007],color3},{Dashed,Thickness[0.007],color4},{Dashed,Thickness[0.007],color5}} ,BaseStyle->basestyleplot,AxesLabel->{sRefull,sImfull},PlotLegends->{"0.24","0.22","0.20","0.18","0.16"},LabelStyle->{FontFamily->"Avenir Next Condensed",FontSize->13},PlotTheme->"Classic"]
(*Export of the image of the first figure.*)
Export[adressfiles<>"contours_Nyquist_N_6.eps",Show[plcontour,ImageSize->300]];
(*Adding the line for the unstable pattern mode.*)
plline = Graphics[{Dotted,Thickness[0.003],Line[{{0.94,-4},{0.94,-0}}]}];
(*Adding an epilog.*)
epilognorm ={Text[Style["(Subscript[ω, 0])=0.94","Avenir Next Condensed",13],{0.955,-1.15},{-1,0}]};
(*Plots of the norms of the contours.*)
plnorm = Show[ListLinePlot[Transpose[{tabω,Log10@MapThread[Function[{x,y},Sqrt[x^2 + y^2]],TabDetq[#,1]//Transpose]}]&/@{ηp1,ηp2,ηp3,ηp4,ηp5},ImageSize->Large,PlotRange->{{xminnorm,xmaxnorm},{yminnorm,ymaxnorm}},PlotStyle->{{Dashed,Thickness[0.007],color1},{Dashed,Thickness[0.007],color2},{Dashed,Thickness[0.007],color3},{Dashed,Thickness[0.007],color4},{Dashed,Thickness[0.007],color5}},BaseStyle->basestyleplot,AxesLabel->{"(Subscript[ω, 0])",sLogfull},PlotLegends->{"0.24","0.22","0.20","0.18","0.16"},LabelStyle->{FontFamily->"Avenir Next Condensed",FontSize->13},PlotTheme->"Classic",AxesOrigin->{xminnorm,yminnorm},Epilog->epilognorm],plline]
(*Export of the image of the second figure.*)
Export[adressfiles<>"norms_Nyquist_N_6.eps",Show[plnorm,ImageSize->300]];