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