1 files changed, 100 insertions, 0 deletions
diff --git a/source/luametatex/source/libraries/libcerf/width.c b/source/luametatex/source/libraries/libcerf/width.c
new file mode 100644
index 000000000..a844377ff
--- /dev/null
+++ b/source/luametatex/source/libraries/libcerf/width.c
@@ -0,0 +1,100 @@
+/* Library libcerf:
+ *   Compute complex error functions, based on a new implementation of
+ *   Faddeeva's w_of_z. Also provide Dawson and Voigt functions.
+ *
+ * File width.c:
+ *   Computate voigt_hwhm, using Newton's iteration.
+ *
+ * Copyright:
+ *   (C) 2018 Forschungszentrum Jülich GmbH
+ *
+ * Licence:
+ *   Permission is hereby granted, free of charge, to any person obtaining
+ *   a copy of this software and associated documentation files (the
+ *   "Software"), to deal in the Software without restriction, including
+ *   without limitation the rights to use, copy, modify, merge, publish,
+ *   distribute, sublicense, and/or sell copies of the Software, and to
+ *   permit persons to whom the Software is furnished to do so, subject to
+ *   the following conditions:
+ *
+ *   The above copyright notice and this permission notice shall be
+ *   included in all copies or substantial portions of the Software.
+ *
+ *   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ *   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ *   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ *   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ *   LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ *   OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ *   WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ *
+ * Authors:
+ *   Joachim Wuttke, Forschungszentrum Jülich, 2018
+ *
+ * Website:
+ *   http://apps.jcns.fz-juelich.de/libcerf
+ *
+ * Revision history:
+ *   ../CHANGELOG
+ *
+ * Man pages:
+ *   voigt_fwhm(3)
+ */
+
+/*
+
+    This file is patched by Hans Hagen for usage in LuaMetaTeX where we don't want to exit on an 
+    error so we intercept it. 
+
+*/
+
+#include "cerf.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+double dvoigt(double x, double sigma, double gamma, double v0)
+{
+    return voigt(x, sigma, gamma)/v0 - .5;
+}
+
+double voigt_hwhm(double sigma, double gamma, int *error)
+{
+    *error = 0;
+    if (sigma == 0 && gamma == 0) {
+        return 0;
+    } else if (isnan(sigma) || isnan(gamma)) {        
+        *error = 1; 
+        return 0; // return NAN;
+    } else {
+        // start from an excellent approximation [Olivero & Longbothum, J Quant Spec Rad Transf 1977]:
+        const double eps = 1e-14;
+        const double hwhm0 = .5*(1.06868*gamma+sqrt(0.86743*gamma*gamma+4*2*log(2)*sigma*sigma));
+        const double del = eps*hwhm0;
+        double ret = hwhm0;
+        const double v0 = voigt(0, sigma, gamma);
+        for (int i=0; i<300; ++i) {
+            double val = dvoigt(ret, sigma, gamma, v0);
+            if (fabs(val) < 1e-15) {
+                return ret;
+            } else {
+                double step = -del/(dvoigt(ret+del, sigma, gamma, v0)/val-1);
+                double nxt = ret + step;
+                if (nxt < ret/3) {
+                    *error = 2; // fprintf(stderr, "voigt_fwhm terminated because of huge deviation from 1st approx\n");
+                    nxt = ret/3;
+                } else if (nxt > 2*ret) {
+                    *error = 2; // fprintf(stderr, "voigt_fwhm terminated because of huge deviation from 1st approx\n");
+                    nxt = 2*ret;
+                }
+                if (fabs(ret-nxt) < del) {
+                    return nxt;
+                } else { 
+                    ret = nxt;
+                }
+            }
+        }
+        *error = 3; // fprintf(stderr, "voigt_fwhm failed: Newton's iteration did not converge with sigma = %f and gamma = %f\n", sigma, gamma);  exit(-1);
+        return 0;
+    }
+}