summaryrefslogtreecommitdiff
path: root/src/silx/math/fit/functions/src/funs.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/silx/math/fit/functions/src/funs.c')
-rw-r--r--src/silx/math/fit/functions/src/funs.c1357
1 files changed, 1357 insertions, 0 deletions
diff --git a/src/silx/math/fit/functions/src/funs.c b/src/silx/math/fit/functions/src/funs.c
new file mode 100644
index 0000000..4b41fce
--- /dev/null
+++ b/src/silx/math/fit/functions/src/funs.c
@@ -0,0 +1,1357 @@
+#/*##########################################################################
+# Copyright (c) 2004-2016 European Synchrotron Radiation Facility
+#
+# 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.
+#
+#############################################################################*/
+/*
+ This file provides fit functions.
+
+ It is adapted from PyMca source file "SpecFitFuns.c". The main difference
+ with the original code is that this code does not handle the python
+ wrapping, which is done elsewhere using cython.
+
+ Authors: V.A. Sole, P. Knobel
+ License: MIT
+ Last modified: 17/06/2016
+*/
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "functions.h"
+
+#ifndef M_PI
+#define M_PI 3.1415926535
+#endif
+
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+#define MAX(x, y) (((x) > (y)) ? (x) : (y))
+
+#if defined(_WIN32)
+#define erf myerf
+#define erfc myerfc
+#endif
+
+#define LOG2 0.69314718055994529
+
+
+int test_params(int len_params,
+ int len_params_one_function,
+ char* fun_name,
+ char* param_names)
+{
+ if (len_params % len_params_one_function) {
+ printf("[%s]Error: Number of parameters must be a multiple of %d.",
+ fun_name, len_params_one_function);
+ printf("\nParameters expected for %s: %s\n",
+ fun_name, param_names);
+ return(1);
+ }
+ if (len_params == 0) {
+ printf("[%s]Error: No parameters specified.", fun_name);
+ printf("\nParameters expected for %s: %s\n",
+ fun_name, param_names);
+ return(1);
+ }
+ return(0);
+}
+
+/* Complementary error function for a single value*/
+double myerfc(double x)
+{
+ double z;
+ double t;
+ double r;
+
+ z=fabs(x);
+ t=1.0/(1.0+0.5*z);
+ r=t * exp(-z * z - 1.26551223 + t * (1.00002368 + t * (0.3740916 +
+ t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 +
+ t * (1.48851587 + t * (-0.82215223+t*0.17087277)))))))));
+ if (x<0)
+ r=2.0-r;
+ return (r);
+}
+
+/* Gauss error function for a single value*/
+double myerf(double x)
+{
+ return (1.0 - myerfc(x));
+}
+
+/* Gauss error function for an array
+ y[i]=erf(x[i])
+ returns status code 0
+*/
+int erf_array(double* x, int len_x, double* y)
+{
+ int j;
+ for (j=0; j<len_x; j++) {
+ y[j] = erf(x[j]);
+ }
+ return(0);
+}
+
+/* Complementary error function for an array
+ y[i]=erfc(x[i])
+ returns status code 0*/
+int erfc_array(double* x, int len_x, double* y)
+{
+ int j;
+ for (j=0; j<len_x; j++) {
+ y[j] = erfc(x[j]);
+ }
+ return(0);
+}
+
+/* Use lookup table for fast exp computation */
+double fastexp(double x)
+{
+ int expindex;
+ static double EXP[5000] = {0.0};
+ int i;
+
+/*initialize */
+ if (EXP[0] < 1){
+ for (i=0;i<5000;i++){
+ EXP[i] = exp(-0.01 * i);
+ }
+ }
+/*calculate*/
+ if (x < 0){
+ x = -x;
+ if (x < 50){
+ expindex = (int) (x * 100);
+ return EXP[expindex]*(1.0 - (x - 0.01 * expindex)) ;
+ }else if (x < 100) {
+ expindex = (int) (x * 10);
+ return pow(EXP[expindex]*(1.0 - (x - 0.1 * expindex)),10) ;
+ }else if (x < 1000){
+ expindex = (int) x;
+ return pow(EXP[expindex]*(1.0 - (x - expindex)),20) ;
+ }else if (x < 10000){
+ expindex = (int) (x * 0.1);
+ return pow(EXP[expindex]*(1.0 - (x - 10.0 * expindex)),30) ;
+ }else{
+ return 0;
+ }
+ }else{
+ if (x < 50){
+ expindex = (int) (x * 100);
+ return 1.0/EXP[expindex]*(1.0 - (x - 0.01 * expindex)) ;
+ }else if (x < 100) {
+ expindex = (int) (x * 10);
+ return pow(EXP[expindex]*(1.0 - (x - 0.1 * expindex)),-10) ;
+ }else{
+ return exp(x);
+ }
+ }
+}
+
+
+/* sum_gauss
+ Sum of gaussian functions, defined by (height, centroid, fwhm)
+
+ *height* is the peak amplitude
+ *centroid* is the peak x-coordinate
+ *fwhm* is the full-width at half maximum
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pvoigt: Array of gaussian parameters:
+ (height1, centroid1, fwhm1, height2, centroid2, fwhm2,...)
+ - len_pgauss: Number of elements in the pgauss array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_gauss(double* x, int len_x, double* pgauss, int len_pgauss, double* y)
+{
+ int i, j;
+ double dhelp, inv_two_sqrt_two_log2, sigma;
+ double fwhm, centroid, height;
+
+ if (test_params(len_pgauss, 3, "sum_gauss", "height, centroid, fwhm")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pgauss/3; i++) {
+ height = pgauss[3*i];
+ centroid = pgauss[3*i+1];
+ fwhm = pgauss[3*i+2];
+
+ sigma = fwhm * inv_two_sqrt_two_log2;
+
+ for (j=0; j<len_x; j++) {
+ dhelp = (x[j] - centroid) / sigma;
+ if (dhelp <= 20) {
+ y[j] += height * exp (-0.5 * dhelp * dhelp);
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_agauss
+ Sum of gaussian functions defined by (area, centroid, fwhm)
+
+ *area* is the area underneath the peak
+ *centroid* is the peak x-coordinate
+ *fwhm* is the full-width at half maximum
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pgauss: Array of gaussian parameters:
+ (area1, centroid1, fwhm1, area2, centroid2, fwhm2,...)
+ - len_pgauss: Number of elements in the pgauss array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_agauss(double* x, int len_x, double* pgauss, int len_pgauss, double* y)
+{
+ int i, j;
+ double dhelp, height, sqrt2PI, sigma, inv_two_sqrt_two_log2;
+ double fwhm, centroid, area;
+
+ if (test_params(len_pgauss, 3, "sum_agauss", "area, centroid, fwhm")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+ sqrt2PI = sqrt(2.0*M_PI);
+
+ for (i=0; i<len_pgauss/3; i++) {
+ area = pgauss[3*i];
+ centroid = pgauss[3*i+1];
+ fwhm = pgauss[3*i+2];
+
+ sigma = fwhm * inv_two_sqrt_two_log2;
+ height = area / (sigma * sqrt2PI);
+
+ for (j=0; j<len_x; j++) {
+ dhelp = (x[j] - centroid)/sigma;
+ if (dhelp <= 35) {
+ y[j] += height * exp (-0.5 * dhelp * dhelp);
+ }
+ }
+ }
+ return(0);
+}
+
+
+/* sum_fastagauss
+ Sum of gaussian functions defined by (area, centroid, fwhm).
+ This implementation uses a lookup table of precalculated exp values
+ and a limited development (exp(-x) = 1 - x for small values of x)
+
+ *area* is the area underneath the peak
+ *centroid* is the peak x-coordinate
+ *fwhm* is the full-width at half maximum
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pgauss: Array of gaussian parameters:
+ (area1, centroid1, fwhm1, area2, centroid2, fwhm2,...)
+ - len_pgauss: Number of elements in the pgauss array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+
+int sum_fastagauss(double* x, int len_x, double* pgauss, int len_pgauss, double* y)
+{
+ int i, j, expindex;
+ double dhelp, height, sqrt2PI, sigma, inv_two_sqrt_two_log2;
+ double fwhm, centroid, area;
+ static double EXP[5000];
+
+ if (test_params(len_pgauss, 3, "sum_fastagauss", "area, centroid, fwhm")) {
+ return(1);
+ }
+
+ if (EXP[0] < 1){
+ for (i=0; i<5000; i++){
+ EXP[i] = exp(-0.01 * i);
+ }
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+ sqrt2PI = sqrt(2.0*M_PI);
+
+ for (i=0; i<len_pgauss/3; i++) {
+ area = pgauss[3*i];
+ centroid = pgauss[3*i+1];
+ fwhm = pgauss[3*i+2];
+
+ sigma = fwhm * inv_two_sqrt_two_log2;
+ height = area / (sigma * sqrt2PI);
+
+ for (j=0; j<len_x; j++) {
+ dhelp = (x[j] - centroid)/sigma;
+ if (dhelp <= 15){
+ dhelp = 0.5 * dhelp * dhelp;
+ if (dhelp < 50){
+ expindex = (int) (dhelp * 100);
+ y[j] += height * EXP[expindex] * (1.0 - (dhelp - 0.01 * expindex));
+ }
+ else if (dhelp < 100) {
+ expindex = (int) (dhelp * 10);
+ y[j] += height * pow(EXP[expindex] * (1.0 - (dhelp - 0.1 * expindex)), 10);
+ }
+ else if (dhelp < 1000){
+ expindex = (int) (dhelp);
+ y[j] += height * pow(EXP[expindex] * (1.0 - (dhelp - expindex)), 20);
+ }
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_splitgauss
+ Sum of split gaussian functions, defined by (height, centroid, fwhm1, fwhm2)
+
+ *height* is the peak amplitude
+ *centroid* is the peak x-coordinate
+ *fwhm1* is the full-width at half maximum of the left half of the curve (x < centroid)
+ *fwhm1* is the full-width at half maximum of the right half of the curve (x > centroid)
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pgauss: Array of gaussian parameters:
+ (height1, centroid1, fwhm11, fwhm21, height2, centroid2, fwhm12, fwhm22,...)
+ - len_pgauss: Number of elements in the pgauss array. Must be
+ a multiple of 4.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_splitgauss(double* x, int len_x, double* pgauss, int len_pgauss, double* y)
+{
+ int i, j;
+ double dhelp, inv_two_sqrt_two_log2, sigma1, sigma2;
+ double fwhm1, fwhm2, centroid, height;
+
+ if (test_params(len_pgauss, 4, "sum_splitgauss", "height, centroid, fwhm1, fwhm2")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pgauss/4; i++) {
+ height = pgauss[4*i];
+ centroid = pgauss[4*i+1];
+ fwhm1 = pgauss[4*i+2];
+ fwhm2 = pgauss[4*i+3];
+
+ sigma1 = fwhm1 * inv_two_sqrt_two_log2;
+ sigma2 = fwhm2 * inv_two_sqrt_two_log2;
+
+ for (j=0; j<len_x; j++) {
+ dhelp = (x[j] - centroid);
+ if (dhelp > 0) {
+ /* Use fwhm2 when x > centroid */
+ dhelp = dhelp / sigma2;
+ }
+ else {
+ /* Use fwhm1 when x < centroid */
+ dhelp = dhelp / sigma1;
+ }
+
+ if (dhelp <= 20) {
+ y[j] += height * exp (-0.5 * dhelp * dhelp);
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_apvoigt
+ Sum of pseudo-Voigt functions, defined by (area, centroid, fwhm, eta).
+
+ The pseudo-Voigt profile PV(x) is an approximation of the Voigt profile
+ using a linear combination of a Gaussian curve G(x) and a Lorentzian curve
+ L(x) instead of their convolution.
+
+ *area* is the area underneath both G(x) and L(x)
+ *centroid* is the peak x-coordinate for both functions
+ *fwhm* is the full-width at half maximum of both functions
+ *eta* is the Lorentzian fraction: PV(x) = eta * L(x) + (1 - eta) * G(x)
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pvoigt: Array of Voigt function parameters:
+ (area1, centroid1, fwhm1, eta1, area2, centroid2, fwhm2, eta2,...)
+ - len_voigt: Number of elements in the pvoigt array. Must be
+ a multiple of 4.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_apvoigt(double* x, int len_x, double* pvoigt, int len_pvoigt, double* y)
+{
+ int i, j;
+ double dhelp, inv_two_sqrt_two_log2, sqrt2PI, sigma, height;
+ double area, centroid, fwhm, eta;
+
+ if (test_params(len_pvoigt, 4, "sum_apvoigt", "area, centroid, fwhm, eta")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+ sqrt2PI = sqrt(2.0*M_PI);
+
+
+ for (i=0; i<len_pvoigt/4; i++) {
+ area = pvoigt[4*i];
+ centroid = pvoigt[4*i+1];
+ fwhm = pvoigt[4*i+2];
+ eta = pvoigt[4*i+3];
+
+ sigma = fwhm * inv_two_sqrt_two_log2;
+ height = area / (sigma * sqrt2PI);
+
+ for (j=0; j<len_x; j++) {
+ /* Lorentzian term */
+ dhelp = (x[j] - centroid) / (0.5 * fwhm);
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += eta * (area / (0.5 * M_PI * fwhm * dhelp));
+
+ /* Gaussian term */
+ dhelp = (x[j] - centroid) / sigma;
+ if (dhelp <= 35) {
+ y[j] += (1.0 - eta) * height * exp (-0.5 * dhelp * dhelp);
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_pvoigt
+ Sum of pseudo-Voigt functions, defined by (height, centroid, fwhm, eta).
+
+ The pseudo-Voigt profile PV(x) is an approximation of the Voigt profile
+ using a linear combination of a Gaussian curve G(x) and a Lorentzian curve
+ L(x) instead of their convolution.
+
+ *height* is the peak amplitude of G(x) and L(x)
+ *centroid* is the peak x-coordinate for both functions
+ *fwhm* is the full-width at half maximum of both functions
+ *eta* is the Lorentzian fraction: PV(x) = eta * L(x) + (1 - eta) * G(x)
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pvoigt: Array of Voigt function parameters:
+ (height1, centroid1, fwhm1, eta1, height2, centroid2, fwhm2, eta2,...)
+ - len_voigt: Number of elements in the pvoigt array. Must be
+ a multiple of 4.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_pvoigt(double* x, int len_x, double* pvoigt, int len_pvoigt, double* y)
+{
+ int i, j;
+ double dhelp, inv_two_sqrt_two_log2, sigma;
+ double height, centroid, fwhm, eta;
+
+ if (test_params(len_pvoigt, 4, "sum_pvoigt", "height, centroid, fwhm, eta")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pvoigt/4; i++) {
+ height = pvoigt[4*i];
+ centroid = pvoigt[4*i+1];
+ fwhm = pvoigt[4*i+2];
+ eta = pvoigt[4*i+3];
+
+ sigma = fwhm * inv_two_sqrt_two_log2;
+
+ for (j=0; j<len_x; j++) {
+ /* Lorentzian term */
+ dhelp = (x[j] - centroid) / (0.5 * fwhm);
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += eta * height / dhelp;
+
+ /* Gaussian term */
+ dhelp = (x[j] - centroid) / sigma;
+ if (dhelp <= 35) {
+ y[j] += (1.0 - eta) * height * exp (-0.5 * dhelp * dhelp);
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_splitpvoigt
+ Sum of split pseudo-Voigt functions, defined by
+ (height, centroid, fwhm1, fwhm2, eta).
+
+ The pseudo-Voigt profile PV(x) is an approximation of the Voigt profile
+ using a linear combination of a Gaussian curve G(x) and a Lorentzian curve
+ L(x) instead of their convolution.
+
+ *height* is the peak amplitude of G(x) and L(x)
+ *centroid* is the peak x-coordinate for both functions
+ *fwhm1* is the full-width at half maximum of both functions for x < centroid
+ *fwhm2* is the full-width at half maximum of both functions for x > centroid
+ *eta* is the Lorentzian fraction: PV(x) = eta * L(x) + (1 - eta) * G(x)
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pvoigt: Array of Voigt function parameters:
+ (height1, centroid1, fwhm11, fwhm21, eta1, ...)
+ - len_voigt: Number of elements in the pvoigt array. Must be
+ a multiple of 5.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_splitpvoigt(double* x, int len_x, double* pvoigt, int len_pvoigt, double* y)
+{
+ int i, j;
+ double dhelp, inv_two_sqrt_two_log2, x_minus_centroid, sigma1, sigma2;
+ double height, centroid, fwhm1, fwhm2, eta;
+
+ if (test_params(len_pvoigt, 5, "sum_splitpvoigt", "height, centroid, fwhm1, fwhm2, eta")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pvoigt/5; i++) {
+ height = pvoigt[5*i];
+ centroid = pvoigt[5*i+1];
+ fwhm1 = pvoigt[5*i+2];
+ fwhm2 = pvoigt[5*i+3];
+ eta = pvoigt[5*i+4];
+
+ sigma1 = fwhm1 * inv_two_sqrt_two_log2;
+ sigma2 = fwhm2 * inv_two_sqrt_two_log2;
+
+ for (j=0; j<len_x; j++) {
+ x_minus_centroid = (x[j] - centroid);
+
+ /* Use fwhm2 when x > centroid */
+ if (x_minus_centroid > 0) {
+ /* Lorentzian term */
+ dhelp = x_minus_centroid / (0.5 * fwhm2);
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += eta * height / dhelp;
+
+ /* Gaussian term */
+ dhelp = x_minus_centroid / sigma2;
+ if (dhelp <= 35) {
+ y[j] += (1.0 - eta) * height * exp (-0.5 * dhelp * dhelp);
+ }
+ }
+ /* Use fwhm1 when x < centroid */
+ else {
+ /* Lorentzian term */
+ dhelp = x_minus_centroid / (0.5 * fwhm1);
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += eta * height / dhelp;
+
+ /* Gaussian term */
+ dhelp = x_minus_centroid / sigma1;
+ if (dhelp <= 35) {
+ y[j] += (1.0 - eta) * height * exp (-0.5 * dhelp * dhelp);
+ }
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_splitpvoigt2
+ Sum of split pseudo-Voigt functions, defined by
+ (height, centroid, fwhm1, fwhm2, eta1, eta2).
+
+ The pseudo-Voigt profile PV(x) is an approximation of the Voigt profile
+ using a linear combination of a Gaussian curve G(x) and a Lorentzian curve
+ L(x) instead of their convolution.
+
+ *height* is the peak amplitude of G(x) and L(x)
+ *centroid* is the peak x-coordinate for both functions
+ *fwhm1* is the full-width at half maximum of both functions for x < centroid
+ *fwhm2* is the full-width at half maximum of both functions for x > centroid
+ *eta1* is the Lorentzian fraction for x < centroid
+ *eta2* is the Lorentzian fraction for x > centroid
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the gaussians are calculated.
+ - len_x: Number of elements in the x array.
+ - pvoigt: Array of Voigt function parameters:
+ (height1, centroid1, fwhm11, fwhm21, eta11, eta21, ...)
+ - len_voigt: Number of elements in the pvoigt array. Must be
+ a multiple of 6.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+*/
+int sum_splitpvoigt2(double* x, int len_x, double* pvoigt, int len_pvoigt, double* y)
+{
+ int i, j;
+ double dhelp, x_minus_centroid, inv_two_sqrt_two_log2, sigma1, sigma2;
+ double height, centroid, fwhm1, fwhm2, eta1, eta2;
+
+ if (test_params(len_pvoigt, 6, "sum_splitpvoigt2", "height, centroid, fwhm1, fwhm2, eta1, eta2")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ inv_two_sqrt_two_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pvoigt/6; i++) {
+ height = pvoigt[6*i];
+ centroid = pvoigt[6*i+1];
+ fwhm1 = pvoigt[6*i+2];
+ fwhm2 = pvoigt[6*i+3];
+ eta1 = pvoigt[6*i+4];
+ eta2 = pvoigt[6*i+5];
+
+ sigma1 = fwhm1 * inv_two_sqrt_two_log2;
+ sigma2 = fwhm2 * inv_two_sqrt_two_log2;
+
+ for (j=0; j<len_x; j++) {
+ x_minus_centroid = (x[j] - centroid);
+
+ /* Use fwhm2 and eta2 when x > centroid */
+ if (x_minus_centroid > 0) {
+ /* Lorentzian term */
+ dhelp = (2.0 * x_minus_centroid) / fwhm2;
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += eta2 * height / dhelp;
+
+ /* Gaussian term */
+ dhelp = x_minus_centroid / sigma2;
+ if (dhelp <= 35) {
+ dhelp = exp(-0.5 * dhelp * dhelp);
+ y[j] += (1 - eta2) * height * dhelp;
+ }
+ }
+ /* Use fwhm1 and eta1 when x < centroid */
+ else {
+ /* Lorentzian term */
+ dhelp = (2.0 * x_minus_centroid) / fwhm1;
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += eta1 * height / dhelp;
+
+ /* Gaussian term */
+ dhelp = x_minus_centroid / sigma1;
+ if (dhelp <= 35) {
+ dhelp = exp(-0.5 * dhelp * dhelp);
+ y[j] += (1 - eta1) * height * dhelp;
+ }
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_lorentz
+ Sum of Lorentz functions, defined by (height, centroid, fwhm).
+
+ *height* is the peak amplitude
+ *centroid* is the peak's x-coordinate
+ *fwhm* is the full-width at half maximum
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the Lorentzians are calculated.
+ - len_x: Number of elements in the x array.
+ - plorentz: Array of lorentz function parameters:
+ (height1, centroid1, fwhm1, ...)
+ - len_lorentz: Number of elements in the plorentz array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_lorentz(double* x, int len_x, double* plorentz, int len_plorentz, double* y)
+{
+ int i, j;
+ double dhelp;
+ double height, centroid, fwhm;
+
+ if (test_params(len_plorentz, 3, "sum_lorentz", "height, centroid, fwhm")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ for (i=0; i<len_plorentz/3; i++) {
+ height = plorentz[3*i];
+ centroid = plorentz[3*i+1];
+ fwhm = plorentz[3*i+2];
+
+ for (j=0; j<len_x; j++) {
+ dhelp = (x[j] - centroid) / (0.5 * fwhm);
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += height / dhelp;
+ }
+ }
+ return(0);
+}
+
+
+/* sum_alorentz
+ Sum of Lorentz functions, defined by (area, centroid, fwhm).
+
+ *area* is the area underneath the peak
+ *centroid* is the peak's x-coordinate
+ *fwhm* is the full-width at half maximum
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the Lorentzians are calculated.
+ - len_x: Number of elements in the x array.
+ - plorentz: Array of lorentz function parameters:
+ (area1, centroid1, fwhm1, ...)
+ - len_lorentz: Number of elements in the plorentz array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_alorentz(double* x, int len_x, double* plorentz, int len_plorentz, double* y)
+{
+ int i, j;
+ double dhelp;
+ double area, centroid, fwhm;
+
+ if (test_params(len_plorentz, 3, "sum_alorentz", "area, centroid, fwhm")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ for (i=0; i<len_plorentz/3; i++) {
+ area = plorentz[3*i];
+ centroid = plorentz[3*i+1];
+ fwhm = plorentz[3*i+2];
+
+ for (j=0; j<len_x; j++) {
+ dhelp = (x[j] - centroid) / (0.5 * fwhm);
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += area / (0.5 * M_PI * fwhm * dhelp);
+ }
+ }
+ return(0);
+}
+
+
+/* sum_splitlorentz
+ Sum of Lorentz functions, defined by (height, centroid, fwhm1, fwhm2).
+
+ *height* is the peak amplitude
+ *centroid* is the peak's x-coordinate
+ *fwhm1* is the full-width at half maximum for x < centroid
+ *fwhm2* is the full-width at half maximum for x > centroid
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the Lorentzians are calculated.
+ - len_x: Number of elements in the x array.
+ - plorentz: Array of lorentz function parameters:
+ (height1, centroid1, fwhm11, fwhm21 ...)
+ - len_lorentz: Number of elements in the plorentz array. Must be
+ a multiple of 4.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_splitlorentz(double* x, int len_x, double* plorentz, int len_plorentz, double* y)
+{
+ int i, j;
+ double dhelp;
+ double height, centroid, fwhm1, fwhm2;
+
+ if (test_params(len_plorentz, 4, "sum_splitlorentz", "height, centroid, fwhm1, fwhm2")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ for (i=0; i<len_plorentz/4; i++) {
+ height = plorentz[4*i];
+ centroid = plorentz[4*i+1];
+ fwhm1 = plorentz[4*i+2];
+ fwhm2 = plorentz[4*i+3];
+
+ for (j=0; j<len_x; j++) {
+ dhelp = (x[j] - centroid);
+ if (dhelp>0) {
+ dhelp = dhelp / (0.5 * fwhm2);
+ }
+ else {
+ dhelp = dhelp / (0.5 * fwhm1);
+ }
+ dhelp = 1.0 + (dhelp * dhelp);
+ y[j] += height / dhelp;
+ }
+ }
+ return(0);
+}
+
+/* sum_stepdown
+ Sum of stepdown functions, defined by (height, centroid, fwhm).
+
+ *height* is the step amplitude
+ *centroid* is the step's x-coordinate
+ *fwhm* is the full-width at half maximum of the derivative
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the stepdown functions are calculated.
+ - len_x: Number of elements in the x array.
+ - pdstep: Array of downstpe function parameters:
+ (height1, centroid1, fwhm1, ...)
+ - len_pdstep: Number of elements in the pdstep array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_stepdown(double* x, int len_x, double* pdstep, int len_pdstep, double* y)
+{
+ int i, j;
+ double dhelp, sqrt2_inv_2_sqrt_two_log2 ;
+ double height, centroid, fwhm;
+
+ if (test_params(len_pdstep, 3, "sum_stepdown", "height, centroid, fwhm")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ sqrt2_inv_2_sqrt_two_log2 = sqrt(2.0) / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pdstep/3; i++) {
+ height = pdstep[3*i];
+ centroid = pdstep[3*i+1];
+ fwhm = pdstep[3*i+2];
+
+ for (j=0; j<len_x; j++) {
+ dhelp = fwhm * sqrt2_inv_2_sqrt_two_log2;
+ dhelp = (x[j] - centroid) / dhelp;
+ y[j] += height * 0.5 * erfc(dhelp);
+ }
+ }
+ return(0);
+}
+
+/* sum_stepup
+ Sum of stepup functions, defined by (height, centroid, fwhm).
+
+ *height* is the step amplitude
+ *centroid* is the step's x-coordinate
+ *fwhm* is the full-width at half maximum of the derivative
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the stepup functions are calculated.
+ - len_x: Number of elements in the x array.
+ - pustep: Array of stepdown function parameters:
+ (height1, centroid1, fwhm1, ...)
+ - len_pustep: Number of elements in the pustep array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_stepup(double* x, int len_x, double* pustep, int len_pustep, double* y)
+{
+ int i, j;
+ double dhelp, sqrt2_inv_2_sqrt_two_log2 ;
+ double height, centroid, fwhm;
+
+ if (test_params(len_pustep, 3, "sum_stepup", "height, centroid, fwhm")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ sqrt2_inv_2_sqrt_two_log2 = sqrt(2.0) / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pustep/3; i++) {
+ height = pustep[3*i];
+ centroid = pustep[3*i+1];
+ fwhm = pustep[3*i+2];
+
+ for (j=0; j<len_x; j++) {
+ dhelp = fwhm * sqrt2_inv_2_sqrt_two_log2;
+ dhelp = (x[j] - centroid) / dhelp;
+ y[j] += height * 0.5 * (1.0 + erf(dhelp));
+ }
+ }
+ return(0);
+}
+
+
+/* sum_slit
+ Sum of slit functions, defined by (height, position, fwhm, beamfwhm).
+
+ *height* is the slit height
+ *position* is the slit's center x-coordinate
+ *fwhm* is the full-width at half maximum of the slit
+ *beamfwhm* is the full-width at half maximum of derivative's peaks
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the slit functions are calculated.
+ - len_x: Number of elements in the x array.
+ - pslit: Array of slit function parameters:
+ (height1, centroid1, fwhm1, beamfwhm1 ...)
+ - len_pslit: Number of elements in the pslit array. Must be
+ a multiple of 3.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_slit(double* x, int len_x, double* pslit, int len_pslit, double* y)
+{
+ int i, j;
+ double dhelp, dhelp1, dhelp2, sqrt2_inv_2_sqrt_two_log2, centroid1, centroid2;
+ double height, position, fwhm, beamfwhm;
+
+ if (test_params(len_pslit, 4, "sum_slit", "height, centroid, fwhm, beamfwhm")) {
+ return(1);
+ }
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ sqrt2_inv_2_sqrt_two_log2 = sqrt(2.0) / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_pslit/4; i++) {
+ height = pslit[4*i];
+ position = pslit[4*i+1];
+ fwhm = pslit[4*i+2];
+ beamfwhm = pslit[4*i+3];
+
+ centroid1 = position - 0.5 * fwhm;
+ centroid2 = position + 0.5 * fwhm;
+
+ for (j=0; j<len_x; j++) {
+ dhelp = beamfwhm * sqrt2_inv_2_sqrt_two_log2;
+ dhelp1 = (x[j] - centroid1) / dhelp;
+ dhelp2 = (x[j] - centroid2) / dhelp;
+ y[j] += height * 0.25 * (1.0 + erf(dhelp1)) * erfc(dhelp2);
+ }
+ }
+ return(0);
+}
+
+
+/* sum_ahypermet
+ Sum of hypermet functions, defined by
+ (area, position, fwhm, st_area_r, st_slope_r, lt_area_r, lt_slope_r, step_height_r).
+
+ - *area* is the area underneath the gaussian peak
+ - *position* is the center of the various peaks and the position of
+ the step down
+ - *fwhm* is the full-width at half maximum of the terms
+ - *st_area_r* is factor between the gaussian area and the area of the
+ short tail term
+ - *st_slope_r* is a parameter related to the slope of the short tail
+ in the low ``x`` values (the lower, the steeper)
+ - *lt_area_r* is factor between the gaussian area and the area of the
+ long tail term
+ - *lt_slope_r* is a parameter related to the slope of the long tail
+ in the low ``x`` values (the lower, the steeper)
+ - *step_height_r* is the factor between the height of the step down
+ and the gaussian height
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the functions are calculated.
+ - len_x: Number of elements in the x array.
+ - phypermet: Array of hypermet function parameters:
+ *(area1, position1, fwhm1, st_area_r1, st_slope_r1, lt_area_r1,
+ lt_slope_r1, step_height_r1, ...)*
+ - len_phypermet: Number of elements in the phypermet array. Must be
+ a multiple of 8.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+ - tail_flags: sum of binary flags to activate the various terms of the
+ function:
+
+ - 1 (b0001): Gaussian term
+ - 2 (b0010): st term
+ - 4 (b0100): lt term
+ - 8 (b1000): step term
+
+ E.g., to activate all termsof the hypermet, use ``tail_flags = 1 + 2 + 4 + 8 = 15``
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_ahypermet(double* x, int len_x, double* phypermet, int len_phypermet, double* y, int tail_flags)
+{
+ int i, j;
+ int g_term_flag, st_term_flag, lt_term_flag, step_term_flag;
+ double c1, c2, sigma, height, sigma_sqrt2, sqrt2PI, inv_2_sqrt_2_log2, x_minus_position, epsilon;
+ double area, position, fwhm, st_area_r, st_slope_r, lt_area_r, lt_slope_r, step_height_r;
+
+ if (test_params(len_phypermet, 8, "sum_hypermet",
+ "height, centroid, fwhm, st_area_r, st_slope_r, lt_area_r, lt_slope_r, step_height_r")) {
+ return(1);
+ }
+
+ g_term_flag = tail_flags & 1;
+ st_term_flag = (tail_flags>>1) & 1;
+ lt_term_flag = (tail_flags>>2) & 1;
+ step_term_flag = (tail_flags>>3) & 1;
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ /* define epsilon to compare floating point values with 0. */
+ epsilon = 0.00000000001;
+
+ sqrt2PI= sqrt(2.0 * M_PI);
+ inv_2_sqrt_2_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_phypermet/8; i++) {
+ area = phypermet[8*i];
+ position = phypermet[8*i+1];
+ fwhm = phypermet[8*i+2];
+ st_area_r = phypermet[8*i+3];
+ st_slope_r = phypermet[8*i+4];
+ lt_area_r = phypermet[8*i+5];
+ lt_slope_r = phypermet[8*i+6];
+ step_height_r = phypermet[8*i+7];
+
+ sigma = fwhm * inv_2_sqrt_2_log2;
+ height = area / (sigma * sqrt2PI);
+
+ /* Prevent division by 0 */
+ if (sigma == 0) {
+ printf("fwhm must not be equal to 0");
+ return(1);
+ }
+ sigma_sqrt2 = sigma * 1.4142135623730950488;
+
+ for (j=0; j<len_x; j++) {
+ x_minus_position = x[j] - position;
+ c2 = (0.5 * x_minus_position * x_minus_position) / (sigma * sigma);
+ /* gaussian term */
+ if (g_term_flag) {
+ y[j] += exp(-c2) * height;
+ }
+
+ /* st term */
+ if (st_term_flag) {
+ if (fabs(st_slope_r) > epsilon) {
+ c1 = st_area_r * 0.5 * \
+ erfc((x_minus_position/sigma_sqrt2) + 0.5 * sigma_sqrt2 / st_slope_r);
+ y[j] += ((area * c1) / st_slope_r) * \
+ exp(0.5 * (sigma / st_slope_r) * (sigma / st_slope_r) + \
+ (x_minus_position / st_slope_r));
+ }
+ }
+
+ /* lt term */
+ if (lt_term_flag) {
+ if (fabs(lt_slope_r) > epsilon) {
+ c1 = lt_area_r * \
+ 0.5 * erfc((x_minus_position/sigma_sqrt2) + 0.5 * sigma_sqrt2 / lt_slope_r);
+ y[j] += ((area * c1) / lt_slope_r) * \
+ exp(0.5 * (sigma / lt_slope_r) * (sigma / lt_slope_r) + \
+ (x_minus_position / lt_slope_r));
+ }
+ }
+
+ /* step term flag */
+ if (step_term_flag) {
+ y[j] += step_height_r * (area / (sigma * sqrt2PI)) * \
+ 0.5 * erfc(x_minus_position / sigma_sqrt2);
+ }
+ }
+ }
+ return(0);
+}
+
+/* sum_fastahypermet
+
+ Sum of hypermet functions, defined by
+ (area, position, fwhm, st_area_r, st_slope_r, lt_area_r, lt_slope_r, step_height_r).
+
+ - *area* is the area underneath the gaussian peak
+ - *position* is the center of the various peaks and the position of
+ the step down
+ - *fwhm* is the full-width at half maximum of the terms
+ - *st_area_r* is factor between the gaussian area and the area of the
+ short tail term
+ - *st_slope_r* is a parameter related to the slope of the short tail
+ in the low ``x`` values (the lower, the steeper)
+ - *lt_area_r* is factor between the gaussian area and the area of the
+ long tail term
+ - *lt_slope_r* is a parameter related to the slope of the long tail
+ in the low ``x`` values (the lower, the steeper)
+ - *step_height_r* is the factor between the height of the step down
+ and the gaussian height
+
+ Parameters:
+ -----------
+
+ - x: Independant variable where the functions are calculated.
+ - len_x: Number of elements in the x array.
+ - phypermet: Array of hypermet function parameters:
+ *(area1, position1, fwhm1, st_area_r1, st_slope_r1, lt_area_r1,
+ lt_slope_r1, step_height_r1, ...)*
+ - len_phypermet: Number of elements in the phypermet array. Must be
+ a multiple of 8.
+ - y: Output array. Must have memory allocated for the same number
+ of elements as x (len_x).
+ - tail_flags: sum of binary flags to activate the various terms of the
+ function:
+
+ - 1 (b0001): Gaussian term
+ - 2 (b0010): st term
+ - 4 (b0100): lt term
+ - 8 (b1000): step term
+
+ E.g., to activate all termsof the hypermet, use ``tail_flags = 1 + 2 + 4 + 8 = 15``
+
+ Adapted from PyMca module SpecFitFuns
+*/
+int sum_fastahypermet(double* x, int len_x, double* phypermet, int len_phypermet, double* y, int tail_flags)
+{
+ int i, j;
+ int g_term_flag, st_term_flag, lt_term_flag, step_term_flag;
+ double c1, c2, sigma, height, sigma_sqrt2, sqrt2PI, inv_2_sqrt_2_log2, x_minus_position, epsilon;
+ double area, position, fwhm, st_area_r, st_slope_r, lt_area_r, lt_slope_r, step_height_r;
+
+ if (test_params(len_phypermet, 8, "sum_hypermet",
+ "height, centroid, fwhm, st_area_r, st_slope_r, lt_area_r, lt_slope_r, step_height_r")) {
+ return(1);
+ }
+
+ g_term_flag = tail_flags & 1;
+ st_term_flag = (tail_flags>>1) & 1;
+ lt_term_flag = (tail_flags>>2) & 1;
+ step_term_flag = (tail_flags>>3) & 1;
+
+ /* Initialize output array */
+ for (j=0; j<len_x; j++) {
+ y[j] = 0.;
+ }
+
+ /* define epsilon to compare floating point values with 0. */
+ epsilon = 0.00000000001;
+
+ sqrt2PI= sqrt(2.0 * M_PI);
+ inv_2_sqrt_2_log2 = 1.0 / (2.0 * sqrt(2.0 * LOG2));
+
+ for (i=0; i<len_phypermet/8; i++) {
+ area = phypermet[8*i];
+ position = phypermet[8*i+1];
+ fwhm = phypermet[8*i+2];
+ st_area_r = phypermet[8*i+3];
+ st_slope_r = phypermet[8*i+4];
+ lt_area_r = phypermet[8*i+5];
+ lt_slope_r = phypermet[8*i+6];
+ step_height_r = phypermet[8*i+7];
+
+ sigma = fwhm * inv_2_sqrt_2_log2;
+ height = area / (sigma * sqrt2PI);
+
+ /* Prevent division by 0 */
+ if (sigma == 0) {
+ printf("fwhm must not be equal to 0");
+ return(1);
+ }
+ sigma_sqrt2 = sigma * 1.4142135623730950488;
+
+ for (j=0; j<len_x; j++) {
+ x_minus_position = x[j] - position;
+ c2 = (0.5 * x_minus_position * x_minus_position) / (sigma * sigma);
+ /* gaussian term */
+ if (g_term_flag && c2 < 100) {
+ y[j] += fastexp(-c2) * height;
+ }
+
+ /* st term */
+ if (st_term_flag && (fabs(st_slope_r) > epsilon) && (x_minus_position / st_slope_r) <= 612) {
+ c1 = st_area_r * 0.5 * \
+ erfc((x_minus_position/sigma_sqrt2) + 0.5 * sigma_sqrt2 / st_slope_r);
+ y[j] += ((area * c1) / st_slope_r) * \
+ fastexp(0.5 * (sigma / st_slope_r) * (sigma / st_slope_r) +\
+ (x_minus_position / st_slope_r));
+ }
+
+ /* lt term */
+ if (lt_term_flag && (fabs(lt_slope_r) > epsilon) && (x_minus_position / lt_slope_r) <= 612) {
+ c1 = lt_area_r * \
+ 0.5 * erfc((x_minus_position/sigma_sqrt2) + 0.5 * sigma_sqrt2 / lt_slope_r);
+ y[j] += ((area * c1) / lt_slope_r) * \
+ fastexp(0.5 * (sigma / lt_slope_r) * (sigma / lt_slope_r) +\
+ (x_minus_position / lt_slope_r));
+
+ }
+
+ /* step term flag */
+ if (step_term_flag) {
+ y[j] += step_height_r * (area / (sigma * sqrt2PI)) *\
+ 0.5 * erfc(x_minus_position / sigma_sqrt2);
+ }
+ }
+ }
+ return(0);
+}
+
+void pileup(double* x, long len_x, double* ret, int input2, double zero, double gain)
+{
+ //int input2=0;
+ //double zero=0.0;
+ //double gain=1.0;
+
+ int i, j, k;
+ double *px, *pret, *pall;
+
+ /* the pointer to the starting position of par data */
+ px = x;
+ pret = ret;
+
+ *pret = 0;
+ k = (int )(zero/gain);
+ for (i=input2; i<len_x; i++){
+ pall = x;
+ if ((i+k) >= 0)
+ {
+ pret = (double *) ret+(i+k);
+ for (j=0; j<len_x-i-k ;j++){
+ *pret += *px * (*pall);
+ pall++;
+ pret++;
+ }
+ }
+ px++;
+ }
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-26 11:33:19 +0000