[Scipy-svn] r2848 - in trunk/Lib/stsci/convolve: . lib src

scipy-svn@scip... scipy-svn@scip...
Wed Mar 14 11:53:04 CDT 2007


Author: chanley
Date: 2007-03-14 11:52:59 -0500 (Wed, 14 Mar 2007)
New Revision: 2848

Added:
   trunk/Lib/stsci/convolve/lib/lineshape.py
   trunk/Lib/stsci/convolve/src/_lineshapemodule.c
Modified:
   trunk/Lib/stsci/convolve/setup.py
Log:
Adding numpy version of numarray's lineshape module from the convolve package.

Added: trunk/Lib/stsci/convolve/lib/lineshape.py
===================================================================
--- trunk/Lib/stsci/convolve/lib/lineshape.py	2007-03-14 12:46:49 UTC (rev 2847)
+++ trunk/Lib/stsci/convolve/lib/lineshape.py	2007-03-14 16:52:59 UTC (rev 2848)
@@ -0,0 +1,102 @@
+# lineshape functors
+# -*- coding: iso-8859-1 -*-
+#
+# Copyright (C) 2002 Jochen Küpper <jochen@jochen-kuepper.de>
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice,
+#    this list of conditions and the following disclaimer.
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+# 3. The name of the author may not be used to endorse or promote products
+#    derived from this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
+# WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+# EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+__doc__     = """Lineshape functors.
+
+The objects defined in this module can be used to calculate numarrays containing
+common lineshape-profiles.
+
+For the *Profile classes the profile is only evaluated once for each object and
+then reused for each call. If you only use a profile once and you are concerned
+about memory consumption you could call the underlying functions directly.
+"""
+
+__author__  = "Jochen Küpper <jochen@jochen-kuepper.de>"
+__date__    = "$Date: 2007/03/14 16:35:57 $"[7:-11]
+__version__ = "$Revision: 1.1 $"[11:-2]
+
+
+import numpy as num
+from convolve._lineshape import *
+
+
+class Profile(object):
+    """An base object to provide a convolution kernel."""
+
+    def __init__(self, x, w, x0=0.0):
+        # call init for all superclasses
+        super(Profile, self).__init__(x, w, x0)
+        self._recalculate(x, w, x0)
+
+    def __call__(self):
+        return self._kernel
+
+    def _recalculate(self, x, w, x0):
+        self._kernel = self._profile(x, w, x0)
+
+
+class GaussProfile(Profile):
+    """An object for Gauss-folding."""
+
+    def __init__(self, x, w, x0=0.0):
+        self._profile = gauss
+        # call init for all superclasses
+        super(GaussProfile, self).__init__(x, w, x0)
+
+
+class LorentzProfile(Profile):
+    """An object for Lorentz-folding."""
+
+    def __init__(self, x, w, x0=0.0):
+        self._profile = lorentz
+        # call init for all superclasses
+        super(LorentzProfile, self).__init__(x, w, x0)
+
+
+
+class VoigtProfile(Profile):
+    """An object for Voigt-folding.
+
+    The constructor takes the following parameter:
+      |x|  Scalar or numarray with values to calculate profile at.
+      |w|  Tuple of Gaussian and Lorentzian linewidth contribution
+      |x0| Center frequency
+    """
+
+    def __init__(self, x, w, x0=0.0):
+        self._profile = voigt
+        # call init for all superclasses
+        super(VoigtProfile, self).__init__(x, w, x0)
+
+
+
+## Local Variables:
+## mode: python
+## mode: auto-fill
+## fill-column: 80
+## End:

Modified: trunk/Lib/stsci/convolve/setup.py
===================================================================
--- trunk/Lib/stsci/convolve/setup.py	2007-03-14 12:46:49 UTC (rev 2847)
+++ trunk/Lib/stsci/convolve/setup.py	2007-03-14 16:52:59 UTC (rev 2848)
@@ -10,6 +10,10 @@
                          sources=["src/_correlatemodule.c"],
                          define_macros = [('NUMPY', '1')],
                          include_dirs = [numpy.get_numarray_include()])
+    config.add_extension('_lineshape',
+                         sources=["src/_lineshapemodule.c"],
+                         define_macros = [('NUMPY', '1')],
+                         include_dirs = [numpy.get_numarray_include()])
     return config
 
 if __name__ == "__main__":

Added: trunk/Lib/stsci/convolve/src/_lineshapemodule.c
===================================================================
--- trunk/Lib/stsci/convolve/src/_lineshapemodule.c	2007-03-14 12:46:49 UTC (rev 2847)
+++ trunk/Lib/stsci/convolve/src/_lineshapemodule.c	2007-03-14 16:52:59 UTC (rev 2848)
@@ -0,0 +1,381 @@
+/* C implementations of various lineshape functions
+ *
+ * Copyright (C) 2002,2003 Jochen Küpper <jochen@jochen-kuepper.de>
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice,
+ *    this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ * 3. The name of the author may not be used to endorse or promote products
+ *    derived from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+ * EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+ * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+ * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "Python.h"
+
+#include <math.h>
+
+#include "numpy/libnumarray.h"
+
+
+#define sqr(x) ((x)*(x))
+
+
+/* These are apparently not defined in MSVC */
+#if !defined(M_PI)
+#define M_PI 3.14159265358979323846
+#endif
+#if !defined(M_LN2)
+#define M_LN2 0.69314718055994530942
+#endif
+
+
+/*** C implementation ***/
+
+static void gauss(size_t n, double *x, double *y, double w, double xc) 
+    /* Evaluate normalized Gauss profile around xc with FWHM w at all x_i,
+       return in y. */
+{
+    int i;
+    for(i=0; i<n; i++)
+        y[i] = 2. * sqrt(M_LN2/M_PI) / w * exp(-4.*M_LN2 * sqr((x[i]-xc)/w));
+}
+
+
+
+static void lorentz(size_t n, double *x, double *y, double w, double xc) 
+    /* Evaluate normalized Lorentz profile around xc with FWHM w at all x_i,
+       return in y. */
+{
+    int i;
+    for(i=0; i<n; i++)
+        y[i] = 2.*w/M_PI / (sqr(w) + 4.*(sqr(x[i]-xc)));
+}
+
+
+
+static double humlicek_v12(double x, double y)
+    /** Approximation of Voigt profile by Humlicek's 12-point formula.
+     *
+     * J. Humlicek, J. Quant. Spectrosc. Radiat. Transfer, 21(1978), 309.
+     *
+     * Voigt-Profil:
+     * V(x, y) = 2/pi^(1.5) * y^2/FWHM_L * \int[-inf,+inf](e^(-y^2)/(x+y)^2+...)
+     */
+{
+    static const double T_v12[6] = {
+	0.314240376254359,    0.947788391240164,    1.597682635152605,
+	2.27950708050106,     3.020637025120890,    3.889724897869782
+    };
+    static const double alpha_v12[6] = {
+       -1.393236997981977,   -0.231152406188676,    0.155351465642094,
+       -6.21836623696556e-3, -9.190829861057113e-5, 6.275259577497896e-7
+    };
+    static const double beta_v12[6] = {
+	1.011728045548831,   -0.751971469674635,    1.255772699323164e-2,
+	1.0022008145159e-2,  -2.42068134815573e-4,  5.008480613664573e-7
+    };
+    static const double  y0_v12 = 1.50;
+    double yp, xp, xm, sum, yp2, xp2, xm2;
+    int k;
+
+    sum = 0.;
+    yp = y + y0_v12;
+    yp2 = yp * yp;
+    if((y > 0.85) || (fabs(x) < (18.1 * y + 1.65))) {
+	/* Bereich I */
+	for(k=0; k<6; k++) {
+	    xp = x + T_v12[k];
+	    xm = x - T_v12[k];
+	    sum += ((alpha_v12[k] * xm + beta_v12[k] * yp) / (xm * xm + yp2)
+                    + (beta_v12[k] * yp - alpha_v12[k] * xp) / (xp * xp + yp2));
+	}
+    } else {
+	/* Bereich II */
+	for(k=0; k<6; k++) {
+	    xp = x + T_v12[k];
+	    xp2 = xp * xp;
+	    xm = x - T_v12[k];
+	    xm2 = xm * xm;
+	    sum += (((beta_v12[k] * (xm2 - y0_v12 * yp) - alpha_v12[k] * xm * (yp + y0_v12))
+                     / ((xm2 + yp2) * (xm2 + y0_v12 * y0_v12)))
+                    + ((beta_v12[k] * (xp2 - y0_v12 * yp) + alpha_v12[k] * xp * (yp + y0_v12))
+                       / ((xp2 + yp2) * (xp2 + y0_v12 * y0_v12))));
+	}
+	if(fabs(x) < 100.)
+	    sum = y * sum + exp(-pow(x, 2));
+	else
+	    sum *= y;
+    }
+    return sum;
+}
+
+
+static void voigt(size_t n, double *x, double *y, double w[2], double xc) 
+    /* Evaluate normalized Voigt profile at x around xc with Gaussian
+     * linewidth contribution w[0] and Lorentzian linewidth
+     * contribution w[1].
+     */
+{
+    /* Transform into reduced coordinates and call Humlicek's 12 point
+     * formula:
+     *     x = 2 \sqrt{\ln2} \frac{\nu-\nu_0}{\Delta\nu_G}
+     *     y = \sqrt{\ln2} \frac{\Delta\nu_L}{\Delta\nu_G}
+     */
+    int i;
+    double yh = sqrt(M_LN2) * w[1] / w[0];
+    for(i=0; i<n; i++) {
+        double xh = 2. * sqrt(M_LN2) * (x[i]-xc) / w[0];
+        y[i] = 2.*sqrt(M_LN2/M_PI)/w[0] * humlicek_v12(xh, yh);
+    }
+}
+
+
+
+/*** Python interface ***/
+
+static PyObject *_Error;
+
+
+static PyObject * 
+_lineshape_gauss(PyObject *self, PyObject *args, PyObject *keywds) 
+{
+    int f;
+    double w, xc = 0.0;
+    static char *kwlist[] = {"x", "w", "xc", "y", NULL};
+    PyObject *ox, *oy=Py_None;
+    PyArrayObject *x, *y;
+
+    if(! PyArg_ParseTupleAndKeywords(args, keywds, "Od|dO", kwlist,
+                                     &ox, &w, &xc, &oy))
+        return PyErr_Format(PyExc_RuntimeError,  "gauss: invalid parameters");
+
+    if((f = PyFloat_Check(ox)) || PyInt_Check(ox)) {
+        /* scalar arguments -- always *return* Float result */
+        double xa[1], ya[1];
+        if(f)
+            xa[0] = PyFloat_AS_DOUBLE(ox);
+        else
+            xa[0] = (double)PyInt_AS_LONG(ox);
+        Py_BEGIN_ALLOW_THREADS;
+        gauss(1, xa, ya, w, xc);
+        Py_END_ALLOW_THREADS;
+        Py_DECREF(ox);
+        return PyFloat_FromDouble(ya[0]);
+    } else {
+        /* array conversion */        
+        if(! ((x = NA_InputArray(ox, tFloat64, C_ARRAY))
+              && (y = NA_OptionalOutputArray(oy, tFloat64, C_ARRAY, x))))
+            return 0;
+        if(x->nd != 1)
+            return PyErr_Format(_Error, "gauss: x must be scalar or 1d array.");
+        if (!NA_ShapeEqual(x, y))
+            return PyErr_Format(_Error, "gauss: x and y numarray must have same length.");
+
+        /* calculate profile */
+	{
+        double *xa = NA_OFFSETDATA(x);
+        double *ya = NA_OFFSETDATA(y);
+        Py_BEGIN_ALLOW_THREADS;
+        gauss(x->dimensions[0], xa, ya, w, xc);
+        Py_END_ALLOW_THREADS;
+	}
+    
+        /* cleanup and return */
+        Py_XDECREF(x);
+        return NA_ReturnOutput(oy, y);
+    }
+}
+
+
+
+static PyObject * 
+_lineshape_lorentz(PyObject *self, PyObject *args, PyObject *keywds) 
+{
+    int f;
+    double w, xc = 0.0;
+    static char *kwlist[] = {"x", "w", "xc", "y", NULL};
+    PyObject *ox, *oy=Py_None;
+    PyArrayObject *x, *y;
+
+    if(! PyArg_ParseTupleAndKeywords(args, keywds, "Od|dO", kwlist,
+                                     &ox, &w, &xc, &oy))
+        return PyErr_Format(PyExc_RuntimeError,  "lorentz: invalid parameters");
+
+    if((f = PyFloat_Check(ox)) || PyInt_Check(ox)) {
+        /* scalar arguments -- always *return* Float result */
+        double xa[1], ya[1];
+        if(f)
+            xa[0] = PyFloat_AS_DOUBLE(ox);
+        else
+            xa[0] = (double)PyInt_AS_LONG(ox);
+        Py_BEGIN_ALLOW_THREADS;
+        lorentz(1, xa, ya, w, xc);
+        Py_END_ALLOW_THREADS;
+        Py_DECREF(ox);
+        return PyFloat_FromDouble(ya[0]);
+    } else {
+        /* array conversion */        
+        if(! ((x = NA_InputArray(ox, tFloat64, C_ARRAY))
+              && (y = NA_OptionalOutputArray(oy, tFloat64, C_ARRAY, x))))
+            return 0;
+        if(x->nd != 1)
+            return PyErr_Format(_Error, "lorentz: x must be scalar or 1d array.");
+        if (!NA_ShapeEqual(x, y))
+            return PyErr_Format(_Error, "lorentz: x and y numarray must have same length.");
+
+        /* calculate profile */
+	{
+        double *xa = NA_OFFSETDATA(x);
+        double *ya = NA_OFFSETDATA(y);
+
+        Py_BEGIN_ALLOW_THREADS;
+        lorentz(x->dimensions[0], xa, ya, w, xc);
+        Py_END_ALLOW_THREADS;
+	}
+
+        /* cleanup and return */
+        Py_XDECREF(x);
+        return NA_ReturnOutput(oy, y);
+    }
+}
+
+
+
+static PyObject * 
+_lineshape_voigt(PyObject *self, PyObject *args, PyObject *keywds) 
+{
+    int f;
+    double w[2], xc = 0.0;
+    static char *kwlist[] = {"x", "w", "xc", "y", NULL};
+    PyObject *wt, *ox, *oy=Py_None;
+    PyArrayObject *x, *y;
+
+    if(! PyArg_ParseTupleAndKeywords(args, keywds, "OO|dO", kwlist,
+                                     &ox, &wt, &xc, &oy))
+        return PyErr_Format(PyExc_RuntimeError,  "voigt: invalid parameters");
+
+    /* parse linewidths tuple */
+    if(! PyArg_ParseTuple(wt, "dd", &(w[0]), &(w[1])))
+        return(0);
+
+    if((f = PyFloat_Check(ox)) || PyInt_Check(ox)) {
+        /* scalar arguments -- always *return* Float result */
+        double xa[1], ya[1];
+        if(f)
+            xa[0] = PyFloat_AS_DOUBLE(ox);
+        else
+            xa[0] = (double)PyInt_AS_LONG(ox);
+        Py_BEGIN_ALLOW_THREADS;
+        voigt(1, xa, ya, w, xc);
+        Py_END_ALLOW_THREADS;
+        Py_DECREF(ox);
+        return PyFloat_FromDouble(ya[0]);
+    } else {
+        /* array conversion */        
+        if(! ((x = NA_InputArray(ox, tFloat64, C_ARRAY))
+              && (y = NA_OptionalOutputArray(oy, tFloat64, C_ARRAY, x))))
+            return 0;
+        if(x->nd != 1)
+            return PyErr_Format(_Error, "voigt: x must be scalar or 1d array.");
+        if (!NA_ShapeEqual(x, y))
+            return PyErr_Format(_Error, "voigt: x and y numarray must have same length.");
+
+        /* calculate profile */
+	{
+        double *xa = NA_OFFSETDATA(x);
+        double *ya = NA_OFFSETDATA(y);
+        Py_BEGIN_ALLOW_THREADS;
+        voigt(x->dimensions[0], xa, ya, w, xc);
+        Py_END_ALLOW_THREADS;
+	}
+
+        /* cleanup and return */
+        Py_XDECREF(x);
+        return NA_ReturnOutput(oy, y);
+    }
+}
+
+
+
+
+/*** table of methods ***/
+
+static PyMethodDef _lineshape_Methods[] = {
+    {"gauss", (PyCFunction)_lineshape_gauss, METH_VARARGS|METH_KEYWORDS,
+     "gauss(x, w, xc=0.0, y=None)\n\n"
+     "Gaussian lineshape function\n\n" \
+     "Calculate normalized Gaussian with full-width at half maximum |w| at |x|,\n" \
+     "optionally specifying the line-center |xc|.\n" \
+     "If, and only if |x| is an array an optional output array |y| can be\n" \
+     "specified.  In this case |x| and |y| must be one-dimensional numarray\n" \
+     "with identical shapes.\n\n" \
+     "If |x| is an scalar the routine always gives the result as scalar\n" \
+     "return value."
+    },
+    {"lorentz", (PyCFunction)_lineshape_lorentz, METH_VARARGS|METH_KEYWORDS,
+     "lorentz(x, w, xc=0.0, y=None)\n\n"
+     "Lorentzian lineshape function\n\n" \
+     "Calculate normalized Lorentzian with full-width at half maximum |w| at |x|,\n" \
+     "optionally specifying the line-center |xc|.\n" \
+     "If, and only if |x| is an array an optional output array |y| can be\n" \
+     "specified.  In this case |x| and |y| must be one-dimensional numarray\n" \
+     "with identical shapes.\n\n" \
+     "If |x| is an scalar the routine always gives the result as scalar\n" \
+     "return value."
+    },
+    {"voigt", (PyCFunction)_lineshape_voigt, METH_VARARGS|METH_KEYWORDS,
+     "voigt(x, w, xc=0.0, y=None)\n\n"
+     "Voigt-lineshape function\n\n" \
+     "Calculate normalized Voigt-profile with Gaussian full-width at half maximum |w[0]| and\n" \
+     "Lorentzian full-width at half maximum |w[1]| at |x|, optionally specifying the line-center\n" \
+     "|xc|.\n" \
+     "If, and only if |x| is an array an optional output array |y| can be\n" \
+     "specified.  In this case |x| and |y| must be one-dimensional numarray\n" \
+     "with identical shapes.\n\n" \
+     "If |x| is an scalar the routine always gives the result as scalar\n" \
+     "return value.\n\n" \
+     "This function uses Humlicek's 12-point formula to approximate the Voigt\n" \
+     "profile (J. Humlicek, J. Quant. Spectrosc. Radiat. Transfer, 21, 309 (1978))."
+    },
+    {NULL, NULL, 0, ""}
+};
+
+
+
+
+/*** module initialization ***/
+
+DL_EXPORT(void) init_lineshape(void)
+{
+    PyObject *m, *d;
+    m = Py_InitModule("_lineshape", _lineshape_Methods);
+    d = PyModule_GetDict(m);
+    _Error = PyErr_NewException("_lineshape.error", NULL, NULL);
+    PyDict_SetItemString(d, "error", _Error);
+    import_libnumarray();
+}
+
+
+
+/*
+ * Local Variables:
+ * mode: c
+ * c-file-style: "Stroustrup"
+ * fill-column: 80
+ * End:
+ */



More information about the Scipy-svn mailing list