home *** CD-ROM | disk | FTP | other *** search
/ gondwana.ecr.mu.oz.au/pub/ / Graphics.tar / Graphics / vopl.tar.Z / vopl.tar / vopl / src / ilinfit.c < prev    next >
C/C++ Source or Header  |  1989-09-11  |  2KB  |  100 lines
  1. #include "vopl.h"
  2. #include "math.h"
  3.  
  4. /*
  5.  * pefit
  6.  *
  7.  *    power equation fit    
  8.  *
  9.  *    ln y = ln a + b.ln x
  10.  *
  11.  */
  12. void
  13. pefit(x, y, n)
  14.     float    x[], y[];
  15.     int    n;
  16. {
  17.     float    xval, xstep, yval, sumx, sumy, sumxy, sumxx, a, b, yave, xave;
  18.     int    i;
  19.     
  20.     sumx = sumy = sumxy = sumxx = 0.0;
  21.  
  22.     for (i = plotdev.startind; i < n; i += plotdev.arrayind) {
  23.         sumx += WhatX(log(x[i]));
  24.         sumy += WhatY(log(y[i]));
  25.         sumxy += WhatX(log(x[i])) * WhatY(log(y[i]));
  26.         sumxx += WhatX(log(x[i])) * WhatY(log(x[i]));
  27.     }
  28.  
  29.     yave = sumy / n;
  30.     xave = sumx / n;
  31.  
  32.     b = (n * sumxy - sumx * sumy) / (n * sumxx - sumx * sumx);
  33.     a = exp(yave - b * xave);
  34.  
  35.     xval = WhatX(plotdev.axes[XIND].min);
  36.     yave = a * pow(xval, b);
  37.  
  38.     move2(xval, yave);
  39.  
  40.     xstep = (WhatX(plotdev.axes[XIND].max) - WhatX(plotdev.axes[XIND].min)) / (2 * n);
  41.     for (i = plotdev.startind; i < 2 * n; i += plotdev.arrayind) {
  42.         xval += xstep;
  43.         yval = a * pow(xval, b); 
  44.         draw2(xval, yval);
  45.     }
  46. }
  47.  
  48. /*
  49.  * sgrfit
  50.  *
  51.  *    saturated growth rate fit
  52.  *
  53.  *    y = c * x / (d + x)
  54.  *
  55.  *      invert 
  56.  *
  57.  *      1/y = d/c . 1/x + 1/c 
  58.  *
  59.  *
  60.  */
  61. void
  62. sgrfit(x, y, n)
  63.     float    x[], y[];
  64.     int    n;
  65. {
  66.     float    xval, xstep, yval, sumx, sumy, sumxy, sumxx, a, b, yave, xave;
  67.     float    a3, b3;
  68.     int    i;
  69.     
  70.     sumx = sumy = sumxy = sumxx = 0.0;
  71.  
  72.     for (i = plotdev.startind; i < n; i += plotdev.arrayind) {
  73.         sumx += 1.0 / x[i];
  74.         sumy += 1.0 / y[i];
  75.         sumxy += (1.0 / x[i]) * (1.0 / y[i]);
  76.         sumxx += (1.0 / x[i]) * (1.0 / x[i]);
  77.     }
  78.  
  79.     yave = sumy / n;
  80.     xave = sumx / n;
  81.  
  82.     b = (n * sumxy - sumx * sumy) / (n * sumxx - sumx * sumx);
  83.     a = yave - b * xave;
  84.  
  85.     a3 = 1.0 / a;
  86.     b3 = b * a3; 
  87.  
  88.     xval = WhatX(plotdev.axes[XIND].min);
  89.     yval = a3 * xval / (b3 + xval);
  90.     move2(xval, yave);
  91.  
  92.     xstep = (WhatX(plotdev.axes[XIND].max) - WhatX(plotdev.axes[XIND].min)) / (2 * n);
  93.  
  94.     for (i = plotdev.startind; i < 2 * n; i += plotdev.arrayind) {
  95.         xval += xstep;
  96.         yval = a3 * xval / (b3 + xval);
  97.         draw2(xval, yval);
  98.     }
  99. }
  100.