unit Unit2;
{
 function  y = ae^(bx)+c
 calculate a,b,c from [x1,y1] , [x2,y2] , ...., [xn,yn]
}

interface

uses sysutils;

type  Tpnt = record
               x,y : double;
              end;
      Tpoints = array[1..10] of Tpnt;

procedure makeExpFunction;

procedure GetParams(var pa,pb,pc : double);
function getXY(nr : byte; var pdx,pdy, pcx, pdq : double) : boolean;

implementation

uses unit1,unit3,unit4;

const FlString : string = '####0.0#####';
      scaleDy = 1;

var Epoints : TPoints;
    a,b,c : double;                       //y = a*e(bx)+c
    dx,dy : array[1..10] of double;
    cx    : array[1..10] of double;       //center between x values
    dq    : array[1..10] of double;       //differential quotient
    maxpnt : byte;                        //number of points

//testmode helpers

procedure GetParams(var pa,pb,pc : double);
begin
 pa := a; pb := b; pc := c;
end;

function getXY(nr : byte; var pdx,pdy, pcx, pdq : double) : boolean;
begin
 result := nr < maxpnt;
 if result then
  begin
   pdx := dx[nr];
   pdy := dy[nr];
   pcx := cx[nr];
   pdq := dq[nr];
  end;
end;

//----

procedure makeXYarrays;
var i : byte;
begin
 for i := 1 to maxpnt-1 do
  begin
   dx[i] := epoints[i+1].x - epoints[i].x;
   dy[i] := epoints[i+1].y - epoints[i].y;
   cx[i] := (epoints[i].x + epoints[i+1].x)/2;
   dq[i] := dy[i]/dx[i];
  end;
end;    

function Deviation(cc : double; var af : boolean) : double;
//calculate variation in growth factors after c asymptote correction
//called for fine tuning c in ae^(bx)+c formula
//cc: y value (asymptote); af : abort flag...abort if intersection
var i : byte;
    avG,ee : double;
    growthfactor : array[1..40] of double;
begin
 result := 0;
 avG := 0;
 for i := 1 to maxpnt-1 do
  begin
   ee := (Epoints[i].y-cc);//denominator
   if ee <> 0 then ee := (Epoints[i+1].y-cc)/ee;
   if ee <= 0 then //test for intersection of asymptote and function
    begin
     af := true;//abort flag
     exit;
    end;
   ee := ln(ee);
   avG := avG + ee;
   growthfactor[i] := exp(ee/dx[i]);//save
  end;
 af := false;
 avG := exp(avG /(Epoints[maxpnt].x-Epoints[1].x));//exp added
 ee := 0;
 for i := 1 to maxpnt-1 do
  ee := ee + abs(avG-growthFactor[i]);
 result :=  ee;
end;

procedure morepointsExp(var code : byte);
//called for >= 3 points
var Sum,deltaC,Bvar1,Bvar2 : double;
    i : byte;
    w : word;
    abort : boolean;
begin
//approximate value of b ....growthfactor e^b
 try
  sum := 0;
  for i := 1 to maxpnt-2 do
  Sum := Sum + ln(dq[i+1]/dq[i])/(cx[i+1]-cx[i]);
  b := Sum / (maxpnt - 2);
  if testmode then testreport(2);
//approximate value of a
  Sum := 0;
  for i := 1 to maxpnt-1 do
   Sum := Sum + dy[i]/(exp(b*epoints[i+1].x) - exp(b*epoints[i].x));
  a := Sum / (maxpnt - 1);
  if testmode then testreport(3);
//approximate value of c
  sum := 0;
  for i := 1 to maxpnt do
   sum := sum + epoints[i].y-a*exp(b*epoints[i].x);
  c := sum / maxpnt;
  if testmode then testreport(4);
 except
  code := 1;
  exit;
 end;
 if a > 0 then deltaC := -scaledY //move away from asymptote first
  else deltaC := scaledY;
 if testmode then
  begin
   testreport(5);
   testreportDC(deltaC);
  end;
 repeat
  Bvar1 := Deviation(c,abort);
  if abort then c := c + deltaC;
  if testmode then testreportDCA(Bvar1,c,deltaC,abort);
 until not abort;
 w := 0;
 deltaC := -deltaC;//towards asymptote first
 if testmode then testreport(6);
 repeat
  inc(w);
  Bvar2 := Deviation(c + deltaC,abort);
  if abort then deltaC := deltaC/2  //half distance in same direction
  else
   if Bvar2 < Bvar1 then
    begin
     Bvar1 := Bvar2;  //if better result
     c := c + deltaC;
    end
   else deltaC := -deltaC/2;
  if testmode then testreportDCA(Bvar1,c,deltaC,abort);
 until (abs(deltaC) <= 1e-8) or (w > 100);
 if abort then
  begin
   code := 1;
   exit;
  end;
 if w > 100 then
  begin
   code := 2;
   exit;
  end;
 if testmode then testreport(7);
//call least square procedure to make best a,b
 for i := 1 to maxpnt do epoints[i].y := ln(abs(epoints[i].y-c));
 buildexp12(epoints,maxpnt);  //call least square unit
 if a < 0 then a := -exp(epoints[1].y)
  else a := exp(epoints[1].y);
 b := epoints[2].y;
 if testmode then testreport(8);
end;

//--- main call ---

procedure makeExpFunction;
var i,code : byte; //0:OK  1,2: error
    s : string;
begin
 maxpnt := pcount;
 for i := 1 to maxpnt do            //get points
  begin
   Epoints[i].x := points[i].x;
   Epoints[i].y := points[i].y;
  end;
//
 makeXYarrays;
 if testmode then testreport(1);
 code := 0;
 morePointsExp(code);
 case code of
  0 : begin
       aa := a; bb := exp(b); cc := c;     //ae^(bx)+c ....>   ab^x + c
       setparamsOK(true);
       s := 'a,b,c calculated';
      end;
  1 : begin
       setparamsOK(false);
       s := 'no exponential function';
      end;
  2 : begin
       setparamsOK(false);
       s := 'points out of range';
      end;
   else s := '';
 end;
 form1.msgtext.Caption := s;
end;

end.