Commit dd24d6c6 authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

Improve voigt_hwhm, and provide test

  - Added widthtest to test voigt_hwhm for numerous ratios gamma/sigma
  - Use assertion in voigt_hwhm to check for impossible situations
      (to avoid fprintf and exit, as suggested by Etham A Merrit)
parent 07aaaa85
Pipeline #38857 passed with stage
in 15 seconds
== Revision history of libcerf, maintained by Joachim Wuttke ==
libcerf-1.15:
- Added widthtest to test voigt_hwhm for numerous ratios gamma/sigma
- Use assertion in voigt_hwhm to check for impossible situations
(to avoid fprintf and exit, as suggested by Etham A Merrit)
libcerf-1.14, released 19oct20:
- Simplified test code
- Homepage moved to https://jugit.fz-juelich.de/mlz/libcerf, 17mar19
......
......@@ -42,9 +42,8 @@
*/
#include "cerf.h"
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
double dvoigt(double x, double sigma, double gamma, double v0)
{
......@@ -63,24 +62,17 @@ double voigt_hwhm(double sigma, double gamma)
const double del = eps*hwhm0;
double ret = hwhm0;
const double v0 = voigt(0, sigma, gamma);
for(int i=0; i<300; ++i) {
for(int i=0; i<30; ++i) {
double val = dvoigt(ret, sigma, gamma, v0);
if (fabs(val)<1e-15)
return ret;
double step = -del/(dvoigt(ret+del, sigma, gamma, v0)/val-1);
double nxt = ret + step;
if (nxt<ret/3) {
fprintf(stderr, "voigt_fwhm terminated because of huge deviation from 1st approx\n");
nxt = ret/3;
} else if (nxt>2*ret) {
fprintf(stderr, "voigt_fwhm terminated because of huge deviation from 1st approx\n");
nxt = 2*ret;
}
assert(nxt>ret/3);
assert(nxt<2*ret);
if (fabs(ret-nxt)<del)
return nxt;
ret = nxt;
}
fprintf(stderr, "voigt_fwhm failed: Newton's iteration did not converge"
" with sigma = %f and gamma = %f\n", sigma, gamma);
exit(-1);
assert(0);
}
/* Library libcerf:
* Compute complex error functions, based on a new implementation of
* Faddeeva's w_of_z. Also provide Dawson and Voigt functions.
*
* File widthtest.c:
* Test function voigt_hwhm
*
* Copyright:
* (C) 2021 Forschungszentrum Jülich GmbH
*
* Licence:
* ../LICENSE
*
* Authors:
* Joachim Wuttke, Forschungszentrum Jülich, 2021
*
* Website:
* http://apps.jcns.fz-juelich.de/libcerf
*
* Revision history:
* ../CHANGELOG
*
* More information:
* man 3 voigt_hwhm
*/
#include "cerf.h"
#include "testtool.h"
#include <math.h>
#include <stdio.h>
// excellent approximation [Olivero & Longbothum, 1977], used as starting value in voigt_hwhm
double hwhm0(double sigma, double gamma)
{
return .5*(1.06868*gamma+sqrt(0.86743*gamma*gamma+4*2*log(2)*sigma*sigma));
}
int main()
{
result_t result;
const int N = 100000;
for (int i = 0; i<N; ++i ) {
const double gamma = pow(10., 180*(i-N/2)/N);
RTEST(result, 1e-2, voigt_hwhm(1., gamma), hwhm0(1.,gamma));
}
printf("%i/%i tests failed", result.failed, result.total);
return result.failed ? 1 : 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment