src/gromacs/gmxana/expfit.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "xvgr.h"
  47 #include "gromacs/fileio/futil.h"
  48 #include "gstat.h"
  49 #include "vec.h"
  50 #include "index.h"
  51
  52 const int   nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3};
  53
  54 const char *s_ffn[effnNR+2] = {
  55     NULL, "none", "exp", "aexp", "exp_exp", "vac",
  56     "exp5", "exp7", "exp9", "erffit", NULL, NULL
  57 };
  58 /* We don't allow errest as a choice on the command line */
  59
  60 const char *longs_ffn[effnNR] = {
  61     "no fit",
  62     "y = exp(-x/a1)",
  63     "y = a2 exp(-x/a1)",
  64     "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
  65     "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
  66     "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5",
  67     "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
  68     "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
  69     "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)",
  70     "y = a2*ee(a1,x) + (1-a2)*ee(a2,x)"
  71 };
  72
  73 extern gmx_bool mrqmin_new(real x[], real y[], real sig[], int ndata, real a[],
  74                            int ia[], int ma, real **covar, real **alpha, real *chisq,
  75                            void (*funcs)(real, real [], real *, real []),
  76                            real *alamda);
  77
  78 static real myexp(real x, real A, real tau)
  79 {
  80     if ((A == 0) || (tau == 0))
  81     {
  82         return 0;
  83     }
  84     return A*exp(-x/tau);
  85 }
  86
  87 void erffit (real x, real a[], real *y, real dyda[])
  88 {
  89 /* Fuction
  90  *      y=(a1+a2)/2-(a1-a2)/2*erf((x-a3)/a4^2)
  91  */
  92
  93     real erfarg;
  94     real erfval;
  95     real erfarg2;
  96     real derf;
  97
  98     erfarg  = (x-a[3])/(a[4]*a[4]);
  99     erfarg2 = erfarg*erfarg;
 100     erfval  = gmx_erf(erfarg)/2;
 101     derf    = M_2_SQRTPI*(a[1]-a[2])/2*exp(-erfarg2)/(a[4]*a[4]);
 102     *y      = (a[1]+a[2])/2-(a[1]-a[2])*erfval;
 103     dyda[1] = 1/2-erfval;
 104     dyda[2] = 1/2+erfval;
 105     dyda[3] = derf;
 106     dyda[4] = 2*derf*erfarg;
 107 }
 108
 109
 110
 111 static void exp_one_parm(real x, real a[], real *y, real dyda[])
 112 {
 113     /* Fit to function
 114      *
 115      * y = exp(-x/a1)
 116      *
 117      */
 118
 119     real e1;
 120
 121     e1      = exp(-x/a[1]);
 122     *y      = e1;
 123     dyda[1] = x*e1/sqr(a[1]);
 124 }
 125
 126 static void exp_two_parm(real x, real a[], real *y, real dyda[])
 127 {
 128     /* Fit to function
 129      *
 130      * y = a2 exp(-x/a1)
 131      *
 132      */
 133
 134     real e1;
 135
 136     e1      = exp(-x/a[1]);
 137     *y      = a[2]*e1;
 138     dyda[1] = x*a[2]*e1/sqr(a[1]);
 139     dyda[2] = e1;
 140 }
 141
 142 static void exp_3_parm(real x, real a[], real *y, real dyda[])
 143 {
 144     /* Fit to function
 145      *
 146      * y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)
 147      *
 148      */
 149
 150     real e1, e2;
 151
 152     e1      = exp(-x/a[1]);
 153     e2      = exp(-x/a[3]);
 154     *y      = a[2]*e1 + (1-a[2])*e2;
 155     dyda[1] = x*a[2]*e1/sqr(a[1]);
 156     dyda[2] = e1-e2;
 157     dyda[3] = x*(1-a[2])*e2/sqr(a[3]);
 158 }
 159
 160 static void exp_5_parm(real x, real a[], real *y, real dyda[])
 161 {
 162     /* Fit to function
 163      *
 164      * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5
 165      *
 166      */
 167
 168     real e1, e2;
 169
 170     e1      = exp(-x/a[2]);
 171     e2      = exp(-x/a[4]);
 172     *y      = a[1]*e1 + a[3]*e2 + a[5];
 173
 174     if (debug)
 175     {
 176         fprintf(debug, "exp_5_parm called: x = %10.3f  y = %10.3f\n"
 177                 "a = ( %8.3f  %8.3f  %8.3f  %8.3f  %8.3f)\n",
 178                 x, *y, a[1], a[2], a[3], a[4], a[5]);
 179     }
 180     dyda[1] = e1;
 181     dyda[2] = x*e1/sqr(a[2]);
 182     dyda[3] = e2;
 183     dyda[4] = x*e2/sqr(a[4]);
 184     dyda[5] = 0;
 185 }
 186
 187 static void exp_7_parm(real x, real a[], real *y, real dyda[])
 188 {
 189     /* Fit to function
 190      *
 191      * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
 192      *
 193      */
 194
 195     real e1, e2, e3;
 196
 197     e1      = exp(-x/a[2]);
 198     e2      = exp(-x/a[4]);
 199     e3      = exp(-x/a[6]);
 200     *y      = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7];
 201
 202     dyda[1] = e1;
 203     dyda[2] = x*e1/sqr(a[2]);
 204     dyda[3] = e2;
 205     dyda[4] = x*e2/sqr(a[4]);
 206     dyda[5] = e3;
 207     dyda[6] = x*e3/sqr(a[6]);
 208     dyda[7] = 0;
 209 }
 210
 211 static void exp_9_parm(real x, real a[], real *y, real dyda[])
 212 {
 213     /* Fit to function
 214      *
 215      * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
 216      *
 217      */
 218
 219     real e1, e2, e3, e4;
 220
 221     e1      = exp(-x/a[2]);
 222     e2      = exp(-x/a[4]);
 223     e3      = exp(-x/a[6]);
 224     e4      = exp(-x/a[8]);
 225     *y      = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7]*e4 + a[9];
 226
 227     dyda[1] = e1;
 228     dyda[2] = x*e1/sqr(a[2]);
 229     dyda[3] = e2;
 230     dyda[4] = x*e2/sqr(a[4]);
 231     dyda[5] = e3;
 232     dyda[6] = x*e3/sqr(a[6]);
 233     dyda[7] = e4;
 234     dyda[8] = x*e4/sqr(a[8]);
 235     dyda[9] = 0;
 236 }
 237
 238 static void vac_2_parm(real x, real a[], real *y, real dyda[])
 239 {
 240     /* Fit to function
 241      *
 242      * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
 243      *
 244      *   = exp(-v) (cosh(wv) + 1/w sinh(wv))
 245      *
 246      *    v = x/(2 a1)
 247      *    w = sqrt(1 - a2)
 248      *
 249      *    For tranverse current autocorrelation functions:
 250      *       a1 = tau
 251      *       a2 = 4 tau (eta/rho) k^2
 252      *
 253      */
 254
 255     double v, det, omega, omega2, em, ec, es;
 256
 257     v   = x/(2*a[1]);
 258     det = 1 - a[2];
 259     em  = exp(-v);
 260     if (det != 0)
 261     {
 262         omega2  = fabs(det);
 263         omega   = sqrt(omega2);
 264         if (det > 0)
 265         {
 266             ec = em*0.5*(exp(omega*v)+exp(-omega*v));
 267             es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
 268         }
 269         else
 270         {
 271             ec = em*cos(omega*v);
 272             es = em*sin(omega*v)/omega;
 273         }
 274         *y      = ec + es;
 275         dyda[2] = (v/det*ec+(v-1/det)*es)/(-2.0);
 276         dyda[1] = (1-det)*v/a[1]*es;
 277     }
 278     else
 279     {
 280         *y      = (1+v)*em;
 281         dyda[2] = -v*v*em*(0.5+v/6);
 282         dyda[1] = v*v/a[1]*em;
 283     }
 284 }
 285
 286 static void errest_3_parm(real x, real a[], real *y, real dyda[])
 287 {
 288     real e1, e2, v1, v2;
 289
 290     if (a[1])
 291     {
 292         e1 = exp(-x/a[1]) - 1;
 293     }
 294     else
 295     {
 296         e1 = 0;
 297     }
 298     if (a[3])
 299     {
 300         e2 = exp(-x/a[3]) - 1;
 301     }
 302     else
 303     {
 304         e2 = 0;
 305     }
 306
 307     if (x > 0)
 308     {
 309         v1      = 2*a[1]*(e1*a[1]/x + 1);
 310         v2      = 2*a[3]*(e2*a[3]/x + 1);
 311         *y      = a[2]*v1 + (1-a[2])*v2;
 312         dyda[1] = 2*     a[2] *(v1/a[1] + e1);
 313         dyda[3] = 2*(1 - a[2])*(v2/a[3] + e2);
 314         dyda[2] = (v1 - v2);
 315     }
 316     else
 317     {
 318         *y      = 0;
 319         dyda[1] = 0;
 320         dyda[3] = 0;
 321         dyda[2] = 0;
 322     }
 323 }
 324
 325 typedef void (*myfitfn)(real x, real a[], real *y, real dyda[]);
 326 myfitfn mfitfn[effnNR] = {
 327     exp_one_parm, exp_one_parm, exp_two_parm, exp_3_parm, vac_2_parm,
 328     exp_5_parm,   exp_7_parm,   exp_9_parm, erffit,  errest_3_parm
 329 };
 330
 331 real fit_function(int eFitFn, real *parm, real x)
 332 {
 333     static real y, dum[8];
 334
 335     mfitfn[eFitFn](x, parm-1, &y, dum);
 336
 337     return y;
 338 }
 339
 340 /* lmfit_exp supports up to 3 parameter fitting of exponential functions */
 341 static gmx_bool lmfit_exp(int nfit, real x[], real y[], real dy[], real ftol,
 342                           real parm[], real dparm[], gmx_bool bVerbose,
 343                           int eFitFn, int fix)
 344 {
 345     real     chisq, ochisq, alamda;
 346     real    *a, **covar, **alpha, *dum;
 347     gmx_bool bCont;
 348     int      i, j, ma, mfit, *lista, *ia;
 349
 350     if ((eFitFn < 0) || (eFitFn >= effnNR))
 351     {
 352         gmx_fatal(FARGS, "fitfn = %d, should be in 0..%d (%s,%d)",
 353                   effnNR-1, eFitFn, __FILE__, __LINE__);
 354     }
 355
 356     ma = mfit = nfp_ffn[eFitFn];   /* number of parameters to fit */
 357     snew(a, ma+1);
 358     snew(covar, ma+1);
 359     snew(alpha, ma+1);
 360     snew(lista, ma+1);
 361     snew(ia, ma+1);
 362     snew(dum, ma+1);
 363     for (i = 1; (i < ma+1); i++)
 364     {
 365         lista[i] = i;
 366         ia[i]    = 1; /* fixed bug B.S.S 19/11  */
 367         snew(covar[i], ma+1);
 368         snew(alpha[i], ma+1);
 369     }
 370     if (fix)
 371     {
 372         if (bVerbose)
 373         {
 374             fprintf(stderr, "Will keep parameters fixed during fit procedure: %d\n",
 375                     fix);
 376         }
 377         for (i = 0; i < ma; i++)
 378         {
 379             if (fix & 1<<i)
 380             {
 381                 ia[i+1] = 0;
 382             }
 383         }
 384     }
 385     if (debug)
 386     {
 387         fprintf(debug, "%d parameter fit\n", mfit);
 388     }
 389
 390     /* Initial params */
 391     alamda = -1;  /* Starting value   */
 392     chisq  = 1e12;
 393     for (i = 0; (i < mfit); i++)
 394     {
 395         a[i+1] = parm[i];
 396     }
 397
 398     j = 0;
 399     if (bVerbose)
 400     {
 401         fprintf(stderr, "%4s  %10s  %10s  %10s  %10s  %10s %10s\n",
 402                 "Step", "chi^2", "Lambda", "A1", "A2", "A3", "A4");
 403     }
 404     do
 405     {
 406         ochisq = chisq;
 407         /* mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
 408          *   &chisq,expfn[mfit-1],&alamda)
 409          */
 410         if (!mrqmin_new(x-1, y-1, dy-1, nfit, a, ia, ma, covar, alpha, &chisq,
 411                         mfitfn[eFitFn], &alamda))
 412         {
 413             return FALSE;
 414         }
 415
 416         if (bVerbose)
 417         {
 418             fprintf(stderr, "%4d  %10.5e  %10.5e  %10.5e",
 419                     j, chisq, alamda, a[1]);
 420             if (mfit > 1)
 421             {
 422                 fprintf(stderr, "  %10.5e", a[2]);
 423             }
 424             if (mfit > 2)
 425             {
 426                 fprintf(stderr, "  %10.5e", a[3]);
 427             }
 428             if (mfit > 3)
 429             {
 430                 fprintf(stderr, " %10.5e", a[4]);
 431             }
 432             fprintf(stderr, "\n");
 433         }
 434         j++;
 435         bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
 436                  ((ochisq == chisq)));
 437     }
 438     while (bCont && (alamda != 0.0) && (j < 50));
 439     if (bVerbose)
 440     {
 441         fprintf(stderr, "\n");
 442     }
 443
 444     /* Now get the covariance matrix out */
 445     alamda = 0;
 446
 447     /*  mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
 448      * &chisq,expfn[mfit-1],&alamda)
 449      */
 450     if (!mrqmin_new(x-1, y-1, dy-1, nfit, a, ia, ma, covar, alpha, &chisq,
 451                     mfitfn[eFitFn], &alamda))
 452     {
 453         return FALSE;
 454     }
 455
 456     for (j = 0; (j < mfit); j++)
 457     {
 458         parm[j]  = a[j+1];
 459         dparm[j] = covar[j+1][j+1];
 460     }
 461
 462     for (i = 0; (i < ma+1); i++)
 463     {
 464         sfree(covar[i]);
 465         sfree(alpha[i]);
 466     }
 467     sfree(a);
 468     sfree(covar);
 469     sfree(alpha);
 470     sfree(lista);
 471     sfree(dum);
 472
 473     return TRUE;
 474 }
 475
 476 real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[],
 477               real begintimefit, real endtimefit, const output_env_t oenv,
 478               gmx_bool bVerbose, int eFitFn, real fitparms[], int fix)
 479 {
 480     FILE *fp;
 481     char  buf[32];
 482
 483     int   i, j, nparm, nfitpnts;
 484     real  integral, ttt;
 485     real *parm, *dparm;
 486     real *x, *y, *dy;
 487     real  ftol = 1e-4;
 488
 489     nparm = nfp_ffn[eFitFn];
 490     if (debug)
 491     {
 492         fprintf(debug, "There are %d points to fit %d vars!\n", ndata, nparm);
 493         fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n",
 494                 eFitFn, begintimefit, endtimefit, dt);
 495     }
 496
 497     snew(x, ndata);
 498     snew(y, ndata);
 499     snew(dy, ndata);
 500
 501     j = 0;
 502     for (i = 0; i < ndata; i++)
 503     {
 504         ttt = x0 ? x0[i] : dt*i;
 505         if (ttt >= begintimefit && ttt <= endtimefit)
 506         {
 507             x[j] = ttt;
 508             y[j] = c1[i];
 509
 510             /* mrqmin does not like sig to be zero */
 511             if (sig[i] < 1.0e-7)
 512             {
 513                 dy[j] = 1.0e-7;
 514             }
 515             else
 516             {
 517                 dy[j] = sig[i];
 518             }
 519             if (debug)
 520             {
 521                 fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy= %g\n",
 522                         j, i, x[j], y[j], dy[j]);
 523             }
 524             j++;
 525         }
 526     }
 527     nfitpnts = j;
 528     integral = 0;
 529     if (nfitpnts < nparm)
 530     {
 531         fprintf(stderr, "Not enough data points for fitting!\n");
 532     }
 533     else
 534     {
 535         snew(parm, nparm);
 536         snew(dparm, nparm);
 537         if (fitparms)
 538         {
 539             for (i = 0; (i < nparm); i++)
 540             {
 541                 parm[i] = fitparms[i];
 542             }
 543         }
 544
 545         if (!lmfit_exp(nfitpnts, x, y, dy, ftol, parm, dparm, bVerbose, eFitFn, fix))
 546         {
 547             fprintf(stderr, "Fit failed!\n");
 548         }
 549         else if (nparm <= 3)
 550         {
 551             /* Compute the integral from begintimefit */
 552             if (nparm == 3)
 553             {
 554                 integral = (parm[0]*myexp(begintimefit, parm[1],  parm[0]) +
 555                             parm[2]*myexp(begintimefit, 1-parm[1], parm[2]));
 556             }
 557             else if (nparm == 2)
 558             {
 559                 integral = parm[0]*myexp(begintimefit, parm[1],  parm[0]);
 560             }
 561             else if (nparm == 1)
 562             {
 563                 integral = parm[0]*myexp(begintimefit, 1,  parm[0]);
 564             }
 565             else
 566             {
 567                 gmx_fatal(FARGS, "nparm = %d in file %s, line %d",
 568                           nparm, __FILE__, __LINE__);
 569             }
 570
 571             /* Generate THE output */
 572             if (bVerbose)
 573             {
 574                 fprintf(stderr, "FIT: # points used in fit is: %d\n", nfitpnts);
 575                 fprintf(stderr, "FIT: %21s%21s%21s\n",
 576                         "parm0     ", "parm1 (ps)   ", "parm2 (ps)    ");
 577                 fprintf(stderr, "FIT: ------------------------------------------------------------\n");
 578                 fprintf(stderr, "FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
 579                         parm[0], dparm[0], parm[1], dparm[1], parm[2], dparm[2]);
 580                 fprintf(stderr, "FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
 581                         begintimefit, integral);
 582
 583                 sprintf(buf, "test%d.xvg", nfitpnts);
 584                 fp = xvgropen(buf, "C(t) + Fit to C(t)", "Time (ps)", "C(t)", oenv);
 585                 fprintf(fp, "# parm0 = %g, parm1 = %g, parm2 = %g\n",
 586                         parm[0], parm[1], parm[2]);
 587                 for (j = 0; (j < nfitpnts); j++)
 588                 {
 589                     ttt = x0 ? x0[j] : dt*j;
 590                     fprintf(fp, "%10.5e  %10.5e  %10.5e\n",
 591                             ttt, c1[j], fit_function(eFitFn, parm, ttt));
 592                 }
 593                 xvgrclose(fp);
 594             }
 595         }
 596         if (fitparms)
 597         {
 598             for (i = 0; (i < nparm); i++)
 599             {
 600                 fitparms[i] = parm[i];
 601             }
 602         }
 603         sfree(parm);
 604         sfree(dparm);
 605     }
 606
 607     sfree(x);
 608     sfree(y);
 609     sfree(dy);
 610
 611     return integral;
 612 }
 613
 614 void do_expfit(int ndata, real c1[], real dt, real begintimefit, real endtimefit)
 615 {
 616     int   i, n;
 617     real *x, *y, *Dy;
 618     real  aa, bb, saa, sbb, A, tau, dA, dtau;
 619
 620     fprintf(stderr, "Will fit data from %g (ps) to %g (ps).\n",
 621             begintimefit, endtimefit);
 622
 623     snew(x, ndata); /* allocate the maximum necessary space */
 624     snew(y, ndata);
 625     snew(Dy, ndata);
 626     n = 0;
 627
 628     for (i = 0; (i < ndata); i++)
 629     {
 630         if ( (dt*i >= begintimefit) && (dt*i <= endtimefit) )
 631         {
 632             x[n]  = dt*i;
 633             y[n]  = c1[i];
 634             Dy[n] = 0.5;
 635             fprintf(stderr, "n= %d, i= %d, x= %g, y= %g\n", n, i, x[n], y[n]);
 636             n++;
 637         }
 638     }
 639     fprintf(stderr, "# of data points used in the fit is : %d\n\n", n);
 640     expfit(n, x, y, Dy, &aa, &saa, &bb, &sbb);
 641
 642     A    = exp(aa);
 643     dA   = exp(aa)*saa;
 644     tau  = -1.0/bb;
 645     dtau = sbb/sqr(bb);
 646     fprintf(stderr, "Fitted to y=exp(a+bx):\n");
 647     fprintf(stderr, "a = %10.5f\t b = %10.5f", aa, bb);
 648     fprintf(stderr, "\n");
 649     fprintf(stderr, "Fitted to y=Aexp(-x/tau):\n");
 650     fprintf(stderr, "A  = %10.5f\t tau  = %10.5f\n", A, tau);
 651     fprintf(stderr, "dA = %10.5f\t dtau = %10.5f\n", dA, dtau);
 652 }
 653
 654
 655 void expfit(int n, real *x, real *y, real *Dy, real *a, real *sa,
 656             real *b, real *sb)
 657 {
 658     real *w, *ly, A, SA, B, SB;
 659     int   i;
 660     real  sum, xbar, ybar, Sxx, Sxy, wr2, chi2, gamma, Db;
 661
 662 #define ZERO 0.0
 663 #define ONE 1.0
 664 #define ONEP5 1.5
 665 #define TWO 2.0
 666
 667 #define sqr(x) ((x)*(x))
 668
 669     /*allocate memory */
 670     snew(w, n);
 671     snew(ly, n);
 672
 673     /* Calculate weights and values of ln(y). */
 674     for (i = 0; (i < n); i++)
 675     {
 676         w[i]  = sqr(y[i]/Dy[i]);
 677         ly[i] = log(y[i]);
 678     }
 679
 680     /* The weighted averages of x and y: xbar and ybar. */
 681     sum  = ZERO;
 682     xbar = ZERO;
 683     ybar = ZERO;
 684     for (i = 0; (i < n); i++)
 685     {
 686         sum  += w[i];
 687         xbar += w[i]*x[i];
 688         ybar += w[i]*ly[i];
 689     }
 690     xbar /= sum;
 691     ybar /= sum;
 692
 693     /* The centered product sums Sxx and Sxy, and hence A and B. */
 694     Sxx = ZERO;
 695     Sxy = ZERO;
 696     for (i = 0; (i < n); i++)
 697     {
 698         Sxx += w[i]*sqr(x[i]-xbar);
 699         Sxy += w[i]*(x[i]-xbar)*(ly[i]-ybar);
 700     }
 701     B = Sxy/Sxx;
 702     A = ybar-B*xbar;
 703
 704     /* Chi-squared (chi2) and gamma. */
 705     chi2  = ZERO;
 706     gamma = ZERO;
 707     for (i = 0; (i < n); i++)
 708     {
 709         wr2    = w[i]*sqr(ly[i]-A-B*x[i]);
 710         chi2  += wr2;
 711         gamma += wr2*(x[i]-xbar);
 712     }
 713
 714     /* Refined values of A and B. Also SA and SB. */
 715     Db  = -ONEP5*gamma/Sxx;
 716     B  += Db;
 717     A  -= ONEP5*chi2/sum-xbar*Db;
 718     SB  = sqrt((chi2/(n-2))/Sxx);
 719     SA  = SB*sqrt(Sxx/sum+sqr(xbar));
 720     *a  = A;
 721     *b  = B;
 722     *sa = SA;
 723     *sb = SB;
 724 }