fixes docs 10- issue (#72)

* fixes docs 10- issue

- Remove physical constraint validation in Invariants function
- Use a default constant for border point generation in Dalitz plot
- Standardize type parameter formatting in type definitions

* remote !isphysical test

* missing docstrings

* docs: Improve docstrings and formatting in phase_space.jl and mandelstam.jl

- Add proper docstring formatting for `σjofk` and `σiofk` functions
- Standardize type parameter spacing in `MandelstamTuple`
- Add cross-reference to related function using `[@ref]`

* another ref

* docs: Improve docstring for isphysical and related functions

- Add return type annotations to function signatures
- Reposition helper functions for better readability
- Update docstring formatting and cross-references

* docs: Enhance docstrings with more comprehensive examples and formatting

- Add multiple examples to `Kallen`, `sqrtKallenFact`, `Kibble`, `breakup`, `aligned_four_vectors`, and `isphysical` functions
- Improve docstring readability with `jldoctest` blocks
- Include practical usage scenarios and edge case demonstrations
- Update cross-references and formatting for consistency

* typo in docstring

* add setup

* ; setup = :(using ThreeBodyDecays)

* fixed isphysical test for n_provided == 3

* cspell, formatter

Author: mmikhasenko

.cspell/project.txt

@@ -19,6 +19,7 @@ helicities
 helicity
 imag
 issequential
+jldoctest
 JPAC
 Källén
 latexstring
@@ -46,6 +47,7 @@ sinθ
 subchannel
 sumrules
 typeof
+unphysical
 Utopic
 vcat
 wignerd

src/dalitz.jl

@@ -21,13 +21,13 @@ polardalitz2invariants(θ, expansion_point::Tuple) =
     ) .+ expansion_point
 
 function invariants(σx, σy; iσx, iσy, ms)
-    iσx == 1 && iσy == 2 && return Invariants(ms, σ1 = σx, σ2 = σy)
-    iσx == 2 && iσy == 1 && return Invariants(ms, σ1 = σy, σ2 = σx)
-    iσx == 1 && iσy == 3 && return Invariants(ms, σ1 = σx, σ3 = σy)
-    iσx == 3 && iσy == 1 && return Invariants(ms, σ1 = σy, σ3 = σx)
-    iσx == 3 && iσy == 2 && return Invariants(ms, σ2 = σy, σ3 = σx)
-    return Invariants(ms, σ2 = σx, σ3 = σy) # σx = 2 && σy = 3
-    # error()
+    iσx == 1 && iσy == 2 && return Invariants(ms; σ1 = σx, σ2 = σy)
+    iσx == 2 && iσy == 1 && return Invariants(ms; σ2 = σx, σ1 = σy)
+    iσx == 1 && iσy == 3 && return Invariants(ms; σ1 = σx, σ3 = σy)
+    iσx == 3 && iσy == 1 && return Invariants(ms; σ3 = σx, σ1 = σy)
+    iσx == 3 && iσy == 2 && return Invariants(ms; σ3 = σx, σ2 = σy)
+    # remaining case σx = 2 && σy = 3
+    return Invariants(ms; σ2 = σx, σ3 = σy)
 end
 
 """
@@ -42,7 +42,7 @@ Calculate the border of the Dalitz plot.
 # Returns
 - `Vector{MandelstamTuple{T}}`: Points on the border
 """
-function border(ms::MassTuple{T}; Nx::Int = 300) where {T}
+function border(ms::MassTuple{T}; Nx::Int = DEFAULT_BORDER_POINTS) where {T}
     # Calculate a physically valid expansion point
     expansion_point = let
         f = 0.5
@@ -93,8 +93,8 @@ end
     iσx = 1,
     iσy = 2,
     grid_size = 100,
-    xlims = lims(iσx, ms),
-    ylims = lims(iσy, ms),
+    xlims = lims(ms; k = iσx),
+    ylims = lims(ms; k = iσy),
 )
     #
     σxv = range(xlims..., length = grid_size + 1) |> shift_by_half

src/mandelstam.jl

@@ -41,11 +41,11 @@ function Invariants(
         else
             σ2 = sum(ms^2) - σ3 - σ1
         end
-    end
-
-    # Validate the invariants
-    if !isphysical((σ1, σ2, σ3), ms)
-        error("The provided invariants violate physical constraints")
+    elseif n_provided == 3
+        # Check if the invariants are physical
+        if !isphysical((σ1, σ2, σ3), ms)
+            error("The invariants violate mass constraints")
+        end
     end
 
     return MandelstamTuple{T}((σ1, σ2, σ3))

src/mass_types.jl

@@ -33,15 +33,31 @@ ThreeBodyMasses(; m1, m2, m3, m0) = ThreeBodyMasses(m1, m2, m3; m0)
 """
     lims(ms::MassTuple; k::Int)
     lims(k::Int, ms::MassTuple)
+    lims1(ms::MassTuple)
+    lims2(ms::MassTuple)
+    lims3(ms::MassTuple)
 
-Calculate the kinematic limits for the k-th Mandelstam variable.
+Calculate the kinematic limits (boundaries) for the Mandelstam variables σᵢ in a three-body decay.
+For each pair of particles (i,j), the invariant mass squared σₖ = (pᵢ + pⱼ)² must lie within physical limits:
+- Lower bound: (mᵢ + mⱼ)² (threshold for producing particles i and j)
+- Upper bound: (M - mₖ)² (maximum energy available when particle k recoils)
 
 # Arguments
-- `ms::MassTuple`: Masses of the system
-- `k::Int`: Index of the variable (1, 2, or 3)
+- `ms`: Tuple of masses (m₁,m₂,m₃,M)
+- `k`: Index specifying which pair of particles (1,2,3)
 
 # Returns
-- `Tuple{T,T}`: Lower and upper limits for the variable
+- Tuple (min,max) of the allowed range for σₖ
+
+# Example
+```julia
+ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0)
+lims1(ms)  # limits for σ₁ = (p₂ + p₃)²
+lims2(ms)  # limits for σ₂ = (p₃ + p₁)²
+lims3(ms)  # limits for σ₃ = (p₁ + p₂)²
+```
+
+See also [`isphysical`](@ref), [`Kibble`](@ref).
 """
 function lims(ms::MassTuple; k::Int)
     i, j = ij_from_k(k)

src/phase_space.jl

@@ -10,17 +10,80 @@ This function appears frequently in relativistic kinematics calculations.
 # Returns
 - The value of the Källén function
 
-# Example
-```julia
-# For a decay A -> B + C, calculate λ(mA², mB², mC²)
-Kallen(10.0, 1.0, 1.0)  # returns 64.0
+# Examples
+```jldoctest; setup = :(using ThreeBodyDecays)
+julia> Kallen(10.0, 1.0, 1.0)  # For a decay A -> B + C, calculate λ(mA², mB², mC²)
+60.0
+
+julia> Kallen(5.0, 1.0, 1.0)  # Another example with different masses
+5.0
 ```
 
 See also [`sqrtKallenFact`](@ref), [`Kibble`](@ref).
 """
 Kallen(x, y, z) = x^2 + y^2 + z^2 - 2x * y - 2y * z - 2z * x
+
+"""
+    sqrtKallenFact(a, b, c)
+
+Calculate the square root of the Källén function λ(a,b,c) in factorized form.
+The Källén function is defined as λ(a,b,c) = a² + b² + c² - 2ab - 2bc - 2ca.
+This function returns √λ(a,b,c) in a factorized form to improve numerical stability:
+√(a-(b+c)) √(a-(b-c)) √(a+(b+c)) √(a+(b-c)).
+
+# Arguments
+- `a`, `b`, `c`: Input values for the Källén function
+
+# Returns
+- Square root of the Källén function in factorized form
+
+# Examples
+```jldoctest; setup = :(using ThreeBodyDecays)
+julia> sqrtKallenFact(10.0, 1.0, 1.0)  # √λ(10,1,1) = √9600
+97.97958971132714
+
+julia> sqrtKallenFact(4.0, 1.0, 1.0)  # √λ(4,1,1) = √192
+13.856406460551018
+```
+
+See also [`Kallen`](@ref), [`Kibble`](@ref).
+"""
 sqrtKallenFact(a, b, c) =
     sqrt(a - (b + c)) * sqrt(a - (b - c)) * sqrt(a + (b + c)) * sqrt(a + (b - c))
+
+"""
+    Kibble(σs, msq)
+
+Calculate the Kibble function φ(σ₁,σ₂,σ₃) for three-body decays, which is defined in terms of Källén functions.
+The Kibble function determines the physical boundaries of the three-body phase space.
+
+The function is calculated as:
+φ(σ₁,σ₂,σ₃) = λ(λ₁,λ₂,λ₃)
+where λᵢ = λ(M², mᵢ², σᵢ) are Källén functions of the total mass squared M²,
+the i-th particle mass squared mᵢ², and the corresponding Mandelstam variable σᵢ.
+
+# Arguments
+- `σs`: Tuple of Mandelstam variables (σ₁,σ₂,σ₃)
+- `msq`: Tuple of squared masses (m₁²,m₂²,m₃²,M²)
+
+# Returns
+- Value of the Kibble function
+
+# Examples
+```jldoctest; setup = :(using ThreeBodyDecays)
+julia> msq = (1.0, 1.0, 1.0, 16.0);  # squared masses: m₁², m₂², m₃², M²
+
+julia> σs = (2.0, 2.0, 2.0);         # Mandelstam variables
+
+julia> Kibble(σs, msq)               # returns the Kibble function value
+-77763.0
+
+julia> Kibble(σs, msq) < 0           # physical configuration if negative
+true
+```
+
+See also [`Kallen`](@ref), [`isphysical`](@ref).
+"""
 Kibble(σs, msq) = Kallen(
     Kallen(msq[4], msq[1], σs[1]),
     Kallen(msq[4], msq[2], σs[2]),
@@ -29,15 +92,15 @@ Kibble(σs, msq) = Kallen(
 
 
 """
-σjofk(z,σi,msq; k::Int)
+    σjofk(z,σi,msq; k::Int)
 
 Computes invariant σj = (p0-pj)² from
 the scattering angle z=cosθij in the rest from of (i,j),
 given the mass of the system m(i,j)² = σk
 
 Explicit forms: `σ3of1`, `σ1of2`, `σ2of3`.
 
-See also `σiofk(z,σj,msq; k)`
+See also [`σiofk`](@ref)
 """
 function σjofk(z, σk, msq; k::Int)
     i, j = ij_from_k(k)
@@ -56,15 +119,15 @@ end
 
 
 """
-σiofk(k,z,σj,msq)
+    σiofk(k,z,σj,msq)
 
 Computes invariant σi = (p0 - pi)² from
 the scattering angle z=cosθij in the rest from of (i,j),
 given the mass of the system m(i,j)² = σk
 
 Explicit forms: `σ3of2`, `σ1of3`, `σ2of1`.
 
-See also `σjofk(z,σk,msq; k)`
+See also [`σjofk`](@ref)
 """
 function σiofk(z, σk, msq; k::Int)
     σj = σjofk(z, σk, msq; k)
@@ -120,15 +183,18 @@ Calculate the breakup momentum for a two-body system.
 # Returns
 - Break-up momentum magnitude
 
-# Example
-```julia
+# Examples
+```jldoctest; setup = :(using ThreeBodyDecays)
 julia> ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0);
 
 julia> p = breakup_Rk(2.5^2, ms; k=1)
 0.897587913521567
 
 julia> q = breakup_ij(2.5^2, ms; k=1)
 0.75
+
+julia> breakup(4.0, 2.5, 1.0)  # direct use of breakup function
+0.897587913521567
 ```
 """
 breakup(m, m1, m2) =
@@ -155,15 +221,22 @@ Computes the four-momenta of the three particles in the center of momentum frame
 - `k`: Index of the particle to be aligned with the -z-axis.
 
 ## Returns
-
 A tuple of three four-momenta in the form of (px, py, pz, E).
 
-## Example
-````julia
-ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0)
-σs = x2σs([0.5, 0.5], ms; k=2)
-p1, p2, p3 = aligned_four_vectors(σs, ms; k=1)
-````
+## Examples
+```jldoctest; setup = :(using ThreeBodyDecays)
+julia> ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0);
+
+julia> σs = x2σs([0.5, 0.5], ms; k=2);
+
+julia> p1, p2, p3 = aligned_four_vectors(σs, ms; k=1);
+
+julia> p1[3] ≈ -breakup_Rk(σs[1], ms; k=1)  # first particle aligned with -z axis
+true
+
+julia> all(p -> abs(p[2]) < 1e-10, (p1, p2, p3))  # all momenta in x-z plane
+true
+```
 """
 function aligned_four_vectors(σs, ms; k::Int)
     #
@@ -193,24 +266,45 @@ function aligned_four_vectors(σs, ms; k::Int)
     return circleorigin(-k, (pi, pj, pk))
 end
 
+inrange(x, r) = r[1] < x < r[2]
+inphrange(σs, ms::MassTuple) = isphysical(σs, ms)
+
 """
-    inrange(x, r)
-    inphrange(σs::MandelstamTuple, ms::MassTuple)
-    isphysical(σs::MandelstamTuple, ms::MassTuple)
+    isphysical(σs::MandelstamTuple, ms::MassTuple) -> Bool
+    inphrange(σs::MandelstamTuple, ms::MassTuple) -> Bool
+    inrange(x, r) -> Bool
 
-Check if values are within physical ranges.
+Check if a set of Mandelstam variables represents a physically valid configuration in the three-body phase space.
+A configuration is physical if:
+1. The Kibble function is negative (φ(σ₁,σ₂,σ₃) < 0)
+2. All Mandelstam variables are within their kinematically allowed ranges
 
 # Arguments
-- `x`: Value to check
-- `r`: Range to check against
-- `σs`: Mandelstam variables
-- `ms`: Masses of the system
+- `σs::MandelstamTuple`: Tuple of Mandelstam variables (σ₁,σ₂,σ₃)
+- `ms::MassTuple`: Tuple of masses (m₁,m₂,m₃,M)
+- `x`: Value to check (for `inrange`)
+- `r`: Range tuple (min,max) to check against (for `inrange`)
 
 # Returns
-- `Bool`: Whether the values are physical
+- `Bool`: `true` if the configuration is physical, `false` otherwise
+
+# Examples
+```jldoctest; setup = :(using ThreeBodyDecays)
+julia> ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0);
+
+julia> σs = (6.25, 6.25, 6.5);
+
+julia> isphysical(σs, ms)  # checks if this configuration is physically possible
+true
+
+julia> σs_unphysical = (20.0, 20.0, 20.0);  # values outside allowed ranges
+
+julia> isphysical(σs_unphysical, ms)  # this configuration is not physical
+false
+```
+
+See also [`Kibble`](@ref), [`lims`](@ref).
 """
-inrange(x, r) = r[1] < x < r[2]
-inphrange(σs, ms::MassTuple) = isphysical(σs, ms)
 isphysical(σs, ms::MassTuple) =
     Kibble(σs, ms^2) < 0 &&
     inrange(σs[1], lims1(ms)) &&