From 9dde876ddb18fe4f37090dbb7e2e4937ecc5b441 Mon Sep 17 00:00:00 2001 From: Tommy Burch Date: Mon, 23 Mar 2026 14:03:19 +0100 Subject: [PATCH] added debye3 function for eventual T-dependent heat capacities --- Modelica/Math/Special.mo | 83 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) diff --git a/Modelica/Math/Special.mo b/Modelica/Math/Special.mo index 7988a2ef71..a10ce59162 100644 --- a/Modelica/Math/Special.mo +++ b/Modelica/Math/Special.mo @@ -426,6 +426,89 @@ sinc(3) // = 0.0470400026866224 ")); end sinc; + + function debye3 "3rd-order Debye function debye3(u) = 3/u^3*Integral_0_u t^3/(exp(t)-1)*dt" + extends Modelica.Icons.Function; + input Real u "Input argument"; + output Real y "= 3/u^3*Integral_0_u t^3/(exp(t)-1)*dt"; + // + function integrand + extends Modelica.Math.Nonlinear.Interfaces.partialScalarFunction; + algorithm + y := u*u*u / (exp(u) - 1.0); + end integrand; + // + protected + Real z; + Real u3; + Real sum; + constant Real eps = 1.e-3; + constant Real eps3 = eps*eps*eps; + constant Real uu = 15.0; + constant Real uu3 = uu*uu*uu; + constant Real rel_err = 1.e-8; + constant Real Imax = 19.4818182068004875; // u^3 debye3(u) / 3 ~ pi^4 / 15 for u-->inf + algorithm + if u < 0 then + // domain error: return debye3(0) result + y := 1.0; + elseif u < eps then + // small u expansion + y := 1.0 - 3.0*u/8.0 + u*u/20.0; + elseif u <= 0.5*uu then + // forward Lobatto quadrature + u3 := u*u*u; + z := Modelica.Math.Nonlinear.quadratureLobatto(function integrand(), eps, u, rel_err); + y := 3.0/u3 * (z + eps3/3.0 * (1.0 - 3.0*eps/8.0 + eps*eps/20.0)); + elseif u <= uu then + // backward Lobatto quadrature + u3 := u*u*u; + z := Modelica.Math.Nonlinear.quadratureLobatto(function integrand(), u, uu, rel_err); + sum := 6.0 + 6.0*uu + 3.0*uu*uu + uu3; + y := 3.0/u3 * (-z + Imax/3.0 - sum*exp(-uu)); + else + // debye3(u->inf) minus large-u-to-inf expansion + u3 := u*u*u; + sum := 6.0 + 6.0*u + 3.0*u*u + u3; + y := ( Imax - 3.0*sum*exp(-u) ) / u3; + end if; + annotation( + Documentation(info = "

Syntax

+
Special.debye3(u);
+
+ +

Description

+

+This function computes the 3rd order Debye function debye3(u) = 3/u^3*Integral_0_u t^3/(exp(t)-1)*dt numerically with a relative precision of 1e-8. The implementation uses the Lobatto quadrature routine available in the MSL, bracketed by approximations via expansions for small and large arguments.

+ +

+For more details, see Wikipedia. +

+ +

Example

+
debye3(0)     // = 1
+debye3(5.e-4) // = 0.999813
+debye3(1)     // = 0.674416
+debye3(100)   // = 1.94818e-5
+ +

See also

+

+erfc. +

+", revisions = " + + + + + +
Date Description
March 22, 2026 + + +
+ Initial version implemented by T. Burch.
+
+")); + end debye3; package Internal "Internal utility functions that should not be directly utilized by the user"