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