unit expfuncUnit;
{
 calculates exponential functions
 of type:
         1 : y = ab^x+c
         2 : y = ae^(bx)+c
         3 : y = e^(ax^2+bx+c)
         4: ab^x
         5: ae^(bx)
}

interface

uses sysutils,xlateUnit,puntenUnit,formUnit,scaleUnit,dialogs;

type puntarray = array[1..40] of TPunt; //punt = point , coordinates of points

procedure expfunc(n : byte); //n= exp function type requested

implementation

uses controlUnit,klkwUnit;

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

var Epoints : puntarray;
    a,b,c : double;       //y = a*e(bx)+c
    ePower : boolean;     //use power of e
    dx,dy : array[1..40] of double;
    dxy : array[1..40] of TPunt;

procedure ExpPlot(s : string);
begin
 controlform.formedit.text := s;
 addform(formNr);
end;

procedure Getpoints;
//sort points by x value, ascending order
var i,j : byte;
    h : Tpunt;
begin
 for i := 1 to maxpunt do     //1 to max. point , maxpunt in other unit
  begin
   Epoints[i].x := punten[i].x; //punten = points, from other unit
   Epoints[i].y := punten[i].y;
  end;
 for i := 1 to maxpunt-1 do
  for j := i+1 to maxpunt do
    if Epoints[j].x < Epoints[i].x then
     begin
      h.x := Epoints[i].x;
      h.y := Epoints[i].y;
      Epoints[i].x := Epoints[j].x;
      Epoints[i].y := Epoints[j].y;
      Epoints[j].x := h.x;
      Epoints[j].y := h.y;
     end;
end;

function InsufficientPoints(n : byte) : boolean;
//error if < n points supplied
begin
 if maxpunt < n then
  begin
   controlform.msgtext.caption := 'minimal '+inttostr(n)+' points needed';
   result := true;
  end
 else result := false;
end;

procedure TwoPointsExp;
//calculate a and b, if 2 points supplied
var sign1,sign2 : boolean;
    s : string;
    i : byte;
begin
 sign1 := Epoints[1].y < 0;
 sign2 := Epoints[2].y < 0;
 if sign1 <> sign2 then
  begin
   controlform.msgtext.caption := 'no exp. function possible: signs unlike';
   exit;
  end;
 if sign1 then
  for i := 1 to 2 do Epoints[i].y := -Epoints[i].y;//a negative, invert
 for i := 1 to 2 do Epoints[i].y := ln(Epoints[i].y);
 b := (Epoints[2].y - Epoints[1].y)/(Epoints[2].x - Epoints[1].x);
 a := Epoints[1].y - b*Epoints[1].x;
 a := exp(a);
 if not Epower then b := exp(b);
 if sign1 then a := -a;
//a,b calculated
 s := 'y = ';
 s := s + formatfloat(flstring,a) + ' * ';
 if Epower then s := s + 'e^(' + formatfloat(Flstring,b) +'x)'
  else s := s + formatfloat(Flstring,b) + '^x';
 expPlot(s);  //to other unit to be plotted
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
 avG := 0;
 for i := 1 to maxpunt-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[maxpunt].x-Epoints[1].x));//exp added
 ee := 0;
 for i := 1 to maxpunt-1 do
  ee := ee + abs(avG-growthFactor[i]);
 result :=  ee;
end;

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

procedure ExpAbort(m : byte);
const abortmessage : array[1..2] of string =
      ('no exp. function',
       'numbers out of reach');
begin
 controlform.msgtext.caption := abortmessage[m];
 controlform.formedit.text := '';
end;

procedure morepointsExp;
var Sum,deltaC,Bvar1,Bvar2 : double;
    i : byte;
    w : word;
    s : string;
    abort : boolean;
begin
 makeXYarrays;
//approximate value of b ....growthfactor e^b
 try
  sum := 0;
  for i := 1 to maxpunt-2 do
   Sum := Sum + ln(dxy[i+1].ydxy[i].y)/(dxy[i+1].x-dxy[i].x);
  b := Sum / (maxpunt - 2);
//approximate value of a
  Sum := 0;
  for i := 1 to maxpunt-1 do
   Sum := Sum + dy[i]/(exp(b*epoints[i+1].x) - exp(b*epoints[i].x));
  a := Sum / (maxpunt - 1);
//approximate value of c
  sum := 0;
  for i := 1 to maxpunt do
   sum := sum + epoints[i].y-a*exp(b*epoints[i].x);
  c := sum / maxpunt;
 except
  expAbort(1);
  exit;
 end;
 if a > 0 then deltaC := -scaledY //move away from asymptote first
  else deltaC := scaledY;         //scaledY = vertical distance between 2 pixels
 repeat                           //..      defined in other unit
  Bvar1 := Deviation(c,abort);
  if abort then c := c + deltaC;
 until not abort;
 w := 0;
 deltaC := -deltaC;//towards asymptote first
 repeat
  inc(w);
  Bvar2 := Deviation(c + deltaC,abort);
  if abort then deltaC := deltaC2  //half distance in same direction
  else
   if Bvar2 < Bvar1 then
    begin
     Bvar1 := Bvar2;  //if better result
     c := c + deltaC;
    end
   else deltaC := -deltaC2;
 until (abs(deltaC) <= 1e-8) or (w > 100);
 if abort then
  begin
   expAbort(1);
   exit;
  end;
 if w > 100 then
  begin
   expAbort(2);
   exit;
  end;  
//call least square procedure to make best a,b
 for i := 1 to maxpunt do
  epoints[i].y := ln(abs(epoints[i].y-c));
 buildexp12(epoints);  //this is call to least square unit
 if a < 0 then a := -exp(epoints[1].y)
 else a := exp(epoints[1].y);
 b := epoints[2].y;
//make formula
 s := 'y = ' + formatfloat(Flstring,a)+' * ';
 if Epower then s := s + 'e^(' + formatfloat(Flstring,b) + 'x)'
 else begin
       b := exp(b);
       s := s + formatfloat(Flstring,b) + '^x';
      end;
 if c < 0 then begin
                c := abs(c);
                s := s + ' - ';
               end
 else s := s + ' + ';              
 s := s + formatfloat(Flstring,c);
 expPlot(s);
end;

procedure makeExpFunction;
begin
 if Insufficientpoints(2) then exit;//msg al gepost
 GetPoints;
 if maxpunt = 2 then
  begin
   c := 0;
   TwoPointsExp;
   exit;
  end;
 morePointsExp;
end;

procedure exp3;
var i : byte;
    s : string;
begin
 if Insufficientpoints(3) then exit;
 getPoints;
 for i := 1 to maxpunt do
  if epoints[i].y <= 0 then
   begin
    controlform.msgtext.caption := 'fout: y <= 0';
    exit;
   end
   else epoints[i].y := ln(epoints[i].y);
 buildexp3(epoints); //call least squares unit,return a,b,c at epoints[1..3].y
 s := 'y = e^(';
 s := s + formatfloat(Flstring,epoints[3].y) +'x^2';
 if epoints[2].y < 0 then
  begin
   s := s + ' - ';
   epoints[2].y := abs(epoints[2].y);
  end
  else s := s + ' + ';
 s := s + formatfloat(Flstring,epoints[2].y) +'x';
 if epoints[1].y < 0 then
  begin
   s := s + ' - ';
   epoints[1].y := abs(epoints[1].y);
  end
  else s := s + ' + ';
 s := s + formatfloat('####0.0#####',epoints[1].y);
 s := s + ')';
 ExpPlot(s);
end;

procedure exp45;
var i : byte;
    s : string;
begin
 if InsufficientPoints(2) then exit;
 getPoints;
 for i := 1 to maxpunt do
  if epoints[i].y <= 0 then
   begin
    controlform.msgtext.caption := 'error: y <= 0';
    exit;
   end
   else epoints[i].y := ln(epoints[i].y);
 buildexp12(epoints);//call least squares unit, return a,b at epoints[1,2].y
 a := exp(epoints[1].y);
 b := epoints[2].y;
 if not Epower then b := exp(b);
//a,b,c calculated
 s := 'y = ';
 s := s + formatfloat(flstring,a) + ' * ';
 if Epower then s := s + 'e^(' + formatfloat(Flstring,b) +'x)'
  else s := s + formatfloat(Flstring,b) + '^x';
 expPlot(s);  //call to plot the function, outside this unit
end;

//--- central call to make exponential function of type n

procedure expFunc(n : byte);
begin
 case n of
  1 : begin
       Epower := false; //no power of e
       makeExpfunction;
      end;
  2 : begin
       Epower := true; //is power of e
       makeExpfunction;
      end;
  3 : exp3;            //build function type 3
  4 : begin
       Epower := false;
       exp45;          //build function type 4 or 5
      end;
  5 : begin
       Epower := true;
       exp45;          //build function type 4 or 5
      end;
 end;//case
end;

end.