src/tools/powerspect.c

   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
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.0
  11  *
  12  * Copyright (c) 1991-2001
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  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.
  20  *
  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.
  27  *
  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.
  30  *
  31  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  32  *
  33  * And Hey:
  34  * Gyas ROwers Mature At Cryogenic Speed
  35  */
  36
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40 #include <smalloc.h>
  41 #include <gmx_fft.h>
  42 #include <gmx_fatal.h>
  43 #include <futil.h>
  44 #include "interf.h"
  45 #include "powerspect.h"
  46
  47 void addtoavgenergy(t_complex *list, real *result, int size, int tsteps){
  48   int i;
  49         for (i=0;i<size;i++){
  50                 result[i]+=cabs2(list[i])/tsteps;
  51         }
  52
  53 }
  54
  55
  56 void powerspectavg(real ***intftab, int tsteps, int xbins, int ybins, char **outfiles){
  57   /*Fourier plans and output;*/
  58   gmx_fft_t  fftp;
  59  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 */
  60   t_complex *ftspect2;
  61   real *pspectavg1; /*power -spectrum 1st interface*/
  62   real *pspectavg2; /* -------------- 2nd interface*/
  63   real *temp;
  64   FILE *datfile1,*datfile2; /*data-files with interface data*/
  65   int n; /*time index*/
  66   int fy=ybins/2+1; /* number of (symmetric) fourier y elements; */
  67   int rfl=xbins*fy; /*length of real - DFT == Symmetric 2D matrix*/
  68   int status;
  69
  70 /*Prepare data structures for FFT, with time averaging of power spectrum*/
  71 if ( (status=gmx_fft_init_2d_real(&fftp,xbins,ybins,GMX_FFT_FLAG_NONE) )!=0)
  72   {
  73         free(fftp);
  74         gmx_fatal(status,__FILE__,__LINE__,"Error allocating FFT");
  75     }
  76
  77 /*Initialize arrays*/
  78 snew(ftspect1,rfl);
  79 snew(ftspect2,rfl);
  80 snew(temp,rfl);
  81 snew(pspectavg1,rfl);
  82 snew(pspectavg2,rfl);
  83
  84 /*Fouriertransform directly (no normalization or anything)*/
  85 /*NB! Check carefully indexes here*/
  86
  87         for (n=0;n<tsteps;n++){
  88                 gmx_fft_2d_real(fftp,GMX_FFT_REAL_TO_COMPLEX,intftab[0][n],ftspect1);
  89                 gmx_fft_2d_real(fftp,GMX_FFT_REAL_TO_COMPLEX,intftab[1][n],ftspect2);
  90
  91                 /*Add to average for interface 1 here*/
  92                 addtoavgenergy(ftspect1,pspectavg1,rfl,tsteps);
  93                 /*Add to average for interface 2 here*/
  94                 addtoavgenergy(ftspect2,pspectavg2,rfl,tsteps);
  95         }
  96         /*Print out average energy-spectrum to outfiles[0] and outfiles[1];*/
  97
  98 datfile1 = ffopen(outfiles[0],"w");
  99 datfile2 = ffopen(outfiles[1],"w");
 100
 101 /*Filling dat files with spectral data*/
 102 fprintf(datfile1,"%s\n","kx\t ky\t\tPower(kx,ky)");
 103 fprintf(datfile2,"%s\n","kx\t ky\t\tPower(kx,ky)");
 104         for(n=0;n<rfl;n++){
 105                 fprintf(datfile1,"%d\t%d\t %8.6f\n",(n / fy),(n % fy),pspectavg1[n]);
 106                 fprintf(datfile2,"%d\t%d\t %8.6f\n",(n /fy),(n % fy),pspectavg2[n]);
 107         }
 108 ffclose(datfile1);
 109 ffclose(datfile2);
 110
 111 free(ftspect1);
 112 free(ftspect2);
 113
 114 }
 115
 116 void powerspectavg_intf(t_interf ***if1, t_interf ***if2, int t, int xb, int yb, char **fnms)
 117 {
 118         real ***surf;
 119
 120         int xy=xb*yb;
 121         int i,n;
 122
 123         snew(surf,2);
 124         snew(surf[0],t);
 125         snew(surf[1],t);
 126         for (n=0;n<t;n++){
 127                 snew(surf[0][n],xy);
 128                 snew(surf[1][n],xy);
 129                 for(i=0;i<xy;i++){
 130                         surf[0][n][i]=if1[n][i]->Z;
 131                         surf[1][n][i]=if2[n][i]->Z;
 132                 }
 133         }
 134         powerspectavg(surf, t,xb,yb,fnms);
 135 }
 136