(* ::Package:: *)

(* ::Text::Closed:: *)
(*INTRODUCTION AND DESCRIPTION*)


(* ::Input:: *)
(*(*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`"];


(* ::Text::Closed:: *)
(*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.";
\[Eta]Mat ::usage = "\[Eta]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.";
\[CapitalDelta]rGridMat ::usage = "\[CapitalDelta]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 \[CapitalDelta]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,\[Psi]]. Energy from pericenter and apocenter.";
J2fromrad::usage="J2fromrad[rperi,rapo,\[Psi]]. J\[Phi] from pericenter and apocenter.";
J1fromrad::usage = "J1fromrad[rperi,rapo,\[Psi]]. Jr from pericenter and apocenter.";
\[CapitalOmega]1fromrad::usage = "\[CapitalOmega]1fromrad[rperi,rapo,\[Psi]]. \[CapitalOmega]1 = \[Kappa] from pericenter and apocenter.";
\[CapitalOmega]2fromrad::usage = "\[CapitalOmega]2fromrad[rperi,rapo,\[Psi]]. \[CapitalOmega]2 = \[CapitalOmega]\[Phi] from pericenter and apocenter";
D\[CapitalOmega]1Drp::usage = "D\[CapitalOmega]1Drp[rperi,rapo,\[Psi]f]. Finite differences estimation on the grid of \[PartialD]\[CapitalOmega]1/\[PartialD]rp.";
D\[CapitalOmega]1Dra::usage = "D\[CapitalOmega]1Dra[rperi,rapo,\[Psi]f]. Finite differences estimation on the grid of \[PartialD]\[CapitalOmega]1/\[PartialD]ra.";
D\[CapitalOmega]2Drp::usage = "D\[CapitalOmega]2Drp[rperi,rapo,\[Psi]f]. Finite differences estimation on the grid of \[PartialD]\[CapitalOmega]2/\[PartialD]rp.";
D\[CapitalOmega]2Dra::usage = "D\[CapitalOmega]2Dra[rperi,rapo,\[Psi]f]. Finite differences estimation on the grid of \[PartialD]\[CapitalOmega]2/\[PartialD]ra.";
(*END USAGE.*)


(*USAGE 2D BASIS KALNAJS (1976).*)
\[Psi]basisKA::usage = "\[Psi]basisKA[l,n,r,\[Phi]]. 2D potential basis elements from Kalnajs (1976).";
UbasisKA::usage = "UbasisKA[l,n,r]. Radial component of the potential basis elements from Kalnajs (1976).";
\[Rho]basisKA::usage = "\[Rho]basisKA[l,n,r,\[Phi]]. 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).";
\[Alpha]CoeffKA::usage = "\[Alpha]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).";
\[Beta]CoeffKA::usage = "\[Beta]CoeffKA[k,l,n,j]. Coefficient in the sum for the surface density basis elements from Kalnajs (1976).";
(*.END USAGE.*)


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


(*USAGE RESPONSE MATRIX 2D.*)
RepMatCoeff2D::usage = "RepMatCoeff2D[lp,np,lq,nq,\[Omega],\[Eta],\[Psi],Ubasis,DFEL]. One coefficient of the 2D response matrix M.";
ZMatCoeff2D::usage = "ZMatCoeff2D[lp,np,nq,\[Omega],\[Eta],\[Psi],Ubasis,DFEL]. One coefficient of the bloc-diagonal 2D response matrix Z.";
agMat2D::usage = "agMat2D[lp,np,nq,m1,m2,rperi,rapo,\[Psi],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,\[Psi],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,\[Psi],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,\[Omega],rperi,rapo,\[Psi]]. Coefficient in the limited development of the denominator hMat2D to compute the 2D response matrix. (Constant contribution).";
bhMat2D::usage = "bhMat2D[m1,m2,rperi,rapo,\[Psi]]. 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,\[Psi]]. 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,\[Psi],Ubasis,DFEL]. Numerator of the integrated function to obtain the 2D response matrix.";
hMat2D::usage = "hMat2D[m1,m2,\[Omega],rperi,rapo,\[Psi]]. Denominator of the integrated function to obtain the 2D response matrix.";
JacELtorpra::usage = "JacELtorpra[rperi,rapo,\[Psi]]. Jacobian d(E,L)/d(rp,ra) of the transformation (E,L) -> (rp,ra) present in gMat.";
DEDrp::usage = "DEDrp[rperi,rapo,\[Psi]]. Partial derivative dE/drperi used in JacELtorpra.";
DEDra::usage = "DEDra[rperi,rapo,\[Psi]]. Partial derivative dE/drapo used in JacELtorpra.";
DLDrp::usage = "DLDrp[rperi,rapo,\[Psi]]. Partial derivative dL/drperi used in JacELtorpra.";
DLDra::usage = "DLDra[rperi,rapo,\[Psi]]. Partial derivative dL/drapo used in JacELtorpra.";
DFGrad::usage = "DFGrad[m1,m2,rperi,rapo,\[Psi],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[\[Omega],\[Eta],\[Psi],Ubasis,DFEL]. 2D response matrix M(\[Omega],\[Eta]).";
ZMat2D::usage = "ZMat2D[lmode,nmode,\[Omega],\[Eta],\[Psi],Ubasis,DFEL]. 2D bloc diagonal matrix ZMat2D.";
DetZMat2D::usage = "DetZMat2D[lmode,nmode,\[Omega],\[Eta],\[Psi],Ubasis,DFEL]. Det(I - Z(\[Omega],\[Eta])) used to seek for modes for 2D systems.";
(*END USAGE.*)


(*USAGE SUB-REGION INTEGRATION MATRIX.*)
\[Aleph]Mat::usage = "\[Aleph]Mat[ag,bg,cg,ah,bh,ch,\[Eta],\[CapitalDelta]r]. Analytical function to compute int_{-\[CapitalDelta]r/2}^{\[CapitalDelta]r/2} dxp int_{-\[CapitalDelta]r/2}^{\[CapitalDelta]r/2} dxa (ag + bg xp + cg xa)/(ah + bh xp + ch xa + \[Eta]), used in the sub-region matrix integration.";
\[Aleph]MatD::usage = "\[Aleph]MatD[b,c,e,f,\[Eta]]. 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.";
\[Aleph]MatAna::usage = "\[Aleph]MatAna[ag,bg,cg,ah,bh,ch,\[Eta],\[CapitalDelta]r]. Analytical function to compute int_{-\[CapitalDelta]r/2}^{\[CapitalDelta]r/2} dxp int_{-\[CapitalDelta]r/2}^{\[CapitalDelta]r/2} dxa (ag + bg xp + cg xa)/(ah + bh xp + ch xa + \[Eta]), used in the sub-region matrix integration, where it always using \[Aleph]MatDAna.";
\[Aleph]MatDAna::usage = "\[Aleph]MatDAna[b,c,e,f,\[Eta]]. Analytical integration of \[Aleph]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,\[Psi],Ubasis]. Pre-compute in parallel the needed values of WMat2D.";
DumpWMat2D::usage = "DumpWMat2D[lmode,nmode,m1bound,rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi],Ubasis]. To dump the downvalues of WMat2D previously calculated with PrepareWMat2D.";
LoadWMat2D::usage = "LoadWMat2D[lmode,nmode,m1bound,rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi],Ubasis]. To load already computed values of WMat2D.";
Prepare\[CapitalOmega]::usage = "Prepare\[CapitalOmega][\[Psi]]. Pre-compute in parallel the needed values of E,J2,J1,\[CapitalOmega]1 and \[CapitalOmega]2.";
Dump\[CapitalOmega]::usage = "Dump\[CapitalOmega][rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi]]. Dump the downvalues of Efromrad, J2fromrad,J1fromrad,\[CapitalOmega]1fromrad,\[CapitalOmega]1fromrad previously calculated with Prepare\[CapitalOmega].";
Load\[CapitalOmega]::usage = "Load\[CapitalOmega][rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi]]. Load already computed values of Efromrad, J2fromrad, J1fromrad,\[CapitalOmega]1fromrad,\[CapitalOmega]1fromrad.";
PrepareJac::usage = "PrepareJac[\[Psi]]. Pre-compute in parallel the needed values of JacELtorpra.";
DumpJac::usage = "DumpJac[rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi]]. Dump the downvalues of JacELtorpra previously calculated with Prepare\[CapitalOmega].";
LoadJac::usage = "LoadJac[rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi]]. Load already computed values of JacELtorpra.";
PrepareDFGrad::usage = "PrepareDFGrad[lmode,m1bound,\[Psi],DFEL]. Pre-compute in parallel the needed values of DFGrad for the 2D case.";
DumpDFGrad::usage = "DumpDFGrad[lmode,m1bound,rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi],DFEL]. Dump the downvalues of DFGrad previously calculated with PrepareDFGrad.";
LoadDFGrad::usage = "LoadDFGrad[lmode,m1bound,rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi],DFEL]. Load already computed values of DFGrad.";
DumpZMatCoeff2D::usage = "DumpZMatCoeff2D[lmode,nmode,m1bound,rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi],Ubasis,DFEL]. To dump the downvalues of ZMatCoeff2D already computed.";
LoadZMatCoeff2D::usage = "LoadZMatCoeff2D[lmode,nmode,m1bound,rGridmin,rGridmax,\[CapitalDelta]rGrid,\[Psi],Ubasis,DFEL]. To load the already computed values of ZMatCoeff2D.";
(*END USAGE.*)



(* ::Text::Closed:: *)
(*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 \[Eta]Mat to have well-defined integral.*)
\[Eta]Mat = 0.01;
Print["Small imaginary part used to compute a well-defined response matrix integral"];
Print["\[Eta]Mat=",\[Eta]Mat];



Print["Potentials already defined :"];
\[Psi]Kep[r_?NumericQ] := -1/r;
\[Psi]harm[r_?NumericQ] := 1/2 r^(2);
\[Psi]Hernquist[r_?NumericQ]:=- 4 \[Pi] G (1)/(2(1 + r));
\[Psi]Jaffe[r_?NumericQ] :=-  G Log[1 + 1/r];
\[Psi]iso[r_?NumericQ] :=- (G)/(1 + Sqrt[1 + r^(2)]);
\[Psi]Plummer[r_?NumericQ] := - (G)/(Sqrt[r^(2) + 1]);
\[Psi]NFW[r_?NumericQ] := - (G) Log[1+r]/(r);
Print["\[Psi]Kep, ","\[Psi]harm, ","\[Psi]Hernquist, ","\[Psi]Jaffe, ","\[Psi]iso, ","\[Psi]Plummer, ","\[Psi]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];
\[CapitalDelta]rGridMat = 0.1;
Print["\[CapitalDelta]rGridMat=",\[CapitalDelta]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 \[CapitalDelta]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,\[CapitalDelta]rGridMat},{rapotab,rperitab+2 \[CapitalDelta]rGridMat,rGridmaxMat,\[CapitalDelta]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 \[CapitalDelta]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 - \[CapitalDelta]rGridMat],choprad[ra]},{choprad[rp+\[CapitalDelta]rGridMat],choprad[ra]},{choprad[rp],choprad[ra]},{choprad[rp],choprad[ra-\[CapitalDelta]rGridMat]},{choprad[rp],choprad[ra+\[CapitalDelta]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| \[LessEqual] lbasismax, 1 \[LessEqual] n \[LessEqual] 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`"];


(* ::Text::Closed:: *)
(*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 \[CapitalDelta]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,\[CapitalDelta]rGridMat},{rapotab,rperitab+2 \[CapitalDelta]rGridMat,rGridmaxMat,\[CapitalDelta]rGridMat}],1];
nbelementsGrid = Length[listGrid];
Print["nbelementsGrid=",nbelementsGrid];
listGridlarge =
listGrid /. {rp_,ra_} -> {{choprad[rp - \[CapitalDelta]rGridMat],choprad[ra]},{choprad[rp+\[CapitalDelta]rGridMat],choprad[ra]},{choprad[rp],choprad[ra]},{choprad[rp],choprad[ra-\[CapitalDelta]rGridMat]},{choprad[rp],choprad[ra+\[CapitalDelta]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 \[CapitalDelta]rGridMat.*)
CreatelistGridlarge[] := 
Block[{},
(*In case nbelementsGrid had not been specified.*)
nbelementsGrid = Length[listGrid];
listGridlarge = listGrid /. {rp_,ra_} -> {{choprad[rp - \[CapitalDelta]rGridMat],choprad[ra]},{choprad[rp+\[CapitalDelta]rGridMat],choprad[ra]},{choprad[rp],choprad[ra]},{choprad[rp],choprad[ra-\[CapitalDelta]rGridMat]},{choprad[rp],choprad[ra+\[CapitalDelta]rGridMat]}}//Flatten[#,1]&;
listGridlarge = DeleteDuplicates[listGridlarge];
nbelementsGridlarge = Length[listGridlarge];
Print["nbelementsGridlarge=",nbelementsGridlarge];];

(*Definition of the functions i \[Rule] rp(i) and i\[Rule] 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\[Rule] 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| \[LessEqual] lbasismax, 0 \[LessEqual] n \[LessEqual] 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\[Rule]l(p),p\[Rule]n(p) .*)
lfromlistbasis2D[p_]:=listbasis2D[[p,1]];
nfromlistbasis2D[p_]:=listbasis2D[[p,2]];


(* ::Text::Closed:: *)
(*EXPLICIT ACTION MAPPING*)


(*In this section, we give the important functions allowing to obtain the important quantities describing an orbit : E, L, Jr, \[CapitalOmega]1, \[CapitalOmega]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 \[Psi]0*)
Efromrad[rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_] :=Efromrad[rperi,rapo,\[Psi]f]=(rapo^(2) \[Psi]f[rapo] - rperi^(2) \[Psi]f[rperi])/(rapo^(2) - rperi^(2));
(*SetSharedFunction[Efromrad];*)

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

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

(*Azimuthal frequency as a function of the peri/apocenters*)
(*DEPENDS ON THE POTENTIAL \[Psi]0,
We take the real part in order not to have any complex components.*)
\[CapitalOmega]2fromrad[rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_]:= \[CapitalOmega]2fromrad[rperi,rapo,\[Psi]f]=Block[{Eloc,J2loc,\[CapitalOmega]1loc},
Eloc = Efromrad[rperi,rapo,\[Psi]f];
J2loc = J2fromrad[rperi,rapo,\[Psi]f];
\[CapitalOmega]1loc = \[CapitalOmega]1fromrad[rperi,rapo,\[Psi]f];
\[CapitalOmega]1loc J2loc /(\[Pi]) NIntegrate[(1)/(rint^(2) Sqrt[2(Eloc-\[Psi]f[rint]) - J2loc^(2)/rint^(2)]),{rint,rperi,rapo},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}] // Re];
(*SetSharedFunction[\[CapitalOmega]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 \[GreaterEqual] rperi+2\[CapitalDelta]rGridMat.
THESE ESTIMATIONS SHOULD BE MADE ONLY ON THE non-largeGRID.*)
(*Finite differences estimation of \[PartialD]\[CapitalOmega]1/\[PartialD]rp.*)
D\[CapitalOmega]1Drp[rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_]:= D\[CapitalOmega]1Drp[rperi,rapo,\[Psi]f] = (\[CapitalOmega]1fromrad[choprad[rperi+\[CapitalDelta]rGridMat],rapo,\[Psi]f] - \[CapitalOmega]1fromrad[choprad[rperi-\[CapitalDelta]rGridMat],rapo,\[Psi]f])/(2 \[CapitalDelta]rGridMat);
(*Finite differences estimation of \[PartialD]\[CapitalOmega]1/\[PartialD]ra.*)
D\[CapitalOmega]1Dra[rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_]:= D\[CapitalOmega]1Dra[rperi,rapo,\[Psi]f] = (\[CapitalOmega]1fromrad[rperi,choprad[rapo+\[CapitalDelta]rGridMat],\[Psi]f] - \[CapitalOmega]1fromrad[rperi,choprad[rapo-\[CapitalDelta]rGridMat],\[Psi]f])/(2 \[CapitalDelta]rGridMat);
(*Finite differences estimation of \[PartialD]\[CapitalOmega]2/\[PartialD]rp.*)
D\[CapitalOmega]2Drp[rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_]:= D\[CapitalOmega]2Drp[rperi,rapo,\[Psi]f] = (\[CapitalOmega]2fromrad[choprad[rperi+\[CapitalDelta]rGridMat],rapo,\[Psi]f] - \[CapitalOmega]2fromrad[choprad[rperi-\[CapitalDelta]rGridMat],rapo,\[Psi]f])/(2 \[CapitalDelta]rGridMat);
(*Finite differences estimation of \[PartialD]\[CapitalOmega]2/\[PartialD]ra.*)
D\[CapitalOmega]2Dra[rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_]:= D\[CapitalOmega]2Dra[rperi,rapo,\[Psi]f] = (\[CapitalOmega]2fromrad[rperi,choprad[rapo+\[CapitalDelta]rGridMat],\[Psi]f] - \[CapitalOmega]2fromrad[rperi,choprad[rapo-\[CapitalDelta]rGridMat],\[Psi]f])/(2 \[CapitalDelta]rGridMat);



(* ::Text::Closed:: *)
(*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\[GreaterEqual]0), n radial number (n\[GreaterEqual]0). DIFFERENT THAN THE 3D CASE.
The polar coordinates are (R,\[Phi]).
To construct this basis, there are two shared parameters : RKA (scale radius) and kKA (index).*)

(*Definition of the potential elements*)
\[Psi]basisKA[l_,n_,r_?NumericQ,\[Phi]_?NumericQ] := Exp[I l \[Phi]] 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[
\[Alpha]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*)
\[Rho]basisKA[l_,n_,r_?NumericQ,\[Phi]_?NumericQ] :=Exp[I l \[Phi]] 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[
\[Beta]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.*)
\[Alpha]CoeffKA[k_,l_,n_,i_,j_] := \[Alpha]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])/(\[Pi] 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.*)
\[Beta]CoeffKA[k_,l_,n_,j_] := \[Beta]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]);




(* ::Text::Closed:: *)
(*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 \[Psi]f.*)
(*DEPENDS ON THE CHOSEN BASIS Ubasis.
The Kronecker delta has been placed such that if m\[NotEqual]m3, the expensive evaluation is not made.*)
\[Psi]basisFT2D[l_,n_,m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_, Ubasis_] :=
Which[(l!=m2),0,
True,
 WMat2D[m1,l,m2,n,rperi,rapo,\[Psi]f,Ubasis]];


(*Definition of the 2D matrix WMat*)
(*DEPENDS ON THE POTENTIAL \[Psi]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\[Equal]0 && m2\[GreaterEqual]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,\[Psi]f_, Ubasis_] :=WMat2D[m1,l,m2,n,rperi,rapo,\[Psi]f,Ubasis]=
Which[(m1 > 0) || (m1==0 && m2 >= 0),
Block[{Eloc,J2loc,\[CapitalOmega]1loc,\[CapitalOmega]2loc},
Eloc=Efromrad[rperi,rapo,\[Psi]f];
J2loc = J2fromrad[rperi,rapo,\[Psi]f];
\[CapitalOmega]1loc =\[CapitalOmega]1fromrad[rperi,rapo,\[Psi]f];
\[CapitalOmega]2loc=\[CapitalOmega]2fromrad[rperi,rapo,\[Psi]f]; 1/(\[Pi]) NIntegrate[rtod\[Theta]1dr[Eloc,J2loc,\[CapitalOmega]1loc,rint,\[Psi]f] Cos[m1 rto\[Theta]1[rperi,Eloc,J2loc,\[CapitalOmega]1loc,rint,\[Psi]f] + m2 rto\[Theta]2minus\[Psi][rperi,Eloc,J2loc,\[CapitalOmega]2loc,rint,\[Psi]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,\[Psi]f,Ubasis]];
(*SetSharedFunction[WMat];*)

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

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

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


(* ::Text::Closed:: *)
(*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 \[Psi]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_,\[Omega]_?NumericQ,\[Eta]f_?NumericQ,\[Psi]f_,Ubasis_,DFEL_] := RepMatCoeff2D[lp,np,lq,nq,\[Omega],\[Eta]f,\[Psi]f,Ubasis,DFEL]=
Which[(lp!=lq),0,
True,
ZMatCoeff2D[lp,np,nq,\[Omega],\[Eta]f,\[Psi]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 \[LessEqual] 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_,\[Omega]_?NumericQ,\[Eta]f_?NumericQ,\[Psi]f_,Ubasis_,DFEL_] := ZMatCoeff2D[lp,np,nq,\[Omega],\[Eta]f,\[Psi]f,Ubasis,DFEL]=
Which[np <= nq,
Block[{rperiloc,rapoloc},
Sum[
rperiloc = rperifromGrid[iGridsum];
rapoloc = rapofromGrid[iGridsum];
Sum[
Which[(m1sum == 0) && (m2sum == 0),0,
True,\[Aleph]Mat[
agMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,\[Psi]f,Ubasis,DFEL],
bgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,\[Psi]f,Ubasis,DFEL],
cgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,\[Psi]f,Ubasis,DFEL],
ahMat2D[m1sum,m2sum,\[Omega],rperiloc,rapoloc,\[Psi]f],
bhMat2D[m1sum,m2sum,rperiloc,rapoloc,\[Psi]f],
chMat2D[m1sum,m2sum,rperiloc,rapoloc,\[Psi]f],
\[Eta]f,
\[CapitalDelta]rGridMat]],{m1sum,-m1sumMatbound,m1sumMatbound,1},{m2sum,lp,lp,1}],{iGridsum,1,nbelementsGrid}]
],True,
ZMatCoeff2D[lp,nq,np,\[Omega],\[Eta]f,\[Psi]f,Ubasis,DFEL]];
(*SetSharedFunction[ZMatCoeff2D];*)


(*Same as ZMatCoeff2D.
ALL THE COMPUTATIONS ARE MADE USING \[Aleph]MatAna.*)
ZMatCoeffAna2D[lp_,np_,nq_,\[Omega]_?NumericQ,\[Eta]f_?NumericQ,\[Psi]f_,Ubasis_,DFEL_] := ZMatCoeffAna2D[lp,np,nq,\[Omega],\[Eta]f,\[Psi]f,Ubasis,DFEL]=
Which[np <= nq,
Block[{rperiloc,rapoloc},
Sum[
rperiloc = rperifromGrid[iGridsum];
rapoloc = rapofromGrid[iGridsum];
Sum[
Which[(m1sum == 0) && (m2sum == 0),0,
True,\[Aleph]MatAna[
agMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,\[Psi]f,Ubasis,DFEL],
bgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,\[Psi]f,Ubasis,DFEL],
cgMat2D[lp,np,nq,m1sum,m2sum,rperiloc,rapoloc,\[Psi]f,Ubasis,DFEL],
ahMat2D[m1sum,m2sum,\[Omega],rperiloc,rapoloc,\[Psi]f],
bhMat2D[m1sum,m2sum,rperiloc,rapoloc,\[Psi]f],
chMat2D[m1sum,m2sum,rperiloc,rapoloc,\[Psi]f],
\[Eta]f,
\[CapitalDelta]rGridMat]],{m1sum,-m1sumMatbound,m1sumMatbound,1},{m2sum,lp,lp,1}],{iGridsum,1,nbelementsGrid}]
],True,
ZMatCoeffAna2D[lp,nq,np,\[Omega],\[Eta]f,\[Psi]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+\[CapitalDelta]r)-f(r0-\[CapitalDelta]r))/(2 \[CapitalDelta]r).
We note that the only function which is not memoized is ahMat2D, because it has an \[Omega]-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,\[Psi]f_,Ubasis_,DFEL_]:=agMat2D[lp,np,nq,m1,m2,rperi,rapo,\[Psi]f,Ubasis,DFEL]=gMat2D[lp,np,nq,m1,m2,rperi,rapo,\[Psi]f,Ubasis,DFEL];

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

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

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

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

(*This function is memoized and does not depend on \[Omega].
It requires to evaluate hMat2D in \[Omega]=0.*)
chMat2D[m1_,m2_,rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_] := 
chMat2D[m1,m2,rperi,rapo,\[Psi]f]=(hMat2D[m1,m2,0.,rperi,choprad[rapo+\[CapitalDelta]rGridMat],\[Psi]f] - hMat2D[m1,m2,0.,rperi,choprad[rapo-\[CapitalDelta]rGridMat],\[Psi]f])/(2 \[CapitalDelta]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,\[Psi]f_,Ubasis_,DFEL_] :=
gMat2D[lp,np,nq,m1,m2,rperi,rapo,\[Psi]f,Ubasis,DFEL]=(2 \[Pi])^(2) JacELtorpra[rperi,rapo,\[Psi]f] (1/\[CapitalOmega]1fromrad[rperi,rapo,\[Psi]f]) DFGrad[m1,m2,rperi,rapo,\[Psi]f,DFEL] WMat2D[m1,lp,m2,np,rperi,rapo,\[Psi]f,Ubasis] WMat2D[m1,lp,m2,nq,rperi,rapo,\[Psi]f,Ubasis]  ;
(*SetSharedFunction[gMat2D];*)

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

(*Jacobian of the transformation (E,L)\[Rule](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,\[Psi]f_] := JacELtorpra[rperi,rapo,\[Psi]f]=Abs[DEDrp[rperi,rapo,\[Psi]f] DLDra[rperi,rapo,\[Psi]f] - DLDrp[rperi,rapo,\[Psi]f] DEDra[rperi,rapo,\[Psi]f]];
(*SetSharedFunction[JacELtorpra];*)

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

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

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

(*Calculation of \partial L/ \partial ra.*)
DLDra[rperi_?NumericQ,rapo_?NumericQ,\[Psi]f_] := Block[{D\[Psi]Draloc},
D\[Psi]Draloc = D[\[Psi]f[rder],{rder,1}] /. rder -> rapo;
(rapo rperi^(2) (2 rperi^(2) (\[Psi]f[rperi] - \[Psi]f[rapo]) + D\[Psi]Draloc (rapo^(3) - rapo rperi^(2))))/(Sqrt[2] (rapo - rperi)^(2) (rapo+rperi)^(2) Sqrt[(rapo^(2) rperi^(2) (\[Psi]f[rapo] - \[Psi]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,\[Psi]f_,DFEL_] := DFGrad[m1,m2,rperi,rapo,\[Psi]f,DFEL]=Block[{Eloc,Lloc,\[CapitalOmega]1loc,\[CapitalOmega]2loc,DFDEloc},
Eloc = Efromrad[rperi,rapo,\[Psi]f];
Lloc = J2fromrad[rperi,rapo,\[Psi]f];
\[CapitalOmega]1loc= \[CapitalOmega]1fromrad[rperi,rapo,\[Psi]f];
\[CapitalOmega]2loc = \[CapitalOmega]2fromrad[rperi,rapo,\[Psi]f];
DFDEloc = DFDE[Eloc,Lloc,DFEL];
m1 \[CapitalOmega]1loc DFDEloc + m2 (\[CapitalOmega]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 \[GreaterEqual] 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[\[Omega]_?NumericQ,\[Eta]f_?NumericQ,\[Psi]f_,Ubasis_,DFEL_] := Block[{halfMloc,Mloc},
halfMloc=Table[RepMatCoeff2D[lfromlistbasis2D[ptab],nfromlistbasis2D[ptab],lfromlistbasis2D[qtab],nfromlistbasis2D[qtab],\[Omega],\[Eta]f,\[Psi]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.
+ \[Omega] : the real part frequency of evaluation.
+ \[Eta] : 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 \[LessEqual] np \[LessEqual] nmode.
*)
ZMat2D[lmodef_,nmodef_,\[Omega]f_?NumericQ,\[Eta]f_?NumericQ,\[Psi]f_,Ubasis_,DFEL_] := Table[ZMatCoeff2D[lmodef,nptab,nqtab,\[Omega]f,\[Eta]f,\[Psi]f,Ubasis,DFEL],{nptab,0,nmodef,1},{nqtab,0,nmodef,1}];

(*Same as ZMat2D.
ALL THE COMPUTATIONS ARE MADE USING \[Aleph]MatAna.
*)
ZMatAna2D[lmodef_,nmodef_,\[Omega]f_?NumericQ,\[Eta]f_?NumericQ,\[Psi]f_,Ubasis_,DFEL_] := Table[ZMatCoeffAna2D[lmodef,nptab,nqtab,\[Omega]f,\[Eta]f,\[Psi]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.
+ \[Omega] : the real part of the frequency of evaluation.
+ \[Eta] : the imaginary part of the frequency of evaluation.*)
DetZMat2D[lmodef_,nmodef_,\[Omega]f_?NumericQ,\[Eta]f_?NumericQ,\[Psi]f_,Ubasis_,DFEL_]:=
Det[IdentityMatrix[nmodef+1] -  ZMat2D[lmodef,nmodef,\[Omega]f,\[Eta]f,\[Psi]f,Ubasis,DFEL]];



(* ::Text::Closed:: *)
(*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.*)
\[Aleph]MatD = Compile[{{bf,_Real},{cf,_Real},{ef,_Real},{ff,_Real},{\[Eta]f,_Real}},
1/(16 ef^2 ff^2) (-I (ff (bf (ef^2-(-2+ff-2 I \[Eta]f)^2)+4 ef (2+ef+2 I \[Eta]f))-cf ef (2+ef+2 I \[Eta]f)^2) (\[Pi]-2 ArcTan[(2+ef-ff)/(2 \[Eta]f)])-I (-ff (bf (-ef^2+(2+ff+2 I \[Eta]f)^2)+4 ef (-2+ef-2 I \[Eta]f))-cf ef (-2+ef-2 I \[Eta]f)^2) (\[Pi]-2 ArcTan[(2-ef+ff)/(2 \[Eta]f)])+I (ff (bf (ef^2-(-2+ff-2 I \[Eta]f)^2)-4 ef (-2+ef-2 I \[Eta]f))-cf ef (-2+ef-2 I \[Eta]f)^2) (\[Pi]+2 ArcTan[(-2+ef+ff)/(2 \[Eta]f)])+I (ff (bf (ef^2-(2+ff+2 I \[Eta]f)^2)+4 ef (2+ef+2 I \[Eta]f))-cf ef (2+ef+2 I \[Eta]f)^2) (\[Pi]-2 ArcTan[(2+ef+ff)/(2 \[Eta]f)])+ff (ff (8 ef+bf (-4+2 ef+ff-4 I \[Eta]f))+cf ef (-4+2 ef-ff-4 I \[Eta]f)+2 (-4+cf) ef ff Log[1/2 (2-ef-ff+2 I \[Eta]f)])-ff (ff (8 ef+bf (-4-2 ef+ff-4 I \[Eta]f))-cf ef (4+2 ef+ff+4 I \[Eta]f)+2 (-4+cf) ef ff Log[1/2 (2+ef-ff+2 I \[Eta]f)])-ff (ff (-8 ef+bf (4-2 ef+ff+4 I \[Eta]f))-cf ef (-4+2 ef+ff-4 I \[Eta]f)+2 (4+cf) ef ff Log[1/2 (2-ef+ff+2 I \[Eta]f)])+ff (ff (-8 ef+bf (4+2 ef+ff+4 I \[Eta]f))+cf ef (4+2 ef-ff+4 I \[Eta]f)+2 (4+cf) ef ff Log[1/2 (2+ef+ff+2 I \[Eta]f)])+(ff (bf (-ef^2+(-2+ff-2 I \[Eta]f)^2)-4 ef (2+ef+2 I \[Eta]f))+cf ef (2+ef+2 I \[Eta]f)^2) Log[1/4 (4+ef^2-2 ef (-2+ff)-4 ff+ff^2+4 \[Eta]f^2)]-(ff (bf (-ef^2+(-2+ff-2 I \[Eta]f)^2)+4 ef (-2+ef-2 I \[Eta]f))+cf ef (-2+ef-2 I \[Eta]f)^2) Log[1/4 (4+ef^2+2 ef (-2+ff)-4 ff+ff^2+4 \[Eta]f^2)]+(ff (bf (-ef^2+(2+ff+2 I \[Eta]f)^2)+4 ef (-2+ef-2 I \[Eta]f))+cf ef (-2+ef-2 I \[Eta]f)^2) Log[1/4 (4+ef^2+4 ff+ff^2-2 ef (2+ff)+4 \[Eta]f^2)]-(ff (bf (-ef^2+(2+ff+2 I \[Eta]f)^2)-4 ef (2+ef+2 I \[Eta]f))+cf ef (2+ef+2 I \[Eta]f)^2) Log[1/4 (4+ef^2+4 ff+ff^2+2 ef (2+ff)+4 \[Eta]f^2)]),CompilationTarget->"C",CompilationOptions->{"ExpressionOptimization" -> True,"InlineCompiledFunctions"->True,"InlineExternalDefinitions"->True},RuntimeOptions->"Speed"];


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


(*Analytical version of \[Aleph]Mat.*)
\[Aleph]MatAna=Compile[{{agf,_Real},{bgf,_Real},{cgf,_Real},{ahf,_Real},{bhf,_Real},{chf,_Real},{\[Eta]f,_Real},{\[CapitalDelta]f,_Real}},
Block[{bloc,cloc,eloc,floc,\[Eta]loc},
bloc = (bgf \[CapitalDelta]f)/(agf);
cloc = (cgf \[CapitalDelta]f)/(agf);
eloc = (bhf \[CapitalDelta]f)/(ahf);
floc = (chf \[CapitalDelta]f)/(ahf);
\[Eta]loc = (\[Eta]f)/(ahf);
(agf)/(ahf) (\[CapitalDelta]f)^(2) \[Aleph]MatDAna[bloc,cloc,eloc,floc,\[Eta]loc]],
{{\[Aleph]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.*)
\[Aleph]MatDAna[bf_?NumericQ,cf_?NumericQ,ef_?NumericQ,ff_?NumericQ,\[Eta]f_?NumericQ] := NIntegrate[(1 + bf xint + cf yint)/(1 + ef xint + ff yint + I \[Eta]f),{xint,-0.5,0.5},{yint,-0.5,0.5},PrecisionGoal->4,AccuracyGoal->4,Method->{Automatic,SymbolicProcessing->0}];


(* ::Text::Closed:: *)
(*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 \[LessEqual] np \[LessEqual] 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,\[CapitalOmega]1,\[CapitalOmega]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,\[Psi],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 \[Rule] We had to set $DistributedContexts = Automatic.*)
PrepareWMat2D[lmodef_,nmodef_,m1boundf_,\[Psi]f_,Ubasisf_]:=Block[{tabValWMat2D,vloc},
tabValWMat2D = ParallelTable[
Table[
If[(m1tab==0)&&(m2tab==0),Null,
{m1tab,ltab,m2tab,nptab,rperifromGridlarge[iGridtab],rapofromGridlarge[iGridtab],\[Psi]f,Ubasisf,WMat2D[m1tab,ltab,m2tab,nptab,rperifromGridlarge[iGridtab],rapofromGridlarge[iGridtab],\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]f_,Ubasisf_] := Block[{namefileloc},
namefileloc = adressfiles<>"DumpWMat2D_lmode_"<>ToString[lmodef]<>"_nmode_"<>ToString[nmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]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,\[CapitalOmega]1,\[CapitalOmega]2.
The only argument needed is the background potential \[Psi]f.*)
Prepare\[CapitalOmega][\[Psi]f_]:=Block[{tab\[CapitalOmega],rploc,raloc,\[Psi]loc,vloc},
tab\[CapitalOmega] = ParallelTable[
rploc = rperifromGridlarge[iGridtab];
raloc =rapofromGridlarge[iGridtab];
{rploc,raloc,\[Psi]f,
Efromrad[rploc,raloc,\[Psi]f],
J2fromrad[rploc,raloc,\[Psi]f],
J1fromrad[rploc,raloc,\[Psi]f],
\[CapitalOmega]1fromrad[rploc,raloc,\[Psi]f],
\[CapitalOmega]2fromrad[rploc,raloc,\[Psi]f]},{iGridtab,1,nbelementsGridlarge}];

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

(*Function which dumps the downvalues of E,J2,J1,\[CapitalOmega]1,\[CapitalOmega]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\[CapitalOmega][rGridminf_,rGridmaxf_,\[CapitalDelta]rGridf_,\[Psi]f_] := Block[{namefileloc},
(*For E.*)
namefileloc = adressfiles<>"DumpE_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
Put[DownValues[Efromrad],namefileloc];
(*For J2.*)
namefileloc = adressfiles<>"DumpJ2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
Put[DownValues[J2fromrad],namefileloc];
(*For J1.*)
namefileloc = adressfiles<>"DumpJ1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
Put[DownValues[J1fromrad],namefileloc];
(*For \[CapitalOmega]1.*)
namefileloc = adressfiles<>"DumpOmega1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
Put[DownValues[\[CapitalOmega]1fromrad],namefileloc];
(*For \[CapitalOmega]2.*)
namefileloc = adressfiles<>"DumpOmega2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
Put[DownValues[\[CapitalOmega]2fromrad],namefileloc];
(*DownValues[\[CapitalOmega]2fromrad] >> namefileloc;*)
];

(*This function allows to load previously calculated values of E,J2,J1,\[CapitalOmega]1,\[CapitalOmega]2.
It depends on the same parameters a Dump\[CapitalOmega].
ATTENTION : In the current implementation, once loaded, it overwrites all previously memoized values of Efromrad,J2fromrad,J1fromrad,\[CapitalOmega]1fromrad,\[CapitalOmega]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\[CapitalOmega][rGridminf_,rGridmaxf_,\[CapitalDelta]rGridf_,\[Psi]f_]:= Block[{namefileloc},
(*For E.*)
namefileloc = adressfiles<>"DumpE_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
DownValues[Efromrad] = ReadList[namefileloc][[1]];
(*For J2.*)
namefileloc = adressfiles<>"DumpJ2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
DownValues[J2fromrad] = ReadList[namefileloc][[1]];
(*For J1.*)
namefileloc = adressfiles<>"DumpJ1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
DownValues[J1fromrad] = ReadList[namefileloc][[1]];
(*For \[CapitalOmega]1.*)
namefileloc = adressfiles<>"DumpOmega1_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
DownValues[\[CapitalOmega]1fromrad] = ReadList[namefileloc][[1]];
(*For \[CapitalOmega]2.*)
namefileloc = adressfiles<>"DumpOmega2_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>".txt";
DownValues[\[CapitalOmega]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 \[Psi]f.*)
PrepareJac[\[Psi]f_]:=Block[{tabJac,rploc,raloc,\[Psi]loc,vloc},
tabJac = ParallelTable[
rploc = rperifromGridlarge[iGridtab];
raloc =rapofromGridlarge[iGridtab];
{rploc,raloc,\[Psi]f,
JacELtorpra[rploc,raloc,\[Psi]f]},{iGridtab,1,nbelementsGridlarge}];

Do[vloc = tabJac[[ind]];
rploc = vloc[[1]];
raloc = vloc[[2]];
\[Psi]loc = vloc[[3]];
JacELtorpra[rploc,raloc,\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]f_] := Block[{namefileloc},
namefileloc = adressfiles<>"DumpJac_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]f_]:= Block[{namefileloc,listread},
namefileloc = adressfiles<>"DumpJac_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]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_,\[Psi]f_,DFELf_]:=Block[{tabDFGrad,rperiloc,rapoloc,vloc},
tabDFGrad = ParallelTable[
rperiloc = rperifromGridlarge[iGridtab];
rapoloc = rapofromGridlarge[iGridtab];
Table[
{m1tab,m2tab,rperiloc,rapoloc,\[Psi]f,DFELf,DFGrad[m1tab,m2tab,rperiloc,rapoloc,\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]f_,DFELf_] := Block[{namefileloc},
namefileloc = adressfiles<>"DumpDFGrad_lmode_"<>ToString[lmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]f_,DFELf_]:= Block[{namefileloc,listread},
namefileloc = adressfiles<>"DumpDFGrad_lmode_"<>ToString[lmodef]<>"_m1bound_"<>ToString[m1boundf]<>"_rmin_"<>ToString[rGridminf]<>"_rmax_"<>ToString[rGridmaxf]<>"_Deltar_"<>ToString[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]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[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]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[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]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_,\[CapitalDelta]rGridf_,\[Psi]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[\[CapitalDelta]rGridf]<>"_pot_"<>StringDrop[ToString[FullForm[\[Psi]f]],1]<>"_basis_"<>ToString[FullForm[Ubasisf]]<>"_DFEL_"<>ToString[FullForm[DFELf]]<>".txt";
listread = ReadList[namefileloc][[1]];
DownValues[ZMatCoeff2D] = listread;
]


End[];
EndPackage[];
