src/contrib/mk6_n.c

   1 /*
   2  * $Id$
   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.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Great Red Owns Many ACres of Sand
  31  */
  32
  33 #include <stdio.h>
  34 #include <stdlib.h>
  35 #include <math.h>
  36
  37 void doit(char *fn,double myexp[],int n,double tabscale)
  38 {
  39   int    i,k;
  40   double myfac[3] = { 1, -1, 1 };
  41   double x,v,v2;
  42   FILE   *fp;
  43
  44   fp = fopen(fn,"w");
  45   for(i=0; (i<=n); i++) {
  46     x   =  i/tabscale;
  47
  48     fprintf(fp,"%10g",x);
  49
  50     for(k=0; (k<3); k++) {
  51       if (x < 0.04) {
  52         /* Avoid very high numbers */
  53         v = v2 = 0;
  54       }
  55       else {
  56         v  =  myfac[k]*pow(x,-myexp[k]);
  57         v2 = (myexp[k]+1)*(myexp[k])*v/(x*x);
  58       }
  59       fprintf(fp,"  %10g  %10g",v,v2);
  60     }
  61     fprintf(fp,"\n");
  62   }
  63   fclose(fp);
  64 }
  65
  66 int main(int argc,char *argv[])
  67 {
  68   double my8[3]  = { 1, 6, 8 };
  69   double my9[3]  = { 1, 6, 9 };
  70   double my10[3] = { 1, 6, 10 };
  71   double my11[3] = { 1, 6, 11 };
  72   double my12[3] = { 1, 6, 12 };
  73 #ifdef DOUBLE
  74   double tabscale = 2000;
  75 #else
  76   double tabscale = 500;
  77 #endif
  78   int n = (int) (3.0*tabscale);
  79
  80   doit("table6-8.xvg",my8,n,tabscale);
  81   doit("table6-9.xvg",my9,n,tabscale);
  82   doit("table6-10.xvg",my10,n,tabscale);
  83   doit("table6-11.xvg",my11,n,tabscale);
  84   doit("table6-12.xvg",my12,n,tabscale);
  85
  86   return 0;
  87 }
  88