INTRODUCTION AND DESCRIPTION

(*This notebook was written so as to perform the angle-action calculations presented in Fouvry&Pichon&Magorrian&Chavanis (2015). If you were to use this library, please cite Fouvry&Pichon&Magorrian&Chavanis (2015), where the details of this implementation are presented.
The aim of this notebook is to create a general framework in which linear analysis may be perfomed.
This notebook contains :
+ Determination of the actions and intrinsic frequencies thanks to the pericenter and apocenter.
+ 2D computation of the Fourier transform of the basis elements.
+ 2D calculation of the response matrix
+ Discretized sub-region integration for the response matrix
+ Preparation of the calculation of the angle-actions transforms in 2D.
*)

(*In order to share the package context.*)
$DistributedContexts = Automatic;
(*We already launch the parallel kernels.*)
LaunchKernels[];

BeginPackage["MatrixMethod2D`"];

USAGES OF THE FUNCTIONS

(*USAGE SHARED VARIABLES*)
G::usage = "G. Gravitational constant.";
adressfiles::usage = "adressfiles. Name of the directory where the files will be saved.";
m1sumMatbound::usage = "m1sumMatbound. Bound for the sum on m1 used to determine the response matrix.";
ηMat ::usage = "ηMat. Small imaginary frequency used to regularize the expression of the response matrix.";
RKA::usage = "RKA. Scale radius for the Kalnajs (1976) 2D basis.";
kKA::usage = "kKA. Index kKA for the Kalnajs (1976) 2D basis.";
rGridminMat ::usage = "rGridminMat. Minimum radius for the (rperi,rapo)-grid used for the response matrix integration.";
rGridmaxMat ::usage = "rGridmaxMat. Maximum radius for the (rperi,rapo)-grid used for the response matrix integration.";
ΔrGridMat ::usage = "ΔrGridMat. Step distance radius for the (rperi,rapo)-grid used for the response matrix integration.";
listGrid ::usage = "listGrid. List containing the (rperi-rapo)-grid used for the response matrix integration.";
nbelementsGrid ::usage = "nbelementsGrid. Number of elements of the (rperi,rapo)-list listGrid used for the response matrix integration.";
listGridlarge::usage = "listGridlarge. List containing the large (rperi-rapo)-grid used in PrepareWMat before performing the response matrix integration.";
nbelementsGridlarge::usage = "nbelementsGridlarge. Number of elements of the large (rperi,rapo)-grid listGridlarge used un PrepareWMat before performing the response matrix integration.";
lbasismax::usage = "lbasismax. Number of l modes considered in the list of basis elements listbasis2D.";
nbasismax ::usage = "nbasismax. Number of n elements considered in the list of basis elements listbasis 2D.";
listbasis2D ::usage = "listbasis2D. List (lp,np) of the basis elements considered.";
nbelementsbasis2D ::usage = "nbelementsbasis2D. Number of elements of the list of basis elements listbasis2D.";
(*END USAGE.*)


(*USAGE GRIDS AND LISTS GENERATIONS.*)
choprad::usage = "choprad[r]. Used to sanitize the grid so that rperi and rapo may be used as arguments for memoized functions, i.e. have 1.2 instead of 1.2000000000000002.";
GeneratelistGrid::usage = "GeneratelistGrid[]. Updates the (rperi,rapo) grid via the variables listGrid, listGridlarge, nbelementsGrid and nbelementsGridlarge.";
CreatelistGridlarge::usage = 
"CreatelistGridlarge[]. Once listGrid, nbelementsGrid and ΔrGridMat have been specified, creates listGridlarge and nbelementsGridlarge."
rperifromGrid::usage = "rperifromGrid[i]. rperi of the i-th element of the listGrid (rperi,rapo)-grid.";
rapofromGrid::usage = "rapofromGrid[i]. rapo of the i-th element of the listGrid (rperi,rapo)-grid.";
rperifromGridlarge::usage = "rperifromGridlarge[i]. rperi of the i-th element of the listGridlarge (rperi,rapo)-grid.";
rapofromGridlarge::usage = "rapofromGridlarge[i]. rapo of the i-th element of the listGridlarge (rperi,rapo)-grid.";
Generatelistbasis2D::usage = "Generatelistbasis2D[]. Updates the list of the basis elements via the variables listbasis2D and nbelementsbasis2D.";
lfromlistbasis2D::usage = "lfromlistbasis2D[p]. l-index of the p-th basis element from listbasis2D.";
nfromlistbasis2D::usage = "nfromlistbasis2D[p]. n-index of the p-th basis element from listbasis2D.";
(*END USAGE.*)


(*USAGE EXPLICIT ACTION MAPPING.*)
Efromrad::usage="Efromrad[rperi,rapo,ψ]. Energy from pericenter and apocenter.";
J2fromrad::usage="J2fromrad[rperi,rapo,ψ]. Jϕ from pericenter and apocenter.";
J1fromrad::usage = "J1fromrad[rperi,rapo,ψ]. Jr from pericenter and apocenter.";
Ω1fromrad::usage = "Ω1fromrad[rperi,rapo,ψ]. Ω1 = κ from pericenter and apocenter.";
Ω2fromrad::usage = "Ω2fromrad[rperi,rapo,ψ]. Ω2 = Ωϕ from pericenter and apocenter";
DΩ1Drp::usage = "DΩ1Drp[rperi,rapo,ψf]. Finite differences estimation on the grid of ∂Ω1/∂rp.";
DΩ1Dra::usage = "DΩ1Dra[rperi,rapo,ψf]. Finite differences estimation on the grid of ∂Ω1/∂ra.";
DΩ2Drp::usage = "DΩ2Drp[rperi,rapo,ψf]. Finite differences estimation on the grid of ∂Ω2/∂rp.";
DΩ2Dra::usage = "DΩ2Dra[rperi,rapo,ψf]. Finite differences estimation on the grid of ∂Ω2/∂ra.";
(*END USAGE.*)


(*USAGE 2D BASIS KALNAJS (1976).*)
ψbasisKA::usage = "ψbasisKA[l,n,r,ϕ]. 2D potential basis elements from Kalnajs (1976).";
UbasisKA::usage = "UbasisKA[l,n,r]. Radial component of the potential basis elements from Kalnajs (1976).";
ρbasisKA::usage = "ρbasisKA[l,n,r,ϕ]. 2D surface density basis elements from Kalnajs (1976).";
DbasisKA::usage = "DbasisKA[l,n,r]. Radial component of the surface density basis elements from Kalnajs (1976).";
PCoeffKA::usage = "PCoeffKA[k,l,n]. Pre-coefficient for the potential basis elements from Kalnajs (1976).";
αCoeffKA::usage = "αCoeffKA[k,l,n,i,j]. Coefficient in the sum for the potential basis elements from Kalnajs (1976).";
SCoeffKA::usage = "SCoeffKA[k,l,n]. Pre-coefficient for the surface density basis elements from Kalnajs (1976).";
βCoeffKA::usage = "βCoeffKA[k,l,n,j]. Coefficient in the sum for the surface density basis elements from Kalnajs (1976).";
(*.END USAGE.*)


(*USAGE FOURIER TRANSFORM 2D.*)
ψbasisFT2D::usage = "ψbasisFT2D[l,n,m1,m2,rperi,rapo,ψ,Ubasis]. Angles Fourier transform of 2D basis elements.";
WMat2D::usage = "WMat2D[m1,l,m2,n,rperi,rapo,ψ,Ubasis]. Matrix WMat2D used to compute the 2D angles Fourier transform. (Expensive to compute).";
rtodθ1dr::usage = "rtodθ1dr[E,J2,Ω1,r,ψ]. Jacobian dθ1/dr(r) used in the estimation of WMat.";
rtoθ1::usage = "rtoθ1[rperi,E,J2,Ω1,r,ψ]. Mapping r -> θ1(r) used in the estimation of WMat.";
rtoθ2minusψ::usage = "rtoθ2minusψ[rperi,E,J2,Ω2,r,ψ]. Mapping r -> (θ2 - ψ)(r) used in the estimation of WMat.";
(*END USAGE.*)


(*USAGE RESPONSE MATRIX 2D.*)
RepMatCoeff2D::usage = "RepMatCoeff2D[lp,np,lq,nq,ω,η,ψ,Ubasis,DFEL]. One coefficient of the 2D response matrix M.";
ZMatCoeff2D::usage = "ZMatCoeff2D[lp,np,nq,ω,η,ψ,Ubasis,DFEL]. One coefficient of the bloc-diagonal 2D response matrix Z.";
agMat2D::usage = "agMat2D[lp,np,nq,m1,m2,rperi,rapo,ψ,Ubasis,DFEL]. Coefficient in the limited development of the numerator gMat2D to compute the 2D response matrix. (Constant contribution).";
bgMat2D::usage = "bgMat2D[lp,np,nq,m1,m2,rperi,rapo,ψ,Ubasis,DFEL]. Coefficient in the limited development of the numerator gMat2D to compute the 2D response matrix. (Contribution prop. to rperi).";
cgMat2D::usage = "cgMat2D[lp,np,nq,m1,m2,rperi,rapo,ψ,Ubasis,DFEL]. Coefficient in the limited development of the numerator gMat2D to compute the 2D response matrix. (Contribution prop. to rapo).";
ahMat2D::usage = "ahMat2D[m1,m2,ω,rperi,rapo,ψ]. Coefficient in the limited development of the denominator hMat2D to compute the 2D response matrix. (Constant contribution).";
bhMat2D::usage = "bhMat2D[m1,m2,rperi,rapo,ψ]. Coefficient in the limited development of the denominator hMat2D to compute the 2D response matrix. (Contribution prop. to rperi).";
chMat2D::usage = "chMat2D[m1,m2,rperi,rapo,ψ]. Coefficient in the limited development of the denominator hMat2D to compute the 2D response matrix. (Contribution prop. to rapo).";
gMat2D::usage = "gMat2D[lp,np,nq,m1,m2,rperi,rapo,ψ,Ubasis,DFEL]. Numerator of the integrated function to obtain the 2D response matrix.";
hMat2D::usage = "hMat2D[m1,m2,ω,rperi,rapo,ψ]. Denominator of the integrated function to obtain the 2D response matrix.";
JacELtorpra::usage = "JacELtorpra[rperi,rapo,ψ]. Jacobian d(E,L)/d(rp,ra) of the transformation (E,L) -> (rp,ra) present in gMat.";
DEDrp::usage = "DEDrp[rperi,rapo,ψ]. Partial derivative dE/drperi used in JacELtorpra.";
DEDra::usage = "DEDra[rperi,rapo,ψ]. Partial derivative dE/drapo used in JacELtorpra.";
DLDrp::usage = "DLDrp[rperi,rapo,ψ]. Partial derivative dL/drperi used in JacELtorpra.";
DLDra::usage = "DLDra[rperi,rapo,ψ]. Partial derivative dL/drapo used in JacELtorpra.";
DFGrad::usage = "DFGrad[m1,m2,rperi,rapo,ψ,DFEL]. Gradient of the DF (m1,m2).(dF/dJr,dF/dL) appearing in gMat.";
DFDE::usage = "DFDE[E,L,DFEL]. Numerical derivative dF/dE of the distribution function F(E,L).";
DFDL::usage = "DFDL[E,L,DFEL]. Numerical derivative dF/dL of the distribution function F(E,L).";
RepMat2D::usage = "RepMat2D[ω,η,ψ,Ubasis,DFEL]. 2D response matrix M(ω,η).";
ZMat2D::usage = "ZMat2D[lmode,nmode,ω,η,ψ,Ubasis,DFEL]. 2D bloc diagonal matrix ZMat2D.";
DetZMat2D::usage = "DetZMat2D[lmode,nmode,ω,η,ψ,Ubasis,DFEL]. Det(I - Z(ω,η)) used to seek for modes for 2D systems.";
(*END USAGE.*)


(*USAGE SUB-REGION INTEGRATION MATRIX.*)
ℵMat::usage = "ℵMat[ag,bg,cg,ah,bh,ch,η,Δr]. Analytical function to compute int_{-Δr/2}^{Δr/2} dxp int_{-Δr/2}^{Δr/2} dxa (ag + bg xp + cg xa)/(ah + bh xp + ch xa + η), used in the sub-region matrix integration.";
ℵMatD::usage = "ℵMatD[b,c,e,f,η]. Analytical function to compute int_{1/2}^{1/2} dx int_{-1/2}^{1/2} dy (1 + b x + c y)/(1 + e x + f y) used in the sub-region matrix integration.";
ℵMatAna::usage = "ℵMatAna[ag,bg,cg,ah,bh,ch,η,Δr]. Analytical function to compute int_{-Δr/2}^{Δr/2} dxp int_{-Δr/2}^{Δr/2} dxa (ag + bg xp + cg xa)/(ah + bh xp + ch xa + η), used in the sub-region matrix integration, where it always using ℵMatDAna.";
ℵMatDAna::usage = "ℵMatDAna[b,c,e,f,η]. Analytical integration of ℵMatD. ";
(*END USAGE.*)


(*USAGE PREPARATION CALCULATIONS 2D.*)
GetnbevalsWMat2D::usage = "GetnbevalsWMat2D[lmode,nmode,m1bound]. To determine the number of evaluations of WMat2D required, and the expected time of computation.";
PrepareWMat2D::usage = "PrepareWMat2D[lmode,nmode,m1bound,ψ,Ubasis]. Pre-compute in parallel the needed values of WMat2D.";
DumpWMat2D::usage = "DumpWMat2D[lmode,nmode,m1bound,rGridmin,rGridmax,ΔrGrid,ψ,Ubasis]. To dump the downvalues of WMat2D previously calculated with PrepareWMat2D.";
LoadWMat2D::usage = "LoadWMat2D[lmode,nmode,m1bound,rGridmin,rGridmax,ΔrGrid,ψ,Ubasis]. To load already computed values of WMat2D.";
PrepareΩ::usage = "PrepareΩ[ψ]. Pre-compute in parallel the needed values of E,J2,J1,Ω1 and Ω2.";
DumpΩ::usage = "DumpΩ[rGridmin,rGridmax,ΔrGrid,ψ]. Dump the downvalues of Efromrad, J2fromrad,J1fromrad,Ω1fromrad,Ω1fromrad previously calculated with PrepareΩ.";
LoadΩ::usage = "LoadΩ[rGridmin,rGridmax,ΔrGrid,ψ]. Load already computed values of Efromrad, J2fromrad, J1fromrad,Ω1fromrad,Ω1fromrad.";
PrepareJac::usage = "PrepareJac[ψ]. Pre-compute in parallel the needed values of JacELtorpra.";
DumpJac::usage = "DumpJac[rGridmin,rGridmax,ΔrGrid,ψ]. Dump the downvalues of JacELtorpra previously calculated with PrepareΩ.";
LoadJac::usage = "LoadJac[rGridmin,rGridmax,ΔrGrid,ψ]. Load already computed values of JacELtorpra.";
PrepareDFGrad::usage = "PrepareDFGrad[lmode,m1bound,ψ,DFEL]. Pre-compute in parallel the needed values of DFGrad for the 2D case.";
DumpDFGrad::usage = "DumpDFGrad[lmode,m1bound,rGridmin,rGridmax,ΔrGrid,ψ,DFEL]. Dump the downvalues of DFGrad previously calculated with PrepareDFGrad.";
LoadDFGrad::usage = "LoadDFGrad[lmode,m1bound,rGridmin,rGridmax,ΔrGrid,ψ,DFEL]. Load already computed values of DFGrad.";
DumpZMatCoeff2D::usage = "DumpZMatCoeff2D[lmode,nmode,m1bound,rGridmin,rGridmax,ΔrGrid,ψ,Ubasis,DFEL]. To dump the downvalues of ZMatCoeff2D already computed.";
LoadZMatCoeff2D::usage = "LoadZMatCoeff2D[lmode,nmode,m1bound,rGridmin,rGridmax,ΔrGrid,ψ,Ubasis,DFEL]. To load the already computed values of ZMatCoeff2D.";
(*END USAGE.*)

SHARED VARIABLES OF THE LIBRARY AND PRINTING

(*In this part of the notebook, we introduce all the variables which will be shared with the test cases and will therefore be modifiable.*)
Print["Linear Theory"];
Print[];
Print["Modifiable shared variables :"];

(*This package is required in order to be able to compute the numerical derivatives.
It also has to be loaded for all parallel kernels.*)
Needs["NumericalCalculus`"]
ParallelNeeds["NumericalCalculus`"]

(*Adress where the data files will be dumped.*)
adressfiles = NotebookDirectory[];
(*adressfiles = "~/Project/Mathematica_Runs/data/Cusp_Core/";*)
Print["adressfiles=",adressfiles];

G=1;
Print["G=",G];


(*Maximum sum index for m1 to compute the response matrix.*)
m1sumMatbound = 2;
Print["Maximum sum index for m1 to compute the response matrix"];
Print["m1sumMatbound=",m1sumMatbound];

(*We introduce a complex number ηMat to have well-defined integral.*)
ηMat = 0.01;
Print["Small imaginary part used to compute a well-defined response matrix integral"];
Print["ηMat=",ηMat];



Print["Potentials already defined :"];
ψKep[r_?NumericQ] := -1/r;
ψharm[r_?NumericQ] := 1/2 r^(2);
ψHernquist[r_?NumericQ]:=- 4 π G (1)/(2(1 + r));
ψJaffe[r_?NumericQ] :=-  G Log[1 + 1/r];
ψiso[r_?NumericQ] :=- (G)/(1 + Sqrt[1 + r^(2)]);
ψPlummer[r_?NumericQ] := - (G)/(Sqrt[r^(2) + 1]);
ψNFW[r_?NumericQ] := - (G) Log[1+r]/(r);
Print["ψKep, ","ψharm, ","ψHernquist, ","ψJaffe, ","ψiso, ","ψPlummer, ","ψNFW."];
Print[];


Print["Density/Potentials already defined :"];
(*Radius and k index for the Kalnajs (1976) 2D basis elements.*)
RKA = 1;
kKA = 4;
Print["2D Kalnajs (1976), RKA=", RKA," ; kKA=",kKA];
Print[];



Print["Interpolating Grid listGrid in (rperi,rapo) defined with the parameters :"];
(*Parameters for the grid in (rp,ra) to estimate the response matrix.*)
(*PARAMETERS TO GENERATE THE GRID IN RADIUS FOR THE MATRIX CALCULATION.*)
rGridminMat = 1.1;
Print["rGridminMat=",rGridminMat];
rGridmaxMat = 1.3;
Print["rGridmaxMat=",rGridmaxMat];
ΔrGridMat = 0.1;
Print["ΔrGridMat=",ΔrGridMat];
(*Initialization of the full grid of the pericenters and apocenters.
We note that the apocenters and pericenters have to be separated by at least 2 ΔrGridMat, so that the finite differences method used to compute a, b and c works.*)
listGrid = Flatten[Table[{choprad[rperitab],choprad[rapotab]},{rperitab,rGridminMat,rGridmaxMat,ΔrGridMat},{rapotab,rperitab+2 ΔrGridMat,rGridmaxMat,ΔrGridMat}],1];
nbelementsGrid = Length[listGrid];
Print["nbelementsGrid=",nbelementsGrid];

(*We define a "large" grid of pericenters and apocenters, which are all the positions on which gMat will effectively be estimated.
As a consequence, the minimal distance between rperi and rapo has been reduced to ΔrGridMat.
Indeed, we use the finite differences method to estimate the derivatives of gMat and hMat, so that the effective grid of evaluation of these functions is larger than the one used to grid the (rperi,rapo)-space.
We note here that the deletion of the duplicates has to be made while using Chop.*)
listGridlarge =
listGrid /. {rp_,ra_} -> {{choprad[rp - ΔrGridMat],choprad[ra]},{choprad[rp+ΔrGridMat],choprad[ra]},{choprad[rp],choprad[ra]},{choprad[rp],choprad[ra-ΔrGridMat]},{choprad[rp],choprad[ra+ΔrGridMat]}}//Flatten[#,1]&;
listGridlarge = DeleteDuplicates[listGridlarge,(Chop[#1[[1]]-#2[[1]]] == 0) && (Chop[#1[[2]]-#2[[2]]] == 0)&];
nbelementsGridlarge = Length[listGridlarge];
Print["nbelementsGridlarge=",nbelementsGridlarge];
Print["The grid can be regenerated with GeneratelistGrid[]"];
Print[];

(*We need to define the basis elements (p) listbasis2D which will be considered.
Definition of the boundaries for the basis elements.
We have the constraints : |l| ≤ lbasismax, 1 ≤ n ≤ nbasismax.
This variable has to be shared between all the kernels.*)
Print["List of the basis elements listbasis2D defined with the parameters :"];
lbasismax = 2;
Print["lbasismax=",lbasismax];
nbasismax = 2;
Print["nbasismax=", nbasismax];
listbasis2D = Flatten[Table[{ltab,ntab},{ltab,-lbasismax,lbasismax,1},{ntab,1,nbasismax,1}],1];
nbelementsbasis2D = Length[listbasis2D];
Print["nbelementsbasis2D=",nbelementsbasis2D];
Print["The 2D basis list can be regenerated with Generatelistbasis2D[]"];
Print[];
(*?MatrixMethod2D`* *)

Begin["Private`"];

GRIDS AND LIST GENERATIONS

(*Functions to generate the grid in pericenters and apocenters.*)

(*In order to be able to use rperi and rapo as good arguments for memoized functions (i.e have 1.2 instead of 1.2000000000000002), one always has to sanitize the radius by using choprad, which gets rid of the irrelevant decimals.
The parameter 10^(8) for the approximation is hard-coded.*)
choprad[r_?NumericQ]:= r 10^(8) // Round // #/10^(8)&//N;

(*It updates the variables listGrid, listGridlarge, nbelementsGrid and nbelementsGridlarge.*)
(*We note that the apocenters and pericenters have to be separated by at least 2 ΔrGridMat, so that the finite differences method used to compute a, b and c works.*)
GeneratelistGrid[] := Block[{},listGrid = Flatten[Table[{choprad[rperitab],choprad[rapotab]},{rperitab,rGridminMat,rGridmaxMat,ΔrGridMat},{rapotab,rperitab+2 ΔrGridMat,rGridmaxMat,ΔrGridMat}],1];
nbelementsGrid = Length[listGrid];
Print["nbelementsGrid=",nbelementsGrid];
listGridlarge =
listGrid /. {rp_,ra_} -> {{choprad[rp - ΔrGridMat],choprad[ra]},{choprad[rp+ΔrGridMat],choprad[ra]},{choprad[rp],choprad[ra]},{choprad[rp],choprad[ra-ΔrGridMat]},{choprad[rp],choprad[ra+ΔrGridMat]}}//Flatten[#,1]&;
listGridlarge = DeleteDuplicates[listGridlarge];
nbelementsGridlarge = Length[listGridlarge];
Print["nbelementsGridlarge=",nbelementsGridlarge];];

(*Function to thicken a list of (rp,ra).
One must have already defined listGrid, nbelementsGrid and ΔrGridMat.*)
CreatelistGridlarge[] := 
Block[{},
(*In case nbelementsGrid had not been specified.*)
nbelementsGrid = Length[listGrid];
listGridlarge = listGrid /. {rp_,ra_} -> {{choprad[rp - ΔrGridMat],choprad[ra]},{choprad[rp+ΔrGridMat],choprad[ra]},{choprad[rp],choprad[ra]},{choprad[rp],choprad[ra-ΔrGridMat]},{choprad[rp],choprad[ra+ΔrGridMat]}}//Flatten[#,1]&;
listGridlarge = DeleteDuplicates[listGridlarge];
nbelementsGridlarge = Length[listGridlarge];
Print["nbelementsGridlarge=",nbelementsGridlarge];];

(*Definition of the functions i -> rp(i) and i-> ra(i) to simplify the coding.*)
rperifromGrid[i_]:= listGrid[[i,1]];
rapofromGrid[i_]:= listGrid[[i,2]];

(*Definition of the functions i ->rp(i) and i-> ra(i) to simplify the coding for the large grid.*)
rperifromGridlarge[i_]:= listGridlarge[[i,1]];
rapofromGridlarge[i_]:= listGridlarge[[i,2]];


(*Function to generate the list of the 2D basis elements.
It updates the variables listbasis2D and nbelementsbasis2D.
We have the constraints : |l| ≤ lbasismax, 0 ≤ n ≤ nbasismax.
ATTENTION : n starts at 0.*)
Generatelistbasis2D[] := 
Block[{},listbasis2D = Flatten[Table[{ltab,ntab},{ltab,-lbasismax,lbasismax,1},{ntab,0,nbasismax,1}],1];
nbelementsbasis2D = Length[listbasis2D];];


(*Definition of the functions p->l(p),p->n(p) .*)
lfromlistbasis2D[p_]:=listbasis2D[[p,1]];
nfromlistbasis2D[p_]:=listbasis2D[[p,2]];

EXPLICIT ACTION MAPPING

(*In this section, we give the important functions allowing to obtain the important quantities describing an orbit : E, L, Jr, Ω1, Ω2
The main simplifying idea is to represents orbits with the quantities (rperi and rapo) instead of (E,L) or (J1,J2).*)
(*In order not to perform these integrations two times, the important functions will have to be memoized.*)

(*Energy as a function of the peri/apocenters*)
(*DEPENDS ON THE POTENTIAL ψ0*)
Efromrad[rperi_?NumericQ,rapo_?NumericQ,ψf_] :=Efromrad[rperi,rapo,ψf]=(rapo^(2) ψf[rapo] - rperi^(2) ψf[rperi])/(rapo^(2) - rperi^(2));
(*SetSharedFunction[Efromrad];*)

(*Angular momentum as a function of the peri/apocenters*)
(*DEPENDS ON THE POTENTIAL ψ0*)
J2fromrad[rperi_?NumericQ,rapo_?NumericQ,ψf_] := J2fromrad[rperi,rapo,ψf]=Sqrt[(2 (ψf[rapo] - ψf[rperi]))/(rperi^(-2) - rapo^(-2))];
(*SetSharedFunction[J2fromrad];*)

(*Radial action as a function of the peri/apocenters*)
(*DEPENDS ON THE POTENTIAL ψ0*)
J1fromrad[rperi_?NumericQ,rapo_?NumericQ,ψf_] := J1fromrad[rperi,rapo,ψf]=Block[{Eloc,J2loc},
Eloc=Efromrad[rperi,rapo,ψf];
J2loc=J2fromrad[rperi,rapo,ψf];
1/(π) NIntegrate[Sqrt[2(Eloc - ψf[rint]) - J2loc^(2)/rint^(2)],{rint,rperi,rapo},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}]// Re];
(*SetSharedFunction[J1fromrad];*)

(*Radial frequency as a function of the peri/apocenters*)
(*DEPENDS ON THE POTENTIAL ψ0.
We take the real part, in order not to have any complex components.*)
Ω1fromrad[rperi_?NumericQ,rapo_?NumericQ,ψf_] := Ω1fromrad[rperi,rapo,ψf]=Block[{Eloc,J2loc},
Eloc=Efromrad[rperi,rapo,ψf];
J2loc=J2fromrad[rperi,rapo,ψf];
π/(NIntegrate[1/(Sqrt[2(Eloc - ψf[rint]) - J2loc^(2)/rint^(2)]),{rint,rperi,rapo},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}])//Re];
(*SetSharedFunction[Ω1fromrad];*)

(*Azimuthal frequency as a function of the peri/apocenters*)
(*DEPENDS ON THE POTENTIAL ψ0,
We take the real part in order not to have any complex components.*)
Ω2fromrad[rperi_?NumericQ,rapo_?NumericQ,ψf_]:= Ω2fromrad[rperi,rapo,ψf]=Block[{Eloc,J2loc,Ω1loc},
Eloc = Efromrad[rperi,rapo,ψf];
J2loc = J2fromrad[rperi,rapo,ψf];
Ω1loc = Ω1fromrad[rperi,rapo,ψf];
Ω1loc J2loc /(π) NIntegrate[(1)/(rint^(2) Sqrt[2(Eloc-ψf[rint]) - J2loc^(2)/rint^(2)]),{rint,rperi,rapo},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}] // Re];
(*SetSharedFunction[Ω2fromrad];*)

(*Calculation of the derivatives of the intrinsic frequencies with respect to rperi and rapo.
This derivation is made on the grid using the finite differences.
It will only work correctly if rapo ≥ rperi+2ΔrGridMat.
THESE ESTIMATIONS SHOULD BE MADE ONLY ON THE non-largeGRID.*)
(*Finite differences estimation of ∂Ω1/∂rp.*)
DΩ1Drp[rperi_?NumericQ,rapo_?NumericQ,ψf_]:= DΩ1Drp[rperi,rapo,ψf] = (Ω1fromrad[choprad[rperi+ΔrGridMat],rapo,ψf] - Ω1fromrad[choprad[rperi-ΔrGridMat],rapo,ψf])/(2 ΔrGridMat);
(*Finite differences estimation of ∂Ω1/∂ra.*)
DΩ1Dra[rperi_?NumericQ,rapo_?NumericQ,ψf_]:= DΩ1Dra[rperi,rapo,ψf] = (Ω1fromrad[rperi,choprad[rapo+ΔrGridMat],ψf] - Ω1fromrad[rperi,choprad[rapo-ΔrGridMat],ψf])/(2 ΔrGridMat);
(*Finite differences estimation of ∂Ω2/∂rp.*)
DΩ2Drp[rperi_?NumericQ,rapo_?NumericQ,ψf_]:= DΩ2Drp[rperi,rapo,ψf] = (Ω2fromrad[choprad[rperi+ΔrGridMat],rapo,ψf] - Ω2fromrad[choprad[rperi-ΔrGridMat],rapo,ψf])/(2 ΔrGridMat);
(*Finite differences estimation of ∂Ω2/∂ra.*)
DΩ2Dra[rperi_?NumericQ,rapo_?NumericQ,ψf_]:= DΩ2Dra[rperi,rapo,ψf] = (Ω2fromrad[rperi,choprad[rapo+ΔrGridMat],ψf] - Ω2fromrad[rperi,choprad[rapo-ΔrGridMat],ψf])/(2 ΔrGridMat);

CONSTRUCTION OF 2D BASIS ELEMENTS - KALNAJS (1976)

(*The basis elements are given by Kalnajs (1976) and copied in Earn&Sellwood (1995).
Each basis element depends on two index : l azimuthal number (l≥0), n radial number (n≥0). DIFFERENT THAN THE 3D CASE.
The polar coordinates are (R,ϕ).
To construct this basis, there are two shared parameters : RKA (scale radius) and kKA (index).*)

(*Definition of the potential elements*)
ψbasisKA[l_,n_,r_?NumericQ,ϕ_?NumericQ] := Exp[I l ϕ] UbasisKA[l,n,r];

(*Radial component of the potential elements*)
UbasisKA[l_,n_,r_?NumericQ] := - (Sqrt[G/RKA]) PCoeffKA[kKA,l,n] (r/RKA)^(l)Sum[
αCoeffKA[kKA,l,n,isum,jsum] (r/RKA)^(2 isum + 2 jsum)
,{isum,0,kKA,1},{jsum,0,n,1}];

(*Definition of the surface density elements*)
ρbasisKA[l_,n_,r_?NumericQ,ϕ_?NumericQ] :=Exp[I l ϕ] DbasisKA[l,n,r];

(*Radial component of the density elements*)
DbasisKA[l_,n_,r_?NumericQ] := (((-1)^(n))/(Sqrt[G] RKA^(3/2))) SCoeffKA[kKA,l,n] (1 - (r/RKA)^(2))^(kKA - 1/2) (r/RKA)^(l) Sum[
βCoeffKA[kKA,l,n,jsum] (1 - (r/RKA)^(2))^(jsum)
,{jsum,0,n,1}];

(*Definition of the various coefficients needed to define the basis elements.
All these coefficients are memoized because they are expensive to compute.*)
(*Pre-coefficient for the potential elements.*)
PCoeffKA[k_,l_,n_]:= PCoeffKA[k,l,n] =Sqrt[((2 kKA+l+2 n + (1/2)) Gamma[2 kKA+l+n+(1/2)] Gamma[l+n+(1/2)])/(Gamma[2 kKA+n+1] (Gamma[l+1])^(2) Gamma[n+1])];
(*Inside coefficients for the potential elements.*)
αCoeffKA[k_,l_,n_,i_,j_] := αCoeffKA[k,l,n,i,j] = (Pochhammer[-k,i] Pochhammer[l+(1/2),i] Pochhammer[2 k+l+n+(1/2),j] Pochhammer[i+l+(1/2),j] Pochhammer[-n,j])/(Pochhammer[l+1,i] Pochhammer[1,i] Pochhammer[l+i+1,j] Pochhammer[l+(1/2),j] Pochhammer[1,j]);
(*Pre-coefficient for the surface density elements.*)
SCoeffKA[k_,l_,n_] := SCoeffKA[k,l,n] = ((Gamma[kKA+1])/(π Gamma[2 kKA+1] Gamma[kKA+(1/2)])) Sqrt[((2 kKA+l+2 n+(1/2)) Gamma[2 kKA+n+1] Gamma[2 kKA+l+n+(1/2)])/(Gamma[l+n+(1/2)] Gamma[n+1])];
(*Inside coefficients for the surface density elements.*)
βCoeffKA[k_,l_,n_,j_] := βCoeffKA[k,l,n,j] = (Pochhammer[2 k+l+n+(1/2),j] Pochhammer[k+1,j] Pochhammer[-n,j])/(Pochhammer[2 k+1,j] Pochhammer[k+(1/2),j] Pochhammer[1,j]);

2D ANGLES FOURIER TRANSFORM

(*This part of the notebook aims at computing the Fourier transform with respect to the angles of the 2D basis elements.*)
(*The actions always correspond to (J1, J2) = (Jr, L , Lz).
The Fourier index are always given by \bm{m} = (m1, m2).
We have an exact mapping between the peri/apocenters (rperi,rapo) and the actions (J1,J2).*)
(*DEPENDS ON THE POTENTIAL ψf.*)
(*DEPENDS ON THE CHOSEN BASIS Ubasis.
The Kronecker delta has been placed such that if m≠m3, the expensive evaluation is not made.*)
ψbasisFT2D[l_,n_,m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_, Ubasis_] :=
Which[(l≠m2),0,
True,
 WMat2D[m1,l,m2,n,rperi,rapo,ψf,Ubasis]];


(*Definition of the 2D matrix WMat*)
(*DEPENDS ON THE POTENTIAL ψ0*)
(*DEPENDS ON THE CHOSEN BASIS Ubasis.
We take the real part in order not to have any complex components.
One can note that WMat[-m1,-m2,..] = WMat[m1,m2], so that we have to effectively only compute one half of the terms.
If m1>0 or (m1==0 && m2≥0), we compute the term WMat[m1,m2]. Otherwise, we request WMat[-m1,-m2].
This function is the most expensive to compute. The values of this function will be evaluated in parallel. Then known values will be given to it. The function is not shared among the processors.
ATTENTION : We assume that the radial component of the basis elements is always real, so that we do not have any conjugate.
ATTENTION : Because we want to compute WMat in parallel, we had to set $DistributedContexts = Automatic.*)
WMat2D[m1_,l_,m2_,n_,rperi_?NumericQ,rapo_?NumericQ,ψf_, Ubasis_] :=WMat2D[m1,l,m2,n,rperi,rapo,ψf,Ubasis]=
Which[(m1 > 0) || (m1==0 && m2 ≥ 0),
Block[{Eloc,J2loc,Ω1loc,Ω2loc},
Eloc=Efromrad[rperi,rapo,ψf];
J2loc = J2fromrad[rperi,rapo,ψf];
Ω1loc =Ω1fromrad[rperi,rapo,ψf];
Ω2loc=Ω2fromrad[rperi,rapo,ψf]; 1/(π) NIntegrate[rtodθ1dr[Eloc,J2loc,Ω1loc,rint,ψf] Cos[m1 rtoθ1[rperi,Eloc,J2loc,Ω1loc,rint,ψf] + m2 rtoθ2minusψ[rperi,Eloc,J2loc,Ω2loc,rint,ψf]] Ubasis[l,n,rint] ,{rint,rperi,rapo},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}]//Re ],
True,
WMat2D[-m1,l,-m2,n,rperi,rapo,ψf,Ubasis]];
(*SetSharedFunction[WMat];*)

(*Jacobian dθ1/dr(r) used in the estimation of the matrix WMat*)
(*DEPENDS ON THE POTENTIAL ψ0.
We take the real part in order not to have any complex components.*)
rtodθ1dr[Ef_?NumericQ,J2_?NumericQ,Ω1_?NumericQ,r_?NumericQ,ψf_] := (Ω1)/(Sqrt[2(Ef - ψf[r] ) - J2^(2)/r^(2)])//Re;

(*Mapping from r to θ1(r)*)
(*DEPENDS ON THE POTENTIAL ψ0.
We take the real part in order not to have any complex components.*)
rtoθ1[rperi_?NumericQ,Ef_?NumericQ,J2_?NumericQ,Ω1_?NumericQ,r_?NumericQ,ψf_] := NIntegrate[(Ω1)/(Sqrt[2(Ef - ψf[rint]) - J2^(2)/rint^(2)]),{rint,rperi,r},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}]//Re;

(*Mapping from r to θ2minusψ(r)*)
(*DEPENDS ON THE POTENTIAL ψ0.
We take the real part in order not to have any complex components.*)
rtoθ2minusψ[rperi_?NumericQ,Ef_?NumericQ,J2_?NumericQ,Ω2_?NumericQ,r_?NumericQ,ψf_] := NIntegrate[(Ω2 - J2/rint^(2))/(Sqrt[2(Ef - ψf[rint])-J2^(2)/rint^(2)]),{rint,rperi,r},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}]//Re;

2D CALCULATION OF THE RESPONSE MATRIX

(*This part of the notebook aims at computing the Fourier transform with respect to the angles of the 2D basis elements.
Some of the identical functions used here are defined in the section for the 2D response matrix.*)
(*To estimate the response matrix, we will perform the integration in the (rp,ra) space, where all the quantities associated to the actions are easier to compute.
In this space, we will consider a grid of values (rpi, rai), arond which linear approximations, so that Riemann sum formula may be used.*) 


(*Calculation of the response matrix coefficients.
It is made on the grid of pericenters and apocenters.*)
(*DEPENDS ON THE POTENTIAL ψ0*)
(*DEPENDS ON THE CHOSEN BASIS Ubasis*)
(*DEPENDS ON THE BACKGROUND DF : DFEL = DFEL(E,L)*)
(*The sum on m1 has to be made on a finite interval. Note the summation boundaries in this expression.
The Kronecker deltas in (lp,lq) and (mp,mq) have been put in a Which, so that in this situation one does not need to evaluate ZMatCoeff2D.*)

RepMatCoeff2D[lp_,np_,lq_,nq_,ω_?NumericQ,ηf_?NumericQ,ψf_,Ubasis_,DFEL_] := RepMatCoeff2D[lp,np,lq,nq,ω,ηf,ψf,Ubasis,DFEL]=
Which[(lp≠lq),0,
True,
ZMatCoeff2D[lp,np,nq,ω,ηf,ψf,Ubasis,DFEL]];
(*SetSharedFunction[RepMatCoeff2D];*)


(*ZMatCoeff2D is the coefficient which only depend on lp, np, nq.
It allows us to reduce as much as possible the computation required to estimate the full diagonal-by bloc response matrix RepMatCoeff2D.
We also have to get rid of the case for m1sum = m2sum for which the numerator is exactly 0.
We use the symmetry property of ZMat2D with respect to the radial n index.
Indeed, we have ZMatCoeff2D[np,nq] = ZMatCoeff2D[nq,np].
If (np ≤ nq), we compute ZMatCoeff2D[np,nq], otherwise, we compute ZMatCoeff2D[nq,np].
We make two successive loops. The first one is on iGridsum, so that one only has to determine once the values of rperiloc and rapoloc.
In the 2D case, one always has m2sum = lp. We kept the useless sum on m2sum to have the same structure for the function as in the 3D case.*)
ZMatCoeff2D[lp_,np_,nq_,ω_?NumericQ,ηf_?NumericQ,ψf_,Ubasis_,DFEL_] := ZMatCoeff2D[lp,np,nq,ω,ηf,ψf,Ubasis,DFEL]=
Which[np ≤ nq,
Block[{rperiloc,rapoloc},
Sum[
rperiloc = rperifromGrid[iGridsum];
rapoloc = rapofromGrid[iGridsum];
Sum[
Which[(m1sum == 0) && (m2sum == 0),0,
True,ℵMat[
agMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,ψf,Ubasis,DFEL],
bgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,ψf,Ubasis,DFEL],
cgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,ψf,Ubasis,DFEL],
ahMat2D[m1sum,m2sum,ω,rperiloc,rapoloc,ψf],
bhMat2D[m1sum,m2sum,rperiloc,rapoloc,ψf],
chMat2D[m1sum,m2sum,rperiloc,rapoloc,ψf],
ηf,
ΔrGridMat]],{m1sum,-m1sumMatbound,m1sumMatbound,1},{m2sum,lp,lp,1}],{iGridsum,1,nbelementsGrid}]
],True,
ZMatCoeff2D[lp,nq,np,ω,ηf,ψf,Ubasis,DFEL]];
(*SetSharedFunction[ZMatCoeff2D];*)


(*Same as ZMatCoeff2D.
ALL THE COMPUTATIONS ARE MADE USING ℵMatAna.*)
ZMatCoeffAna2D[lp_,np_,nq_,ω_?NumericQ,ηf_?NumericQ,ψf_,Ubasis_,DFEL_] := ZMatCoeffAna2D[lp,np,nq,ω,ηf,ψf,Ubasis,DFEL]=
Which[np ≤ nq,
Block[{rperiloc,rapoloc},
Sum[
rperiloc = rperifromGrid[iGridsum];
rapoloc = rapofromGrid[iGridsum];
Sum[
Which[(m1sum == 0) && (m2sum == 0),0,
True,ℵMatAna[
agMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,ψf,Ubasis,DFEL],
bgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,ψf,Ubasis,DFEL],
cgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,ψf,Ubasis,DFEL],
ahMat2D[m1sum,m2sum,ω,rperiloc,rapoloc,ψf],
bhMat2D[m1sum,m2sum,rperiloc,rapoloc,ψf],
chMat2D[m1sum,m2sum,rperiloc,rapoloc,ψf],
ηf,
ΔrGridMat]],{m1sum,-m1sumMatbound,m1sumMatbound,1},{m2sum,lp,lp,1}],{iGridsum,1,nbelementsGrid}]
],True,
ZMatCoeffAna2D[lp,nq,np,ω,ηf,ψf,Ubasis,DFEL]];


(*Definition of the coefficients of limited development to calculate the response matrix
To calculate these coefficients, we generate the grid for the functions g and h, and then by finite differences we can determine a, b and c.
The finite-difference formula that we use is : fprime(r0) = (f(r0+Δr)-f(r0-Δr))/(2 Δr).
We note that the only function which is not memoized is ahMat2D, because it has an ω-dependence.
We use choprad to get sanitized values for rperi and rapo which will be in listGridlarge.
Some index are useless in gMat2D, agMat2D, bgMat2D, cgMat2D. I kept it so that the fonctions have the same structure as the ones used in the 3D case.*)
agMat2D[lp_,np_,nq_,m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_,Ubasis_,DFEL_]:=agMat2D[lp,np,nq,m1,m2,rperi,rapo,ψf,Ubasis,DFEL]=gMat2D[lp,np,nq,m1,m2,rperi,rapo,ψf,Ubasis,DFEL];

bgMat2D[lp_,np_,nq_,m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_,Ubasis_,DFEL_] :=bgMat2D[lp,np,nq,m1,m2,rperi,rapo,ψf,Ubasis,DFEL]=(gMat2D[lp,np,nq,m1,m2,choprad[rperi+ΔrGridMat],rapo,ψf,Ubasis,DFEL] - gMat2D[lp,np,nq,m1,m2,choprad[rperi-ΔrGridMat],rapo,ψf,Ubasis,DFEL])/(2 ΔrGridMat);

cgMat2D[lp_,np_,nq_,m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_,Ubasis_,DFEL_] :=cgMat2D[lp,np,nq,m1,m2,rperi,rapo,ψf,Ubasis,DFEL]=(gMat2D[lp,np,nq,m1,m2,rperi,choprad[rapo+ΔrGridMat],ψf,Ubasis,DFEL]-gMat2D[lp,np,nq,m1,m2,rperi,choprad[rapo-ΔrGridMat],ψf,Ubasis,DFEL])/(2 ΔrGridMat);

ahMat2D[m1_,m2_,ω_?NumericQ,rperi_?NumericQ,rapo_?NumericQ,ψf_] := 
ahMat2D[m1,m2,ω,rperi,rapo,ψf]=hMat2D[m1,m2,ω,rperi,rapo,ψf];

(*This function is memoized and does not depend on ω.
It requires to evaluate hMat2D in ω=0.*)
bhMat2D[m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_] :=  bhMat2D[m1,m2,rperi,rapo,ψf]=(hMat2D[m1,m2,0.,choprad[rperi+ΔrGridMat],rapo,ψf] - hMat2D[m1,m2,0.,choprad[rperi-ΔrGridMat],rapo,ψf])/(2 ΔrGridMat);

(*This function is memoized and does not depend on ω.
It requires to evaluate hMat2D in ω=0.*)
chMat2D[m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_] := 
chMat2D[m1,m2,rperi,rapo,ψf]=(hMat2D[m1,m2,0.,rperi,choprad[rapo+ΔrGridMat],ψf] - hMat2D[m1,m2,0.,rperi,choprad[rapo-ΔrGridMat],ψf])/(2 ΔrGridMat);

(*Definition of the function g and h, which are at the heart of the calculation of RepMat2D.
These functions are memoized so that the calculation of the coefficients a,b and c, via the finite differences method can be made efficiently.
ATTENTION : We dropped the conjugate on WMat2D[np] because it is always real.*)
gMat2D[lp_,np_,nq_,m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_,Ubasis_,DFEL_] :=
gMat2D[lp,np,nq,m1,m2,rperi,rapo,ψf,Ubasis,DFEL]=(2 π)^(2) JacELtorpra[rperi,rapo,ψf] (1/Ω1fromrad[rperi,rapo,ψf]) DFGrad[m1,m2,rperi,rapo,ψf,DFEL] WMat2D[m1,lp,m2,np,rperi,rapo,ψf,Ubasis] WMat2D[m1,lp,m2,nq,rperi,rapo,ψf,Ubasis]  ;
(*SetSharedFunction[gMat2D];*)

hMat2D[m1_,m2_,ω_?NumericQ,rperi_?NumericQ,rapo_?NumericQ,ψf_] :=
hMat2D[m1,m2,ω,rperi,rapo,ψf]=ω - m1 Ω1fromrad[rperi,rapo,ψf] -m2 Ω2fromrad[rperi,rapo,ψf];
(*SetSharedFunction[hMat2D];*)

(*Jacobian of the transformation (E,L)->(rp,ra), necessary to calculate g
DEDrp is \partial E/ \partial rp.
DEDra is \partial E/ \partial ra.
DLDrp is \partial L/ \partial rp.
DLDra is \partial L/ \partial ra.*) 
JacELtorpra[rperi_?NumericQ,rapo_?NumericQ,ψf_] := JacELtorpra[rperi,rapo,ψf]=Abs[DEDrp[rperi,rapo,ψf] DLDra[rperi,rapo,ψf] - DLDrp[rperi,rapo,ψf] DEDra[rperi,rapo,ψf]];
(*SetSharedFunction[JacELtorpra];*)

(*Calculation of \partial E/ \partial rp.*)
DEDrp[rperi_?NumericQ,rapo_?NumericQ,ψf_]:=Block[{DψDrploc},
DψDrploc=D[ψf[rder],{rder,1}] /. rder -> rperi;
(rperi (2 rapo^(2) (ψf[rapo] - ψf[rperi]) + rperi (rperi^(2) - rapo^(2)) DψDrploc ))/((rapo^(2) - rperi^(2))^(2))];

(*Calculation of \partial E/ \partial ra.*)
DEDra[rperi_?NumericQ,rapo_?NumericQ,ψf_]:= Block[{DψDraloc},
DψDraloc=D[ψf[rder],{rder,1}] /. rder -> rapo;
(rapo (2 rperi^(2) (ψf[rperi] - ψf[rapo]) + rapo (rapo^(2) - rperi^(2)) DψDraloc))/((rapo^(2) - rperi^(2))^(2))];

(*Calculation of \partial L/ \partial rp.*)
DLDrp[rperi_?NumericQ,rapo_?NumericQ,ψf_] := Block[{DψDrploc},
DψDrploc = D[ψf[rder],{rder,1}] /. rder -> rperi;
-(rapo^(2) rperi (2 rapo^(2) (ψf[rperi] - ψf[rapo]) + DψDrploc(rapo^(2) rperi - rperi^(3))))/(Sqrt[2] (rapo - rperi)^(2) (rapo + rperi)^(2) Sqrt[(rapo^(2) rperi^(2) (ψf[rapo] - ψf[rperi]))/(rapo^(2) - rperi^(2))])];

(*Calculation of \partial L/ \partial ra.*)
DLDra[rperi_?NumericQ,rapo_?NumericQ,ψf_] := Block[{DψDraloc},
DψDraloc = D[ψf[rder],{rder,1}] /. rder -> rapo;
(rapo rperi^(2) (2 rperi^(2) (ψf[rperi] - ψf[rapo]) + DψDraloc (rapo^(3) - rapo rperi^(2))))/(Sqrt[2] (rapo - rperi)^(2) (rapo+rperi)^(2) Sqrt[(rapo^(2) rperi^(2) (ψf[rapo] - ψf[rperi]))/(rapo^(2) - rperi^(2))])];


(*Gradient of the distribution function used for the calculation of RepMat2D
DFDE is \partial F/ \partial E.
DFDL is \partial F/ \partial L.
DFEL is F (E,L).*)
DFGrad[m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,ψf_,DFEL_] := DFGrad[m1,m2,rperi,rapo,ψf,DFEL]=Block[{Eloc,Lloc,Ω1loc,Ω2loc,DFDEloc},
Eloc = Efromrad[rperi,rapo,ψf];
Lloc = J2fromrad[rperi,rapo,ψf];
Ω1loc= Ω1fromrad[rperi,rapo,ψf];
Ω2loc = Ω2fromrad[rperi,rapo,ψf];
DFDEloc = DFDE[Eloc,Lloc,DFEL];
m1 Ω1loc DFDEloc + m2 (Ω2loc DFDEloc + DFDL[Eloc,Lloc,DFEL])]
(*SetSharedFunction[DFGrad];*)

(*Calculation of \partial F/ \partial E.
This function is memoized for faster access.
In order not to try to compute for E ≥ 0, we add a scale argument, so that ND does not go too far.
HERE, WE SHOULD GET RID OF THE NUMERICAL DERIVATIVES.*)
DFDE[Ef_?NumericQ,Lf_?NumericQ,DFEL_] := DFDE[Ef,Lf,DFEL]=ND[DFEL[Eder,Lf],Eder,Ef,Scale -> 0.01] ;
(*SetSharedFunction[DFDE];*)

(*Calculation of \partial F/ \partial L.
This function is memoized for faster access.
In order to always have L > 0, we add a scale argument, so that ND does not go too far.
HERE, WE SHOULD GET RID OF THE NUMERICAL DERIVATIVES.*)
DFDL[Ef_?NumericQ, Lf_?NumericQ,DFEL_] := DFDL[Ef,Lf,DFEL]=ND[DFEL[Ef,Lder],Lder,Lf, Scale-> 0.01] ;
(*SetSharedFunction[DFDL];*)


(*Construction of the full response matrix.
It is at this stage that we might use ParallelMap, and RepMatCoeff2D uses Sum, instead of ParallelSum.
In order to only compute one half of the coefficients, we rely on the symmetry of M, so that M^t = M.*)
RepMat2D[ω_?NumericQ,ηf_?NumericQ,ψf_,Ubasis_,DFEL_] := Block[{halfMloc,Mloc},
halfMloc=Table[RepMatCoeff2D[lfromlistbasis2D[ptab],nfromlistbasis2D[ptab],lfromlistbasis2D[qtab],nfromlistbasis2D[qtab],ω,ηf,ψf,Ubasis,DFEL],{ptab,1,nbelementsbasis2D,1},{qtab,ptab,nbelementsbasis2D,1}];
Mloc = #+Transpose@UpperTriangularize[#,1]&@PadLeft@halfMloc;
Mloc
];

(*Construction of the Z response matrix used to determine the modes of a system.
It depends on :
+ lmode : the mode number considered.
+ nmode : the number of radial elements considered. ATTENTION IT STARTS AT 0.
+ ω : the real part frequency of evaluation.
+ η : the imaginary part of the frequency of evaluation.
In order not to have to use the list of the basis elements considered, we only have to provide the number of radial elements considered (nmode).
As a consequence, ZMat2D is a matrix of size nmode x nmode.
The radial numbers np and nq belong to the interval 1 ≤ np ≤ nmode.
*)
ZMat2D[lmodef_,nmodef_,ωf_?NumericQ,ηf_?NumericQ,ψf_,Ubasis_,DFEL_] := Table[ZMatCoeff2D[lmodef,nptab,nqtab,ωf,ηf,ψf,Ubasis,DFEL],{nptab,0,nmodef,1},{nqtab,0,nmodef,1}];

(*Same as ZMat2D.
ALL THE COMPUTATIONS ARE MADE USING ℵMatAna.
*)
ZMatAna2D[lmodef_,nmodef_,ωf_?NumericQ,ηf_?NumericQ,ψf_,Ubasis_,DFEL_] := Table[ZMatCoeffAna2D[lmodef,nptab,nqtab,ωf,ηf,ψf,Ubasis,DFEL],{nptab,0,nmodef,1},{nqtab,0,nmodef,1}];


(*We define the determinant to determine the presence of unstable modes.
DetZMat2D[..] = det(I - Z[..]).
The first evaluation of DetZMat2D is long because the intrinsic frequencies have to determined on the grid to be memoized on the master kernel.
With the function DetZMat2D, one can make Nyquist diagram.
The parameters of the functions are :
+ lmode : the mode number considered.
+ nmode : the number of radial elements considered, i.e. the size of the matrices.
+ ω : the real part of the frequency of evaluation.
+ η : the imaginary part of the frequency of evaluation.*)
DetZMat2D[lmodef_,nmodef_,ωf_?NumericQ,ηf_?NumericQ,ψf_,Ubasis_,DFEL_]:=
Det[IdentityMatrix[nmodef+1] -  ZMat2D[lmodef,nmodef,ωf,ηf,ψf,Ubasis,DFEL]];

SUB - REGION INTEGRATION FOR THE RESPONSE MATRIX

(*In this part of the notebook, we build up the wrapped analytical function, which will allow to calculate each contribution from the subregions to the response matrix.*)
(*THE FUNCTIONS IN THIS SECTION HAVE BEEN IN A REVERSE ORDER IN ORDER TO BE CORRECTLY COMPILED.*)


(*New version of the sub-region matrix integration.
ALL THE CASES ARE HANDLED IN THE SAME WAY.*)
ℵMatD = Compile[{{bf,_Real},{cf,_Real},{ef,_Real},{ff,_Real},{ηf,_Real}},
(1/(16 (ef^2) (ff^2)))(-I (ff (bf ((ef^2)-(((-2+ff-2 I ηf))^2))+4 ef (2+ef+2 I ηf))-cf ef (((2+ef+2 I ηf))^2)) (π-2 ArcTan[((2+ef-ff)/(2 ηf))])-I (-ff (bf (-(ef^2)+(((2+ff+2 I ηf))^2))+4 ef (-2+ef-2 I ηf))-cf ef (((-2+ef-2 I ηf))^2)) (π-2 ArcTan[((2-ef+ff)/(2 ηf))])+I (ff (bf ((ef^2)-(((-2+ff-2 I ηf))^2))-4 ef (-2+ef-2 I ηf))-cf ef (((-2+ef-2 I ηf))^2)) (π+2 ArcTan[((-2+ef+ff)/(2 ηf))])+I (ff (bf ((ef^2)-(((2+ff+2 I ηf))^2))+4 ef (2+ef+2 I ηf))-cf ef (((2+ef+2 I ηf))^2)) (π-2 ArcTan[((2+ef+ff)/(2 ηf))])+ff (ff (8 ef+bf (-4+2 ef+ff-4 I ηf))+cf ef (-4+2 ef-ff-4 I ηf)+2 (-4+cf) ef ff Log[(1/2) (2-ef-ff+2 I ηf)])-ff (ff (8 ef+bf (-4-2 ef+ff-4 I ηf))-cf ef (4+2 ef+ff+4 I ηf)+2 (-4+cf) ef ff Log[(1/2) (2+ef-ff+2 I ηf)])-ff (ff (-8 ef+bf (4-2 ef+ff+4 I ηf))-cf ef (-4+2 ef+ff-4 I ηf)+2 (4+cf) ef ff Log[(1/2) (2-ef+ff+2 I ηf)])+ff (ff (-8 ef+bf (4+2 ef+ff+4 I ηf))+cf ef (4+2 ef-ff+4 I ηf)+2 (4+cf) ef ff Log[(1/2) (2+ef+ff+2 I ηf)])+(ff (bf (-(ef^2)+(((-2+ff-2 I ηf))^2))-4 ef (2+ef+2 I ηf))+cf ef (((2+ef+2 I ηf))^2)) Log[(1/4) (4+(ef^2)-2 ef (-2+ff)-4 ff+(ff^2)+4 ((ηf)^2))]-(ff (bf (-(ef^2)+(((-2+ff-2 I ηf))^2))+4 ef (-2+ef-2 I ηf))+cf ef (((-2+ef-2 I ηf))^2)) Log[(1/4) (4+(ef^2)+2 ef (-2+ff)-4 ff+(ff^2)+4 ((ηf)^2))]+(ff (bf (-(ef^2)+(((2+ff+2 I ηf))^2))+4 ef (-2+ef-2 I ηf))+cf ef (((-2+ef-2 I ηf))^2)) Log[(1/4) (4+(ef^2)+4 ff+(ff^2)-2 ef (2+ff)+4 ((ηf)^2))]-(ff (bf (-(ef^2)+(((2+ff+2 I ηf))^2))-4 ef (2+ef+2 I ηf))+cf ef (((2+ef+2 I ηf))^2)) Log[(1/4) (4+(ef^2)+4 ff+(ff^2)+2 ef (2+ff)+4 ((ηf)^2))]),CompilationTarget->"C",CompilationOptions->{"ExpressionOptimization" -> True,"InlineCompiledFunctions"->True,"InlineExternalDefinitions"->True},RuntimeOptions->"Speed"];


(*Definition of the analytical function ℵMat which is given by
ℵMat[agf,bgf,cgf,ahf,bhf,chf,ηf,Δf] = int_{-Δf/2}^{Δf/2} dxp int_{-Δf/2}^{Δf/2} dxa (agf + bgf xp + cgf xa)/(ahf + bhf xp + chf xa + ηf).
This function, despite the external call has been compiled.
In the case where agf or ahf are zero, this function breaks down.*)
ℵMat=Compile[{{agf,_Real},{bgf,_Real},{cgf,_Real},{ahf,_Real},{bhf,_Real},{chf,_Real},{ηf,_Real},{Δf,_Real}},
Block[{bloc,cloc,eloc,floc,ηloc},
bloc = (bgf Δf)/(agf);
cloc = (cgf Δf)/(agf);
eloc = (bhf Δf)/(ahf);
floc = (chf Δf)/(ahf);
ηloc = (ηf)/(ahf);
(agf)/(ahf) (Δf)^(2) ℵMatD[bloc,cloc,eloc,floc,ηloc]],
{{ℵMatD[_,_,_,_,_],_Complex}},CompilationTarget->"C",CompilationOptions->{"ExpressionOptimization" -> True,"InlineCompiledFunctions"->True,"InlineExternalDefinitions"->True},RuntimeOptions->"Speed"];


(*Analytical version of ℵMat.*)
ℵMatAna=Compile[{{agf,_Real},{bgf,_Real},{cgf,_Real},{ahf,_Real},{bhf,_Real},{chf,_Real},{ηf,_Real},{Δf,_Real}},
Block[{bloc,cloc,eloc,floc,ηloc},
bloc = (bgf Δf)/(agf);
cloc = (cgf Δf)/(agf);
eloc = (bhf Δf)/(ahf);
floc = (chf Δf)/(ahf);
ηloc = (ηf)/(ahf);
(agf)/(ahf) (Δf)^(2) ℵMatDAna[bloc,cloc,eloc,floc,ηloc]],
{{ℵMatDAna[_,_,_,_,_],_Complex}},CompilationTarget->"C",CompilationOptions->{"ExpressionOptimization" -> True,"InlineCompiledFunctions"->True,"InlineExternalDefinitions"->True},RuntimeOptions->"Speed"];

(*Definition of the exact analytical expression.
This function has been used to check the sanity of the approximated expressions.*)
ℵMatDAna[bf_?NumericQ,cf_?NumericQ,ef_?NumericQ,ff_?NumericQ,ηf_?NumericQ] := NIntegrate[(1 + bf xint + cf yint)/(1 + ef xint + ff yint + I ηf),{xint,-0.5,0.5},{yint,-0.5,0.5},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}];

2D PREPARATION OF THE CALCULATIONS

(*All the functions present in this block intend to prepare the parallel calculations of the expensive functions, which will then be memoized.*)

(*We perform once for all the evaluation of all the WMat2D coefficients required.
We first determine the number of evaluation of WMat2D which will be needed to estimate the time of calculation. We divide this number by 2 because of parity.
One evaluation of WMat2D takes at the moment a bit less than 1s.
ATTENTION nptab starts at 0.*)
GetnbevalsWMat2D[lmodef_,nmodef_,m1boundf_]:= Block[{nbevalsWMat2Dloc},nbevalsWMat2Dloc = Table[If[(m1tab==0)&&(m2tab==0),0,1],{m1tab,-m1boundf,m1boundf,1},{ltab,lmodef,lmodef,1},{m2tab,ltab,ltab,1},{nptab,0,nmodef},{iGridtab,1,nbelementsGridlarge}] // Flatten // Total // #&;
Print["nbevalsWMat2D = ",nbevalsWMat2Dloc]];


(*The aim of this function is do once for all the computations needed to estimate the WMat2D coefficients.
The parameters needed to determine the density of the grid in WMat2D are :
+ m1sumMatbound : bounds.
+ lmode : l-number of the mode.
+ nmode : number of radial elements 1 ≤ np ≤ nmode
+ All the parameters used to generate the Gridlarge of rperi and rapo.
We perform once for all the evaluation of all the WMat2D coefficients required.
It is important that the function WMat2D is set as a SetShared Function because the values of WMat2D will be computed in parallel, and the result need to be returned to stored in the main kernel.
WMat2D depends on the following quantities : 
+ m1 : One index of the sum in RepMat2D expression.
+ l : Index of the basis elements. If one considers one specific mode, l = lmode = cst.
+ m2 : One index of the sum in RepMat2D expression.
+ i : Index in the (rperi,rapo) grid.
+ n : Index of the radial mode.
In order not to slow down the computation by sharing too many functions, while still using the parallel tools, one can proceed as follows. Each sub-kernel will be given one index i, i.e. one position (rperi, rapo) in action-space. Each of this sub-kernel will keep its local collection of the quantities depending only on i, i.e. E,L,Jr,Ω1,Ω2. Then each of the sub-kernel will perform the loop on remaining index.
*)
(*We first have to manually build up the values of WMat2D which are computed in parallel.
Each element of the table is of the form {m1tab,ltab,m2tab,nptab,rperi,rapo,ψ,UbasisHE}.
The table has the form of a list (cf Flatten).
We also have to remove the Null elements.*)
(*Once all the required values of WMat2D have been computed and stored in a table, they can be loaded in the definition of WMat2D.
One should pay attention to the order of the arguments and not get mixed up.
ATTENTION nptab starts at 0.*)
(*ATTENTION THIS FUNCTION USES PARALLELTABLE -> We had to set $DistributedContexts = Automatic.*)
PrepareWMat2D[lmodef_,nmodef_,m1boundf_,ψf_,Ubasisf_]:=Block[{tabValWMat2D,vloc},
tabValWMat2D = ParallelTable[
Table[
If[(m1tab==0)&&(m2tab==0),Null,
{m1tab,ltab,m2tab,nptab,rperifromGridlarge[iGridtab],rapofromGridlarge[iGridtab],ψf,Ubasisf,WMat2D[m1tab,ltab,m2tab,nptab,rperifromGridlarge[iGridtab],rapofromGridlarge[iGridtab],ψf,Ubasisf]}],
{m1tab,-m1sumMatbound,m1sumMatbound,1},{ltab,lmodef,lmodef,1},{m2tab,ltab,ltab,1},{nptab,0,nmodef}],{iGridtab,1,nbelementsGridlarge}]//Flatten[#,4]& // DeleteCases[#,Null]&;

Do[vloc = tabValWMat2D[[ind]];
WMat2D[vloc[[1]],vloc[[2]],vloc[[3]],vloc[[4]],vloc[[5]],vloc[[6]],vloc[[7]],vloc[[8]]]= vloc[[9]],{ind,1,tabValWMat2D//Length}]
];

(*After having calculated all the needed values of WMat2D, we dump them in a data file, so that they may be reused.
ATTENTION : This dump will overwrite any existing file with the same name.
This is made using the Downvalues, which are dumped in a data file.
Then, one can load these DownValues.
For some strange reason, inside the function the shortened syntax >> does not work, whereas Put[] works.
In order not to have bad characters, 
adressfiles has to be changed depending if you are in local or on remote servers.*)
DumpWMat2D[lmodef_,nmodef_,m1boundf_,rGridminf_,rGridmaxf_,ΔrGridf_,ψf_,Ubasisf_] := Block[{namefileloc},
namefileloc = adressfiles<>"DumpWMat2D_lmode_"<>ToString[lmodef]<>"_nmode_"<>ToString[nmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>"_basis_"<>ToString[FullForm[Ubasisf]]<>".txt";
Put[DownValues[WMat2D],namefileloc];
(*DownValues[WMat2D] >> namefileloc;*)
];

(*Function which calculates in parallel the intrinsic frequencies to allow a faster access afterwards.
This function computes on listGridlarge : E,J2,J1,Ω1,Ω2.
The only argument needed is the background potential ψf.*)
PrepareΩ[ψf_]:=Block[{tabΩ,rploc,raloc,ψloc,vloc},
tabΩ = ParallelTable[
rploc = rperifromGridlarge[iGridtab];
raloc =rapofromGridlarge[iGridtab];
{rploc,raloc,ψf,
Efromrad[rploc,raloc,ψf],
J2fromrad[rploc,raloc,ψf],
J1fromrad[rploc,raloc,ψf],
Ω1fromrad[rploc,raloc,ψf],
Ω2fromrad[rploc,raloc,ψf]},{iGridtab,1,nbelementsGridlarge}];

Do[vloc = tabΩ[[ind]];
rploc = vloc[[1]];
raloc = vloc[[2]];
ψloc = vloc[[3]];
Efromrad[rploc,raloc,ψloc]=vloc[[4]];
J2fromrad[rploc,raloc,ψloc]=vloc[[5]];
J1fromrad[rploc,raloc,ψloc]=vloc[[6]];
Ω1fromrad[rploc,raloc,ψloc] =vloc[[7]];
Ω2fromrad[rploc,raloc,ψloc]=vloc[[8]];
,{ind,1,tabΩ//Length}]
];

(*Function which dumps the downvalues of E,J2,J1,Ω1,Ω2 already computed.
ATTENTION : This dump will overwrite any existing file with the same name.
This is made using the Downvalues, which are dumped in a data file.
Then, one can load these DownValues.
For some strange reason, inside the function the shortened syntax >> does not work, whereas Put[] works.
In order not to have bad characters, 
adressfiles has to be changed depending if you are in local or on remote servers.*)
DumpΩ[rGridminf_,rGridmaxf_,ΔrGridf_,ψf_] := Block[{namefileloc},
(*For E.*)
namefileloc = adressfiles<>"DumpE_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
Put[DownValues[Efromrad],namefileloc];
(*For J2.*)
namefileloc = adressfiles<>"DumpJ2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
Put[DownValues[J2fromrad],namefileloc];
(*For J1.*)
namefileloc = adressfiles<>"DumpJ1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
Put[DownValues[J1fromrad],namefileloc];
(*For Ω1.*)
namefileloc = adressfiles<>"DumpOmega1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
Put[DownValues[Ω1fromrad],namefileloc];
(*For Ω2.*)
namefileloc = adressfiles<>"DumpOmega2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
Put[DownValues[Ω2fromrad],namefileloc];
(*DownValues[Ω2fromrad] >> namefileloc;*)
];

(*This function allows to load previously calculated values of E,J2,J1,Ω1,Ω2.
It depends on the same parameters a DumpΩ.
ATTENTION : In the current implementation, once loaded, it overwrites all previously memoized values of Efromrad,J2fromrad,J1fromrad,Ω1fromrad,Ω2fromrad, and also changes the definition of these functions using the one present in the dump file.
adressfiles has to be changed depending if you are in local or on remote servers.*)
LoadΩ[rGridminf_,rGridmaxf_,ΔrGridf_,ψf_]:= Block[{namefileloc},
(*For E.*)
namefileloc = adressfiles<>"DumpE_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
DownValues[Efromrad] = ReadList[namefileloc][[1]];
(*For J2.*)
namefileloc = adressfiles<>"DumpJ2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
DownValues[J2fromrad] = ReadList[namefileloc][[1]];
(*For J1.*)
namefileloc = adressfiles<>"DumpJ1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
DownValues[J1fromrad] = ReadList[namefileloc][[1]];
(*For Ω1.*)
namefileloc = adressfiles<>"DumpOmega1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
DownValues[Ω1fromrad] = ReadList[namefileloc][[1]];
(*For Ω2.*)
namefileloc = adressfiles<>"DumpOmega2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
DownValues[Ω2fromrad] = ReadList[namefileloc][[1]];
];


(*Function which calculates in parallel the intrinsic frequencies to allow a faster access afterwards.
This function computes on listGridlarge : JacELtorpra.
The only argument needed is the background potential ψf.*)
PrepareJac[ψf_]:=Block[{tabJac,rploc,raloc,ψloc,vloc},
tabJac = ParallelTable[
rploc = rperifromGridlarge[iGridtab];
raloc =rapofromGridlarge[iGridtab];
{rploc,raloc,ψf,
JacELtorpra[rploc,raloc,ψf]},{iGridtab,1,nbelementsGridlarge}];

Do[vloc = tabJac[[ind]];
rploc = vloc[[1]];
raloc = vloc[[2]];
ψloc = vloc[[3]];
JacELtorpra[rploc,raloc,ψloc]=vloc[[4]];
,{ind,1,tabJac//Length}]
];

(*Function which dumps the downvalues of JacELtorpra already computed.
ATTENTION : This dump will overwrite any existing file with the same name.
This is made using the Downvalues, which are dumped in a data file.
Then, one can load these DownValues.
For some strange reason, inside the function the shortened syntax >> does not work, whereas Put[] works.
In order not to have bad characters, 
adressfiles has to be changed depending if you are in local or on remote servers.*)
DumpJac[rGridminf_,rGridmaxf_,ΔrGridf_,ψf_] := Block[{namefileloc},
namefileloc = adressfiles<>"DumpJac_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
Put[DownValues[JacELtorpra],namefileloc];
(*DownValues[JacELtorpra] >> namefileloc;*)
];

(*This function allows to load previously calculated values of JacELtorpra.
It depends on the same parameters a DumpJac.
ATTENTION : In the current implementation, once loaded, it overwrites all previously memoized values of JacELtorpra, and also changes the definition of these functions using the one present in the dump file.
adressfiles has to be changed depending if you are in local or on remote servers.*)
LoadJac[rGridminf_,rGridmaxf_,ΔrGridf_,ψf_]:= Block[{namefileloc,listread},
namefileloc = adressfiles<>"DumpJac_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>".txt";
listread =  ReadList[namefileloc][[1]];
DownValues[JacELtorpra] =listread;
];



(*To prepare the computations of DFGrad in the 2D case.
For a 2D case, we compute a bit too many values, see the loop on m2tab.*)
PrepareDFGrad[lmodef_,m1boundf_,ψf_,DFELf_]:=Block[{tabDFGrad,rperiloc,rapoloc,vloc},
tabDFGrad = ParallelTable[
rperiloc = rperifromGridlarge[iGridtab];
rapoloc = rapofromGridlarge[iGridtab];
Table[
{m1tab,m2tab,rperiloc,rapoloc,ψf,DFELf,DFGrad[m1tab,m2tab,rperiloc,rapoloc,ψf,DFELf]}
,{m1tab,-m1sumMatbound,m1sumMatbound,1},{ltab,lmodef,lmodef,1},{m2tab,-ltab,ltab,2}],{iGridtab,1,nbelementsGridlarge}]//Flatten[#,3]& ;

Do[vloc = tabDFGrad[[ind]];
DFGrad[vloc[[1]],vloc[[2]],vloc[[3]],vloc[[4]],vloc[[5]],vloc[[6]]]= vloc[[7]],{ind,1,tabDFGrad//Length}]
];

(*To dump the already computed values of DFGrad in the 2D case.*)
DumpDFGrad[lmodef_,m1boundf_,rGridminf_,rGridmaxf_,ΔrGridf_,ψf_,DFELf_] := Block[{namefileloc},
namefileloc = adressfiles<>"DumpDFGrad_lmode_"<>ToString[lmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>"_DFEL_"<>ToString[FullForm[DFELf]]<>".txt";
Put[DownValues[DFGrad],namefileloc];
(*DownValues[DFGrad] >> namefileloc;*)
];

(*To load the already computed values of DFGrad in the 2D case.*)
LoadDFGrad[lmodef_,m1boundf_,rGridminf_,rGridmaxf_,ΔrGridf_,ψf_,DFELf_]:= Block[{namefileloc,listread},
namefileloc = adressfiles<>"DumpDFGrad_lmode_"<>ToString[lmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>"_DFEL_"<>ToString[FullForm[DFELf]]<>".txt";
listread =  ReadList[namefileloc][[1]];
DownValues[DFGrad] =listread;
];



(*This function allows to load previously calculated values of WMat2D.
It depends on the same parameters a DumpWMat2D.
ATTENTION : In the current implementation, once loaded, it overwrites all previously memoized values of WMat2D, and also changes the definition of WMat2D using the one present in the dump file.
adressfiles has to be changed depending if you are in local or on remote servers.*)
LoadWMat2D[lmodef_,nmodef_,m1boundf_,rGridminf_,rGridmaxf_,ΔrGridf_,ψf_,Ubasisf_]:= Block[{namefileloc,listread},
namefileloc = adressfiles<>"DumpWMat2D_lmode_"<>ToString[lmodef]<>"_nmode_"<>ToString[nmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>"_basis_"<>ToString[FullForm[Ubasisf]]<>".txt";
listread = ReadList[namefileloc][[1]];
DownValues[WMat2D] = listread;
];

(*Function to dump the values already computed of ZMatCoeff2D.*)
DumpZMatCoeff2D[lmodef_,nmodef_,m1boundf_,rGridminf_,rGridmaxf_,ΔrGridf_,ψf_,Ubasisf_,DFELf_]:=
Block[{namefileloc},
namefileloc = adressfiles<>"DumpZMatCoeff2D_lmode_"<>ToString[lmodef]<>"_nmode_"<>ToString[nmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>"_basis_"<>ToString[FullForm[Ubasisf]]<>"_DFEL_"<>ToString[FullForm[DFELf]]<>".txt";
Put[DownValues[ZMatCoeff2D],namefileloc];
];

(*Function to load the already computed values of ZMatCoeff2D.
All the previously computed values of ZMatCoeff2D are overwritten.*)
LoadZMatCoeff2D[lmodef_,nmodef_,m1boundf_,rGridminf_,rGridmaxf_,ΔrGridf_,ψf_,Ubasisf_,DFELf_]:= Block[{namefileloc,listread},
namefileloc = adressfiles<>"DumpZMatCoeff2D_lmode_"<>ToString[lmodef]<>"_nmode_"<>ToString[nmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[ΔrGridf]<>"_pot_"<>StringDrop[ToString[FullForm[ψf]],1]<>"_basis_"<>ToString[FullForm[Ubasisf]]<>"_DFEL_"<>ToString[FullForm[DFELf]]<>".txt";
listread = ReadList[namefileloc][[1]];
DownValues[ZMatCoeff2D] = listread;
]

End[];
EndPackage[];