src/tools/lsq.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_lsq_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "gstat.h"
  33 #include "vec.h"
  34
  35 void init_lsq(t_lsq *lsq)
  36 {
  37   lsq->yy=lsq->yx=lsq->xx=lsq->sx=lsq->sy=0.0;
  38   lsq->np=0;
  39 }
  40
  41 int npoints_lsq(t_lsq *lsq)
  42 {
  43   return lsq->np;
  44 }
  45
  46 void done_lsq(t_lsq *lsq)
  47 {
  48   init_lsq(lsq);
  49 }
  50
  51 void add_lsq(t_lsq *lsq,real x,real y)
  52 {
  53   lsq->yy+=y*y;
  54   lsq->yx+=y*x;
  55   lsq->xx+=x*x;
  56   lsq->sx+=x;
  57   lsq->sy+=y;
  58   lsq->np++;
  59 }
  60
  61 void get_lsq_ab(t_lsq *lsq,real *a,real *b)
  62 {
  63   double yx,xx,sx,sy;
  64
  65   yx=lsq->yx/lsq->np;
  66   xx=lsq->xx/lsq->np;
  67   sx=lsq->sx/lsq->np;
  68   sy=lsq->sy/lsq->np;
  69
  70   (*a)=(yx-sx*sy)/(xx-sx*sx);
  71   (*b)=(sy)-(*a)*(sx);
  72 }
  73
  74 real aver_lsq(t_lsq *lsq)
  75 {
  76   if (lsq->np == 0)
  77     fatal_error(0,"No points in distribution\n");
  78
  79   return (lsq->sy/lsq->np);
  80 }
  81
  82 real sigma_lsq(t_lsq *lsq)
  83 {
  84   if (lsq->np == 0)
  85     fatal_error(0,"No points in distribution\n");
  86
  87   return sqrt(lsq->yy/lsq->np - sqr(lsq->sy/lsq->np));
  88 }
  89
  90 real error_lsq(t_lsq *lsq)
  91 {
  92   return sigma_lsq(lsq)/sqrt(lsq->np);
  93 }
  94