1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 * $Id: densorder.c,v 0.9
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-2001
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
34 * Gyas ROwers Mature At Cryogenic Speed
41 * Routines for Filters and convolutions
46 #include "dens_filter.h"
51 #define EXP(x) (exp(x))
53 #define EXP(x) (expf(x))
56 /* Modulo function assuming y>0 but with arbitrary INTEGER x */
57 static int MOD(int x, int y)
67 gmx_bool convolution(int dataSize, real *x, int kernelSize, real* kernel)
72 /* check validity of params */
77 if (dataSize <= 0 || kernelSize <= 0)
82 /* start convolution from out[kernelSize-1] to out[dataSize-1] (last) */
83 for (i = kernelSize-1; i < dataSize; ++i)
85 for (j = i, k = 0; k < kernelSize; --j, ++k)
87 out[i] += x[j] * kernel[k];
91 /* convolution from out[0] to out[kernelSize-2] */
92 for (i = 0; i < kernelSize - 1; ++i)
94 for (j = i, k = 0; j >= 0; --j, ++k)
96 out[i] += x[j] * kernel[k];
100 for (i = 0; i < dataSize; i++)
108 /* Assuming kernel is shorter than x */
110 gmx_bool periodic_convolution(int datasize, real *x, int kernelsize,
120 if (kernelsize <= 0 || datasize <= 0 || kernelsize > datasize)
125 snew(filtered, datasize);
127 for (i = 0; (i < datasize); i++)
129 for (j = 0; (j < kernelsize); j++)
131 filtered[i] += kernel[j]*x[MOD((i-j), datasize)];
134 for (i = 0; i < datasize; i++)
144 /* returns discrete gaussian kernel of size n in *out, where n=2k+1=3,5,7,9,11 and k=1,2,3 is the order
145 * NO checks are performed
147 void gausskernel(real *out, int n, real var)
153 for (i = -k; i <= k; i++)
156 tot += out[j++] = EXP(-arg);
158 for (i = 0; i < j; i++)