diff --git a/docs/Project.toml b/docs/Project.toml
index dfe53045d..0685a35b9 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38"
 
 [compat]
 Documenter = "0.27"
diff --git a/docs/make.jl b/docs/make.jl
index 35586b7f4..ab4b224d7 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -6,7 +6,8 @@ makedocs(;
     sitename="GeoParams.jl",
     authors="Boris Kaus and contributors",
     modules=[GeoParams],
-    format=Documenter.HTML(; prettyurls=get(ENV, "CI", nothing) == "true"), # easier local build
+    format=Documenter.HTML(; assets = ["assets/favicon.ico"],
+    prettyurls=get(ENV, "CI", nothing) == "true"), # easier local build
     pages=[
         "Home" => "index.md",
         "User Guide" => Any[
diff --git a/docs/src/assets/favicon.ico b/docs/src/assets/favicon.ico
new file mode 100644
index 000000000..9b9764f45
Binary files /dev/null and b/docs/src/assets/favicon.ico differ
diff --git a/docs/src/man/plotting.md b/docs/src/man/plotting.md
index 913ec7f71..667995669 100644
--- a/docs/src/man/plotting.md
+++ b/docs/src/man/plotting.md
@@ -1,7 +1,7 @@
 # Plotting
 
 We provide a number of plotting routines. Note that these plotting routines only become available once the `GLMakie` package is loaded.
-### Plot CreepLaws 
+### Plot CreepLaws
 
 ```@docs
 GeoParams.PlotStressStrainrate_CreepLaw
diff --git a/docs/src/man/strengthenvelope.md b/docs/src/man/strengthenvelope.md
index 204a3a34f..8a28b8a79 100644
--- a/docs/src/man/strengthenvelope.md
+++ b/docs/src/man/strengthenvelope.md
@@ -10,5 +10,5 @@ GeoParams.CompTempStruct
 GeoParams.LithPres
 GeoParams.StrengthEnvelopeComp
 GeoParams.StrengthEnvelopePlot
-![XMasTree](../assets/XMasTree.png)
 ```
+![XMasTree](../assets/XMasTree.png)
diff --git a/src/MeltFraction/MeltingParameterization.jl b/src/MeltFraction/MeltingParameterization.jl
index ae271389b..1dfae5698 100644
--- a/src/MeltFraction/MeltingParameterization.jl
+++ b/src/MeltFraction/MeltingParameterization.jl
@@ -7,7 +7,7 @@ using Parameters, LaTeXStrings, Unitful, ForwardDiff
 using ..Units
 using GeoParams:
     AbstractMaterialParam, PhaseDiagram_LookupTable, AbstractMaterialParamsStruct
-import Base.show, GeoParams.param_info 
+import Base.show, GeoParams.param_info
 using ..MaterialParameters: MaterialParamsInfo
 using Setfield # allows modifying fields in immutable struct
 
@@ -105,7 +105,7 @@ Implements the a smooth 3rd order T-dependent melting parameterisation (as used
 ```math
     \\phi_{melt} = {1.0 \\over (1.0 + e^\\theta)}
 ```
-Note that T is in Kelvin. 
+Note that T is in Kelvin.
 
 As default parameters we employ:
 ```math
@@ -115,7 +115,7 @@ which gives a reasonable fit to experimental data for basalt.
 
 Data for rhyolite are:
 ```math
-a=3043.0,  b=-10552.0, c=12204.9, d = -4709.0 
+a=3043.0,  b=-10552.0, c=12204.9, d = -4709.0
 ```
 
 ![MeltingParam_Smooth3rdOrder](./assets/img/MeltingParam_Melnik.png)
@@ -125,28 +125,28 @@ References
 ====
 """
 @with_kw_noshow struct MeltingParam_Smooth3rdOrder{T,U,U1} <: AbstractMeltingParam{T}
-    a::GeoUnit{T,U} =  517.9NoUnits
+    a::GeoUnit{T,U} = 517.9NoUnits
     b::GeoUnit{T,U} = -1619.0NoUnits
     c::GeoUnit{T,U} = 1699.0NoUnits
     d::GeoUnit{T,U} = -597.4NoUnits
-    T0::GeoUnit{T,U1} = 273.15K 
+    T0::GeoUnit{T,U1} = 273.15K
     Tchar::GeoUnit{T,U1} = 1000K # normalization
     apply_bounds::Bool = true
 end
-MeltingParam_Smooth3rdOrder(args...) = MeltingParam_Smooth3rdOrder(convert.(GeoUnit, args)...)
+function MeltingParam_Smooth3rdOrder(args...)
+    return MeltingParam_Smooth3rdOrder(convert.(GeoUnit, args)...)
+end
 
 function param_info(s::MeltingParam_Smooth3rdOrder) # info about the struct
-    return MaterialParamsInfo(;
-        Equation=L"\phi = f(T)"
-    )
+    return MaterialParamsInfo(; Equation=L"\phi = f(T)")
 end
 
 # Calculation routines
 function (p::MeltingParam_Smooth3rdOrder)(; T, kwargs...)
     @unpack_val a, b, c, d, Tchar, T0 = p
-    x = (T - T0)/ Tchar
-    
-    θ = min(evalpoly(x, (a, b, c, d)),200.0)
+    x = (T - T0) / Tchar
+
+    θ = min(evalpoly(x, (a, b, c, d)), 200.0)
 
     ϕ = inv(1.0 + exp(θ))
     return ϕ
@@ -154,20 +154,16 @@ end
 
 function compute_dϕdT(p::MeltingParam_Smooth3rdOrder; T, kwargs...)
     @unpack_val a, b, c, d, Tchar, T0 = p
-    
-    x =  (T - T0) / Tchar
-    θ = min(evalpoly(x, (a, b, c, d)),200.0)
 
-    dϕdT = -exp(θ)*(1 / (1.0 + exp(θ))^2 )*(b / Tchar + (3d*x^2) / Tchar + (2(T - T0)*c) / (Tchar^2))
+    x = (T - T0) / Tchar
+    θ = min(evalpoly(x, (a, b, c, d)), 200.0)
 
-    return dϕdT
-end
+    dϕdT =
+        -exp(θ) *
+        (1 / (1.0 + exp(θ))^2) *
+        (b / Tchar + (3d * x^2) / Tchar + (2(T - T0) * c) / (Tchar^2))
 
-function foo(T,T0,Tchar,a,b,c,d)
-    x =  (T - T0) / Tchar
-    θ = a + b * x + c * x^2 + d * x^3;
-    ϕ = inv(1.0 + exp(θ))
-    return x^2
+    return dϕdT
 end
 
 # Print info
@@ -217,7 +213,6 @@ function param_info(s::MeltingParam_5thOrder) # info about the struct
     return MaterialParamsInfo(; Equation=L"\phi = aT^5 + bT^4 + cT^3 + dT^2 + eT + f")
 end
 
-
 # Calculation routines
 function (p::MeltingParam_5thOrder)(; T, kwargs...)
     @unpack_val a, b, c, d, e, f, T_s, T_l = p
@@ -235,14 +230,14 @@ function (p::MeltingParam_5thOrder)(; T, kwargs...)
     return ϕ
 end
 
-compute_dϕdT(p::MeltingParam_5thOrder, T, kwargs...) = compute_dϕdT(p; T, kwargs...) 
+compute_dϕdT(p::MeltingParam_5thOrder, T, kwargs...) = compute_dϕdT(p; T, kwargs...)
 
 function compute_dϕdT(p::MeltingParam_5thOrder; T, kwargs...)
     @unpack_val a, b, c, d, e, T_s, T_l = p
 
-    coeffs = e, 2*d, 3*c, 4*b, 5*a
+    coeffs = e, 2 * d, 3 * c, 4 * b, 5 * a
     dϕdT = evalpoly(T, coeffs)
-    
+
     if p.apply_bounds && (T < T_s || T > T_l)
         dϕdT = 0.0
     end
@@ -317,7 +312,7 @@ end
 function compute_dϕdT(p::MeltingParam_4thOrder; T, kwargs...)
     @unpack_val b, c, d, e, T_s, T_l = p
 
-    coeffs = e, 2*d, 3*c, 4*b
+    coeffs = e, 2 * d, 3 * c, 4 * b
     dϕdT = evalpoly(T, coeffs)
     if p.apply_bounds && (T < T_s || T > T_l)
         dϕdT = 0.0
@@ -392,7 +387,7 @@ function (p::MeltingParam_Quadratic)(; T, kwargs...)
 end
 
 function compute_dϕdT(p::MeltingParam_Quadratic; T, kwargs...)
-     @unpack_val T_s, T_l = p
+    @unpack_val T_s, T_l = p
 
     dϕdT = (2T_l - 2T) / ((T_l - T_s)^2)
     if p.apply_bounds && (T > T_l || T < T_s)
@@ -525,10 +520,10 @@ end
 
 Stores a vector with melt fraction that can be retrieved by providing an `index`
 """
-struct Vector_MeltingParam{_T, V <: AbstractVector} <: AbstractMeltingParam{_T}
+struct Vector_MeltingParam{_T,V<:AbstractVector} <: AbstractMeltingParam{_T}
     ϕ::V       # melt fraction
 end
-Vector_MeltingParam(; ϕ=Vector{Float64}()) = Vector_MeltingParam{eltype(ϕ), typeof(ϕ)}(ϕ)
+Vector_MeltingParam(; ϕ=Vector{Float64}()) = Vector_MeltingParam{eltype(ϕ),typeof(ϕ)}(ϕ)
 
 function param_info(s::Vector_MeltingParam) # info about the struct
     return MaterialParamsInfo(; Equation=L"\phi from a precomputed vector")
@@ -549,7 +544,6 @@ function show(io::IO, g::Vector_MeltingParam)
 end
 #-------------------------------------------------------------------------
 
-
 # Smooth melting function ------------------------------------------------
 
 """
@@ -627,7 +621,7 @@ function (param::SmoothMelting)(; T, kwargs...)
     ϕ = param.p(; T, kwargs...)     # Melt fraction computed in usual manner
 
     T_s = param.p.T_s
-    H_s =inv(1.0 + exp(-2 * k_sol * (T - T_s - (2 / k_sol))))
+    H_s = inv(1.0 + exp(-2 * k_sol * (T - T_s - (2 / k_sol))))
 
     T_l = param.p.T_l
     H_l = 1.0 - inv(1.0 + exp(-2 * k_liq * (T - T_l + (2 / k_liq))))
@@ -638,7 +632,6 @@ function (param::SmoothMelting)(; T, kwargs...)
     return ϕ
 end
 
-
 function compute_dϕdT(param::SmoothMelting; T, kwargs...)
     @unpack_val k_sol, k_liq = param
 
@@ -674,7 +667,6 @@ function show(io::IO, g::SmoothMelting)
 end
 #-------------------------------------------------------------------------
 
-
 """
     compute_meltfraction(P,T, p::AbstractPhaseDiagramsStruct)
 
@@ -743,7 +735,7 @@ for myType in (
     :MeltingParam_Quadratic,
     :MeltingParam_Assimilation,
     :Vector_MeltingParam,
-    :SmoothMelting, 
+    :SmoothMelting,
 )
     @eval begin
         (p::$(myType))(args) = p(; args...)
@@ -774,14 +766,16 @@ end
 
 Computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
 """
-compute_meltfraction(args::Vararg{Any, N}) where N = compute_param(compute_meltfraction, args...)
+compute_meltfraction(args::Vararg{Any,N}) where {N} =
+    clamp(compute_param(compute_meltfraction, args...), 0e0, 1e0)
 
 """
     compute_meltfraction(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})
 
 In-place computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
 """
-compute_meltfraction!(args::Vararg{Any, N}) where N = compute_param!(compute_meltfraction, args...)
+compute_meltfraction!(args::Vararg{Any,N}) where {N} =
+    compute_param!(compute_meltfraction, args...)
 
 """
     compute_meltfraction(ϕ::AbstractArray{<:AbstractFloat}, PhaseRatios::Union{NTuple{N,T}, SVector{N,T}}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})
@@ -789,7 +783,8 @@ compute_meltfraction!(args::Vararg{Any, N}) where N = compute_param!(compute_mel
 Computation of melt fraction ϕ for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
 This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)
 """
-compute_meltfraction_ratio(args::Vararg{Any, N}) where N = compute_param_times_frac(compute_meltfraction, args...)
+compute_meltfraction_ratio(args::Vararg{Any,N}) where {N} =
+    compute_param_times_frac(compute_meltfraction, args...)
 
 """
     ϕ = compute_dϕdT(Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})
@@ -797,7 +792,7 @@ compute_meltfraction_ratio(args::Vararg{Any, N}) where N = compute_param_times_f
 Computates the derivative of melt fraction ϕ versus temperature `T` for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
 This is employed in computing latent heat terms in an implicit manner, for example
 """
-compute_dϕdT(args::Vararg{Any, N}) where N = compute_param(compute_dϕdT, args...)
+compute_dϕdT(args::Vararg{Any,N}) where {N} = compute_param(compute_dϕdT, args...)
 
 """
     compute_dϕdT!(ϕ::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},T::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})
@@ -805,6 +800,6 @@ compute_dϕdT(args::Vararg{Any, N}) where N = compute_param(compute_dϕdT, args.
 Computates the derivative of melt fraction `ϕ` versus temperature `T`, ``{\\partial \\phi} \\over {\\partial T}`` for the whole domain and all phases, in case an array with phase properties `MatParam` is provided, along with `P` and `T` arrays.
 This is employed, for example, in computing latent heat terms in an implicit manner.
 """
-compute_dϕdT!(args::Vararg{Any, N}) where N = compute_param!(compute_dϕdT, args...)
+compute_dϕdT!(args::Vararg{Any,N}) where {N} = compute_param!(compute_dϕdT, args...)
 
 end