update for OrbitalElement v2dev

examples/IsochroneE/runExampleIsochrone.jl

@@ -20,18 +20,17 @@ basis = AstroBasis.CB73Basis(lmax=lmax, nradial=nradial,G=G,rb=rb)
 
 
 # Model Potential
-const modelname = "IsochroneE"
+const modelname = "IsochroneE2"
 const bc, M = 1.,1. # G is defined above: must agree with basis!
-model = OrbitalElements.IsochronePotential()
-
+model = OrbitalElements.NumericalIsochrone()
 
 rmin = 0.0
 rmax = Inf
 
 
 dfname = "roi1.0"
 
-function ndFdJ(n1::Int64,n2::Int64,E::Float64,L::Float64,ndotOmega::Float64;bc::Float64=1.,M::Float64=1.,astronomicalG::Float64=1.,Ra::Float64=1.)
+function ndFdJ(n1::Int64,n2::Int64,E::Float64,L::Float64,ndotOmega::Float64;bc::Float64=1.,M::Float64=1.,astronomicalG::Float64=1.,Ra::Float64=1.0)
 
     Q = OrbitalElements.isochroneQROI(E,L,Ra,bc,M,astronomicalG)
 
@@ -79,18 +78,17 @@ Etamax   = 0.1
 
 
 VERBOSE   = 2
-OVERWRITE = false
+OVERWRITE = true
 VMAPN     = 1
 ADAPTIVEKW= false
 
-OEparams = OrbitalElements.OrbitalParameters(Ω₀=OrbitalElements.Ω₀(model),
-                                             EDGE=OrbitalElements.DEFAULT_EDGE,TOLECC=OrbitalElements.DEFAULT_TOLECC,TOLA=OrbitalElements.DEFAULT_TOLA,
+OEparams = OrbitalElements.OrbitalParameters(EDGE=OrbitalElements.DEFAULT_EDGE,TOLECC=OrbitalElements.DEFAULT_TOLECC,TOLA=OrbitalElements.DEFAULT_TOLA,
                                              NINT=OrbitalElements.DEFAULT_NINT,
                                              da=OrbitalElements.DEFAULT_DA,de=OrbitalElements.DEFAULT_DE,
                                              ITERMAX=OrbitalElements.DEFAULT_ITERMAX,invε=OrbitalElements.DEFAULT_TOL)
 
 
-Parameters = LinearResponse.LinearParameters(basis,Orbitalparams=OEparams,Ku=Ku,Kv=Kv,Kw=Kw,
+Parameters = LinearResponse.LinearParameters(basis,Orbitalparams=OEparams,Ω₀=OrbitalElements.frequency_scale(model),Ku=Ku,Kv=Kv,Kw=Kw,
                                              modelname=modelname,dfname=dfname,
                                              wmatdir=wmatdir,gfuncdir=gfuncdir,modedir=modedir,axidir=modedir,
                                              lharmonic=lharmonic,n1max=n1max,
@@ -130,7 +128,7 @@ epsilon = abs.(reshape(tabdet,nEta,nOmega))
 contour(tabOmega,tabEta,log10.(epsilon), levels=10, color=:black, #levels=[-2.0, -1.5, -1.0, -0.5, -0.25, 0.0], 
         xlabel="Re[ω]", ylabel="Im[ω]", xlims=(Omegamin,Omegamax), ylims=(Etamin,Etamax),
         clims=(-2, 0), aspect_ratio=:equal, legend=false)
-savefig("ROIdeterminant.png")
+savefig("ROIdeterminant2.png")
 
 
 # find a pole by using gradient descent
@@ -147,5 +145,5 @@ if isfinite(bestomg)
         nmode = 100
         ModeRadius,ModePotentialShape,ModeDensityShape = LinearResponse.GetModeShape(basis,modeRmin,modeRmax,nmode,EM,Parameters)
 else
-        println("runExamplePlummer.jl: no mode found.")
+        println("runExampleIsochrone.jl: no mode found.")
 end

examples/PlummerE/runExamplePlummer.jl

@@ -24,12 +24,12 @@ const bc, M = 1.,1. # G is defined above: must agree with basis!
 model = OrbitalElements.PlummerPotential()
 
 # Model Distribution Function
-dfname = "roi0.95"
+dfname = "roi0.75"
 
 function ndFdJ(n1::Int64,n2::Int64,
                E::Float64,L::Float64,
                ndotOmega::Float64;
-               bc::Float64=1.,M::Float64=1.,astronomicalG::Float64=1.,Ra::Float64=0.95)
+               bc::Float64=1.,M::Float64=1.,astronomicalG::Float64=1.,Ra::Float64=0.75)
 
     return OrbitalElements.plummer_ROI_ndFdJ(n1,n2,E,L,ndotOmega,bc,M,astronomicalG,Ra)
 
@@ -69,14 +69,17 @@ OVERWRITE = false
 VMAPN     = 1
 ADAPTIVEKW= false
 
-OEparams = OrbitalElements.OrbitalParameters(Ω₀=OrbitalElements.Ω₀(model),
+RMIN = 0.0
+RMAX = Inf
+
+OEparams = OrbitalElements.OrbitalParameters(rmin=RMIN,rmax=RMAX,
                                              EDGE=OrbitalElements.DEFAULT_EDGE,TOLECC=OrbitalElements.DEFAULT_TOLECC,TOLA=OrbitalElements.DEFAULT_TOLA,
                                              NINT=OrbitalElements.DEFAULT_NINT,
                                              da=OrbitalElements.DEFAULT_DA,de=OrbitalElements.DEFAULT_DE,
                                              ITERMAX=OrbitalElements.DEFAULT_ITERMAX,invε=OrbitalElements.DEFAULT_TOL)
 
 
-Parameters = LinearResponse.LinearParameters(basis,Orbitalparams=OEparams,Ku=Ku,Kv=Kv,Kw=Kw,
+Parameters = LinearResponse.LinearParameters(basis,Orbitalparams=OEparams,Ω₀=frequency_scale(model),Ku=Ku,Kv=Kv,Kw=Kw,
                                              modelname=modelname,dfname=dfname,
                                              wmatdir=wmatdir,gfuncdir=gfuncdir,modedir=modedir,axidir=modedir,
                                              lharmonic=lharmonic,n1max=n1max,
@@ -101,6 +104,7 @@ Parameters = LinearResponse.LinearParameters(basis,Orbitalparams=OEparams,Ku=Ku,
 
 MMat, tabaMcoef, tabωminωmax = LinearResponse.PrepareM(Parameters)
 
+
 # construct a grid of frequencies to probe
 tabω = LinearResponse.gridomega(Omegamin,Omegamax,nOmega,Etamin,Etamax,nEta)
 @time tabRMreal, tabRMimag = LinearResponse.RunMatrices(tabω,FHT,Parameters)

src/GFunc.jl

@@ -12,7 +12,7 @@ function to compute G(u)
 """
 function MakeGu(ndFdJ::Function,
                 n1::Int64,n2::Int64,
-                Wdata::WMatdataType,
+                Wdata::FourierTransformedBasisData,
                 tabu::Vector{Float64},
                 params::LinearParameters)
 
@@ -44,7 +44,7 @@ function MakeGu(ndFdJ::Function,
     δvp = 1.0/params.Kv
 
     # remove dimensionality from Ω mapping
-    dimensionl = (1.0/params.Orbitalparams.Ω₀)
+    dimensionl = (1.0/Wdata.Ω₀)
 
     # Integration step volume
     δvol = δvp * dimensionl
@@ -80,7 +80,8 @@ function MakeGu(ndFdJ::Function,
 
             # (u,v) -> (α,β).
             # Renormalized. (2/(ωmax-ωmin) * |∂(α,β)/∂(u,v)|)
-            RenormalizedJacαβ = OrbitalElements.RenormalizedJacUVToαβ(n1,n2,uval,vval)
+            resonance = OrbitalElements.Resonance(n1,n2,Wdata.ωmin,Wdata.ωmax)
+            RenormalizedJacαβ = (2/(Wdata.ωmax-Wdata.ωmin))*OrbitalElements.uv_to_αβ_jacobian(uval,vval,resonance)
 
             # (α,β) -> (E,L): this is the most expensive function here,
             # so we have pre-tabulated it
@@ -163,7 +164,7 @@ function RunGfunc(ndFdJ::Function,
         # Check if enough basis element in this file (throw error if not)
         (params.nradial <= read(file,"LinearParameters/nradial")) || error("Not enough basis element in WMat file for ($n1,$n2) resonance.")
         # Construct the Wdata structure for the file
-        Wdata      = WMatdataType(read(file,"omgmin"),read(file,"omgmax"),read(file,"tabvminmax"),                  # Mapping parameters
+        Wdata      = FourierTransformedBasisData(read(file,"omgmin"),read(file,"omgmax"),read(file,"Ω₀"),read(file,"tabvminmax"),                  # Mapping parameters
                                   read(file,"wmat"),                                                                # Basis FT
                                   read(file,"UVmat"),read(file,"Omgmat"),read(file,"AEmat"),read(file,"ELmat"),     # Mappings
                                   read(file,"jELABmat"))                                                            # Jacobians

src/Isochrone/WMatIsochrone.jl

@@ -217,7 +217,7 @@ function WBasisFTIsochrone(a::Float64,e::Float64,
                   Ω1::Float64,Ω2::Float64,
                   n1::Int64,n2::Int64,
                   bc::Float64,M::Float64,G::Float64,
-                  basisFT::BasisFTtype,
+                  basisFT::FourierTransformedBasis,
                   params::LinearParameters)
 
     # Basis FT
@@ -230,7 +230,7 @@ without Ω1, Ω2
 function WBasisFTIsochrone(a::Float64,e::Float64,
                   n1::Int64,n2::Int64,
                   bc::Float64,M::Float64,G::Float64,
-                  basisFT::BasisFTtype,
+                  basisFT::FourierTransformedBasis,
                   params::LinearParameters)
 
     # Frequencies
@@ -242,15 +242,15 @@ end
 
 
 """
-    MakeWmatUVIsochrone(n1::Int64, n2::Int64, tabu::Vector{Float64}, basisFT::BasisFTtype, bc::Float64, M::Float64, G::Float64, params::LinearParameters)
+    MakeWmatUVIsochrone(n1::Int64, n2::Int64, tabu::Vector{Float64}, basisFT::FourierTransformedBasis, bc::Float64, M::Float64, G::Float64, params::LinearParameters)
 
 Construct the matrix `Wdata` based on the given parameters for an isochrone potential.
 
 # Arguments
 - `n1::Int64`: Integer representing some parameter.
 - `n2::Int64`: Integer representing some parameter.
 - `tabu::Vector{Float64}`: Vector of `Float64` representing values of `u`.
-- `basisFT::BasisFTtype`: Type representing some basis Fourier transform.
+- `basisFT::FourierTransformedBasis`: Type representing some basis Fourier transform.
 - `bc::Float64`: Float representing some parameter.
 - `M::Float64`: Float representing mass.
 - `G::Float64`: Float representing gravitational constant.
@@ -265,7 +265,7 @@ Construct the matrix `Wdata` based on the given parameters for an isochrone pote
 """
 function MakeWmatUVIsochrone(n1::Int64,n2::Int64,
                              tabu::Vector{Float64},
-                             basisFT::BasisFTtype,
+                             basisFT::FourierTransformedBasis,
                              bc::Float64,M::Float64,G::Float64,
                              params::LinearParameters)
 

src/ModeComputation/FindPole.jl

@@ -36,7 +36,7 @@ function GoStep(omgval::ComplexF64,
     # curpos = [real(detXival);imag(detXival)]
     # increment = jacobian \ (-curpos)
 
-    increment1, increment2 = OrbitalElements.inverse2Dlinear(dXirir,dXiupr,dXirii,dXiupi,-real(detXival),-imag(detXival))
+    increment1, increment2 = OrbitalElements._inverse_2d(dXirir,dXiupr,dXirii,dXiupi,-real(detXival),-imag(detXival))
 
     # omgval += increment[1] + increment[2]*1im
     omgval += increment1 + increment2*1im

src/ResponseMatrix.jl

@@ -48,7 +48,7 @@ function tabM!(ω::ComplexF64,
     fill!(tabM,0.0 + 0.0*im)
 
     # Rescale to get dimensionless frequency
-    ωnodim = ω/params.Orbitalparams.Ω₀
+    ωnodim = ω/params.Ω₀
 
     # loop over the resonances: no threading here because we parallelise over frequencies
     for nres = 1:nbResVec
@@ -60,7 +60,7 @@ function tabM!(ω::ComplexF64,
         ωmin, ωmax = tabωminωmax[1,nres], tabωminωmax[2,nres]
 
         # get the rescaled frequency
-        ϖ = OrbitalElements.Getϖ(ωnodim,ωmin,ωmax)
+        ϖ = OrbitalElements.rescaled_ϖ(ωnodim,ωmin,ωmax)
 
         # get the integration values
         FiniteHilbertTransform.GettabD!(ϖ,FHT)

src/Utils/ParameterStructure.jl

@@ -9,6 +9,9 @@ struct LinearParameters
     # Orbital Elements parameters
     Orbitalparams::OrbitalElements.OrbitalParameters
 
+    # frequency scale
+    Ω₀::Float64
+
     # Basis parameters
     # dimension and nradial are repeated (for allocation issues - intensively called)
     dimension::Int64
@@ -51,7 +54,7 @@ end
     LinearParameters(basis;Orbitalparams,Ku,Kv,Kw,VMAPN,ADAPTIVEKW,KuTruncation,modelname,dfname,wmatdir,gfuncdir,modedir,OVERWRITE,lharmonic,n1max,VERBOSE)
 """
 function LinearParameters(basis::AstroBasis.AbstractAstroBasis;
-                          Orbitalparams::OrbitalElements.OrbitalParameters=OrbitalElements.OrbitalParameters(),
+                          Orbitalparams::OrbitalElements.OrbitalParameters=OrbitalElements.OrbitalParameters(),Ω₀::Float64=1.0,
                           Ku::Int64=200,Kv::Int64=200,Kw::Int64=200,
                           VMAPN::Int64=1,ADAPTIVEKW::Bool=false,KuTruncation::Int64=10000,
                           modelname::String="model",dfname::String="df",
@@ -69,7 +72,7 @@ function LinearParameters(basis::AstroBasis.AbstractAstroBasis;
     # Resonance vectors
     nbResVec, tabResVec = MakeTabResVec(lharmonic,n1max,dimension)
 
-    return LinearParameters(Orbitalparams,dimension,nradial,Basisparams,
+    return LinearParameters(Orbitalparams,Ω₀,dimension,nradial,Basisparams,
                             Ku,Kv,Kw,VMAPN,ADAPTIVEKW,KuTruncation,
                             modelname,dfname,
                             wmatdir,gfuncdir,axidir,modedir,OVERWRITE,

src/WMat.jl

@@ -8,23 +8,24 @@
 """
 structure to store Fourier Transform values of basis elements
 """
-struct BasisFTtype{BT<:AstroBasis.AbstractAstroBasis}
+struct FourierTransformedBasis{BT<:AstroBasis.AbstractAstroBasis}
     basis::BT
     UFT::Array{Float64}
 end
 
 function BasisFTcreate(basis::BT) where {BT<:AstroBasis.AbstractAstroBasis}
 
-    return BasisFTtype(basis,zeros(Float64,basis.nradial))
+    return FourierTransformedBasis(basis,zeros(Float64,basis.nradial))
 end
 
 
 """
 structure to store the W matrix computation results
 """
-struct WMatdataType
+struct FourierTransformedBasisData
     ωmin::Float64
     ωmax::Float64
+    Ω₀::Float64
     tabvminmax::Array{Float64,2}
 
     tabW::Array{Float64,3}
@@ -43,13 +44,13 @@ function WMatdataCreate(model::OrbitalElements.Potential,
                         params::LinearParameters)
 
     # compute the frequency scaling factors for this resonance
-    ωmin, ωmax = OrbitalElements.Findωminωmax(n1,n2,model,params.Orbitalparams)
+    ωmin, ωmax = OrbitalElements.frequency_extrema(n1,n2,model,params.Orbitalparams)
 
     # Useful parameters
     nradial = params.nradial
     Ku, Kv = params.Ku, params.Kv
 
-    return WMatdataType(ωmin,ωmax,zeros(Float64,2,Ku),
+    return FourierTransformedBasisData(ωmin,ωmax,OrbitalElements.frequency_scale(model),zeros(Float64,2,Ku),
                         zeros(Float64,nradial,Kv,Ku),
                         zeros(Float64,2,Kv,Ku),zeros(Float64,2,Kv,Ku),zeros(Float64,2,Kv,Ku),zeros(Float64,2,Kv,Ku), # Orbital mappings
                         zeros(Float64,Kv,Ku))
@@ -104,10 +105,10 @@ function Wintegrand(model::OrbitalElements.Potential,
 
 
     # Current location of the radius, r=r(w)
-    rval = OrbitalElements.ru(w,a,e)
+    rval = OrbitalElements.radius_from_anomaly(w,a,e)
 
     # Current value of the radial frequency integrand (almost dθ/dw)
-    gval = OrbitalElements.ΘAE(model,w,a,e,params.Orbitalparams)
+    gval = OrbitalElements.Θ(w,a,e,model,params.Orbitalparams)
 
 
     # collect the basis elements (in place!)
@@ -141,7 +142,7 @@ function WBasisFT(model::OrbitalElements.Potential,
     dw = -(2.0)/(Kwp)
 
     # need angular momentum
-    Lval = OrbitalElements.LFromAE(model,a,e,params.Orbitalparams)
+    _,Lval = OrbitalElements.EL_from_ae(a,e,model,params.Orbitalparams)
 
     # Initialise the state vectors: w, θ1, (θ2-psi)
     # Reverse integration, starting at apocenter
@@ -259,7 +260,7 @@ function WBasisFT(model::OrbitalElements.Potential,
                   a::Float64,e::Float64,
                   Ω1::Float64,Ω2::Float64,
                   n1::Int64,n2::Int64,
-                  basisFT::BasisFTtype,
+                  basisFT::FourierTransformedBasis,
                   params::LinearParameters)
 
     # Basis FT
@@ -272,11 +273,11 @@ without Ω1, Ω2
 function WBasisFT(model::OrbitalElements.Potential,
                   a::Float64,e::Float64,
                   n1::Int64,n2::Int64,
-                  basisFT::BasisFTtype,
+                  basisFT::FourierTransformedBasis,
                   params::LinearParameters)
 
     # Frequencies
-    Ω1, Ω2 = OrbitalElements.ComputeFrequenciesAE(model,a,e,params.Orbitalparams)
+    Ω1, Ω2 = OrbitalElements.frequencies_from_ae(a,e,model,params.Orbitalparams)
 
     # Basis FT
 
@@ -296,13 +297,13 @@ end
 function MakeWmatUV(model::OrbitalElements.Potential,
                     n1::Int64,n2::Int64,
                     tabu::Vector{Float64},
-                    basisFT::BasisFTtype,
+                    basisFT::FourierTransformedBasis,
                     params::LinearParameters)
 
     @assert length(tabu) == params.Ku "LinearResponse.WMat.MakeWmatUV: tabu length is not Ku."
 
     Orbitalparams = params.Orbitalparams
-    Ω₀ = Orbitalparams.Ω₀
+    Ω₀ = OrbitalElements.frequency_scale(model)
 
     # allocate the results matrices
     Wdata = WMatdataCreate(model,n1,n2,params)
@@ -317,7 +318,8 @@ function MakeWmatUV(model::OrbitalElements.Potential,
         (params.VERBOSE > 2) && println("LinearResponse.WMat.MakeWMat: on step $kuval of $Ku: u=$uval.")
 
         # get the corresponding v boundary values
-        vmin,vmax = OrbitalElements.FindVminVmax(uval,n1,n2,model,ωmin,ωmax,Orbitalparams)
+        resonance = OrbitalElements.Resonance(n1,n2,model,Orbitalparams)
+        vmin,vmax = OrbitalElements.v_boundaries(uval,resonance,model,Orbitalparams)
 
         # saving them
         Wdata.tabvminmax[1,kuval], Wdata.tabvminmax[2,kuval] = vmin, vmax
@@ -336,11 +338,11 @@ function MakeWmatUV(model::OrbitalElements.Potential,
             # (u,v) -> (a,e)
             ####
             # (u,v) -> (α,β)
-            α,β = OrbitalElements.αβFromUV(uval,vval,n1,n2,ωmin,ωmax)
+            α,β = OrbitalElements.αβ_from_uv(uval,vval,resonance)
             # (α,β) -> (Ω1,Ω2)
-            Ω₁, Ω₂ = OrbitalElements.FrequenciesFromαβ(α,β,Ω₀)
+            Ω₁, Ω₂ = OrbitalElements.frequencies_from_αβ(α,β,Ω₀)
             # (Ω1,Ω2) -> (a,e)
-            a,e = OrbitalElements.ComputeAEFromFrequencies(model,Ω₁,Ω₂,Orbitalparams)
+            a,e = OrbitalElements.ae_from_frequencies(Ω₁,Ω₂,model,Orbitalparams)
 
             (params.VERBOSE > 2) && print("v=$kvval,o1=$Ω1,o2=$Ω2;")
 
@@ -351,10 +353,10 @@ function MakeWmatUV(model::OrbitalElements.Potential,
             # save (a,e) values for later
             Wdata.tabAE[1,kvval,kuval], Wdata.tabAE[2,kvval,kuval] = a, e
             # save (E,L) values for later
-            Wdata.tabEL[1,kvval,kuval], Wdata.tabEL[2,kvval,kuval] = OrbitalElements.ELFromAE(model,a,e,Orbitalparams)
+            Wdata.tabEL[1,kvval,kuval], Wdata.tabEL[2,kvval,kuval] = OrbitalElements.EL_from_ae(a,e,model,Orbitalparams)
 
             # compute the Jacobian of the (α,β) ↦ (E,L) mapping here. a little more expensive, but savings in the long run
-            Wdata.tabJ[kvval,kuval] = OrbitalElements.JacαβToELAE(model,a,e,Orbitalparams)
+            Wdata.tabJ[kvval,kuval] = OrbitalElements.ae_to_EL_jacobian(a,e,model,Orbitalparams)/OrbitalElements.ae_to_αβ_jacobian(a,e,model,Orbitalparams)
 
             # Compute W(u,v) for every basis element using RK4 scheme
             WBasisFT(model,a,e,Ω₁,Ω₂,n1,n2,basisFT,params)
@@ -427,6 +429,7 @@ function RunWmat(model::OrbitalElements.Potential,
             # Mappings parameters
             write(file, "omgmin",Wdata.ωmin)
             write(file, "omgmax",Wdata.ωmax)
+            write(file,"Ω₀",Wdata.Ω₀)
             write(file, "tabvminmax",Wdata.tabvminmax)
             # Mappings
             write(file, "UVmat",Wdata.tabUV)