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
42 #include "gromacs/fft/fft.h"
43 #include "gmx_fatal.h"
44 #include "gromacs/fileio/futil.h"
46 #include "powerspect.h"
48 void addtoavgenergy(t_complex *list, real *result, int size, int tsteps)
51 for (i = 0; i < size; i++)
53 result[i] += cabs2(list[i])/tsteps;
59 void powerspectavg(real ***intftab, int tsteps, int xbins, int ybins, char **outfiles)
61 /*Fourier plans and output;*/
63 t_complex *ftspect1; /* Spatial FFT of interface for each time frame and interface ftint[time,xycoord][0], ftintf[time,xycoord][1] for interface 1 and 2 respectively */
65 real *pspectavg1; /*power -spectrum 1st interface*/
66 real *pspectavg2; /* -------------- 2nd interface*/
68 FILE *datfile1, *datfile2; /*data-files with interface data*/
70 int fy = ybins/2+1; /* number of (symmetric) fourier y elements; */
71 int rfl = xbins*fy; /*length of real - DFT == Symmetric 2D matrix*/
74 /*Prepare data structures for FFT, with time averaging of power spectrum*/
75 if ( (status = gmx_fft_init_2d_real(&fftp, xbins, ybins, GMX_FFT_FLAG_NONE) ) != 0)
78 gmx_fatal(status, __FILE__, __LINE__, "Error allocating FFT");
85 snew(pspectavg1, rfl);
86 snew(pspectavg2, rfl);
88 /*Fouriertransform directly (no normalization or anything)*/
89 /*NB! Check carefully indexes here*/
91 for (n = 0; n < tsteps; n++)
93 gmx_fft_2d_real(fftp, GMX_FFT_REAL_TO_COMPLEX, intftab[0][n], ftspect1);
94 gmx_fft_2d_real(fftp, GMX_FFT_REAL_TO_COMPLEX, intftab[1][n], ftspect2);
96 /*Add to average for interface 1 here*/
97 addtoavgenergy(ftspect1, pspectavg1, rfl, tsteps);
98 /*Add to average for interface 2 here*/
99 addtoavgenergy(ftspect2, pspectavg2, rfl, tsteps);
101 /*Print out average energy-spectrum to outfiles[0] and outfiles[1];*/
103 datfile1 = ffopen(outfiles[0], "w");
104 datfile2 = ffopen(outfiles[1], "w");
106 /*Filling dat files with spectral data*/
107 fprintf(datfile1, "%s\n", "kx\t ky\t\tPower(kx,ky)");
108 fprintf(datfile2, "%s\n", "kx\t ky\t\tPower(kx,ky)");
109 for (n = 0; n < rfl; n++)
111 fprintf(datfile1, "%d\t%d\t %8.6f\n", (n / fy), (n % fy), pspectavg1[n]);
112 fprintf(datfile2, "%d\t%d\t %8.6f\n", (n /fy), (n % fy), pspectavg2[n]);
122 void powerspectavg_intf(t_interf ***if1, t_interf ***if2, int t, int xb, int yb, char **fnms)
132 for (n = 0; n < t; n++)
134 snew(surf[0][n], xy);
135 snew(surf[1][n], xy);
136 for (i = 0; i < xy; i++)
138 surf[0][n][i] = if1[n][i]->Z;
139 surf[1][n][i] = if2[n][i]->Z;
142 powerspectavg(surf, t, xb, yb, fnms);