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){
63 gmx_bool convolution(int dataSize,real *x, int kernelSize, real* kernel)
68 /* check validity of params */
69 if(!x || !kernel) return FALSE;
70 if(dataSize <=0 || kernelSize <= 0) return FALSE;
72 /* start convolution from out[kernelSize-1] to out[dataSize-1] (last) */
73 for(i = kernelSize-1; i < dataSize; ++i)
75 for(j = i, k = 0; k < kernelSize; --j, ++k)
76 out[i] += x[j] * kernel[k];
79 /* convolution from out[0] to out[kernelSize-2] */
80 for(i = 0; i < kernelSize - 1; ++i)
82 for(j = i, k = 0; j >= 0; --j, ++k)
83 out[i] += x[j] * kernel[k];
86 for(i=0;i<dataSize;i++) x[i]=out[i];
91 /* Assuming kernel is shorter than x */
93 gmx_bool periodic_convolution(int datasize, real *x, int kernelsize,
99 if (!x || !kernel) return FALSE;
100 if (kernelsize<=0|| datasize<=0|| kernelsize > datasize) return FALSE;
102 snew(filtered,datasize);
104 for(i=0;(i<datasize); i++){
105 for(j=0; (j<kernelsize); j++){
106 filtered[i] += kernel[j]*x[MOD((i-j),datasize)];
109 for(i=0; i<datasize; i++) x[i]=filtered[i];
116 /* 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
117 * NO checks are performed
119 void gausskernel(real *out, int n, real var){
124 for(i=-k; i<=k; i++){
126 tot+=out[j++]=EXP(-arg);