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