Skip to content

added debye3 function for eventual T-dependent heat capacities #4769

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
tburch-LSS wants to merge 1 commit into
base: master
Choose a base branch
from
+83 −0
Open

added debye3 function for eventual T-dependent heat capacities #4769

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
83 changes: 83 additions & 0 deletions Modelica/Math/Special.mo
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -426,6 +426,89 @@ sinc(3) // = 0.0470400026866224
</pre></blockquote>
</html>"));
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 = "<html><head></head><body><h4>Syntax</h4>
<blockquote><pre>Special.<b>debye3</b>(u);
</pre></blockquote>

<h4>Description</h4>
<p>
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.</p>

<p>
For more details, see <a href=\"http://en.wikipedia.org/wiki/Debye_function\">Wikipedia</a>.
</p>

<h4>Example</h4>
<blockquote><pre>debye3(0) // = 1
debye3(5.e-4) // = 0.999813
debye3(1) // = 0.674416
debye3(100) // = 1.94818e-5</pre></blockquote>

<h4>See also</h4>
<p>
<a href=\"modelica://Modelica.Math.Special.erfc\">erfc</a>.
</p>
</body></html>", revisions = "<html><head></head><body><table border=\"1\" cellspacing=\"0\" cellpadding=\"2\">
<tbody><tr><th>Date</th> <th align=\"left\">Description</th></tr>

<tr><td>March 22, 2026</td>
<td>

<table border=\"0\">
<tbody><tr><td>
Initial version implemented by T. Burch.<br></td></tr></tbody></table>
</td></tr>

</tbody></table>
</body></html>"));
end debye3;

package Internal
"Internal utility functions that should not be directly utilized by the user"
Expand Down