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