Unit Rekenen;

INTERFACE

Uses Math, tabellen;
 
type
  punt = record
    x,y : double;
  end;
  punten = array of punt;

procedure uitkijken(jaar,licht : lijst; var A,B : double);
procedure kijkuit(jaar,licht : lijst; var A,B,R : double);
procedure NewtonHdB(jaar,licht : lijst; var A,B : double; var eps : punt);
procedure Leibovici(jaar,licht : lijst; var A,B : double; var eps : double);
procedure Jacquelin(jaar,licht : lijst; var A,B : double);
procedure Raaklijn(jaar,licht,w : lijst; var A,B,dA : double);

IMPLEMENTATION

procedure uitkijken(jaar,licht : lijst; var A,B : double);
{
  Least Squares best fit
  to hyperbola y=B/(x-A)
  with linear regression
}
var
  N,k : integer;
  D : double;
  xx,yy : array of double;
  M11,M22,M12,R1,R2 : double;
begin
  N := Length(jaar);
  SetLength(xx,N);
  SetLength(yy,N);
  for k := 0 to N-1 do
  begin
    xx[k] := jaar[k];
    yy[k] := licht[k];
  end;
  M11 := 0; M12 := 0; M22 := 0;
  R1 := 0; R2 := 0;
  for k := 0 to N-1 do
  begin
    if yy[k] <= 0 then Continue;
    M11 := M11 + yy[k]*yy[k]; 
    M12 := M12 + yy[k];
    M22 := M22 + 1;
    R1 := R1 + xx[k]*yy[k]*yy[k];
    R2 := R2 + xx[k]*yy[k];
  end;
  D := M11*M22-M12*M12;
  A := (M22*R1-M12*R2)/D;
  B := (M11*R2-M12*R1)/D;
end;

procedure kijkuit(jaar,licht : lijst; var A,B,R : double);
{
  Least Squares best fit
  to hyperbola y=B/(x-A)
  with linear regression
}
var
  N,k : integer;
  D : double;
  xx,yy : array of double;
  M11,M22,M12,R1,R2 : double;
begin
  N := Length(jaar);
  SetLength(xx,N);
  SetLength(yy,N);
  for k := 0 to N-1 do
  begin
    xx[k] := jaar[k];
    yy[k] := licht[k];
  end;
  M11 := 0; M12 := 0; M22 := 0;
  R1 := 0; R2 := 0;
  for k := 0 to N-1 do
  begin
    if yy[k] <= 0 then Continue;
    M11 := M11 + 1;
    M12 := M12 + 1/yy[k];
    M22 := M22 + 1/(yy[k]*yy[k]); 
    R1 := R1 + xx[k];
    R2 := R2 + xx[k]/yy[k];
  end;
  D := M11*M22-M12*M12;
  A := (M22*R1-M12*R2)/D;
  B := (M11*R2-M12*R1)/D;
{ Residual }
  R := 0;
  for k := 0 to N-1 do
  begin
    if yy[k] <= 0 then Continue;
    R := R + sqr(xx[k]-A-B/yy[k]);
  end;
  R := sqrt(R/M11);
end;

procedure NewtonHdB(jaar,licht : lijst; var A,B : double; var eps : punt);
{
  Least Squares best fit
  to hyperbola y=B/(x-A)
  non-linear regression
}
var
  x,y,D : double;
  M,N : array[0..1,0..1] of double;
  u,v : array[0..1] of double;
  L,k : integer;
begin
  L := Length(jaar);
  M[0,0] := 0; M[0,1] := 0;
  M[1,1] := 0;
  u[0] := 0; u[1] := 0;
  for k := 0 to L-1 do
  begin
    x := jaar[k]; y := licht[k];
    M[0,0] := M[0,0] -2*B*(-3*B+2*x*y-2*y*A)/sqr(sqr(x-A));
    M[0,1] := M[0,1] -2*(-2*B+y*x-y*A)/(sqr(x-A)*(x-A));
    M[1,1] := M[1,1] +2/sqr(x-A);
    u[0] := u[0] -2*(y-B/(x-A))*B/sqr(x-A);
    u[1] := u[1] -2*(y-B/(x-A))/(x-A);
  end;
  M[1,0] := M[0,1];
  D := M[0,0]*M[1,1] - M[0,1]*M[1,0];
  N[0,0] := + M[1,1]/D;
  N[1,1] := + M[0,0]/D;
  N[0,1] := - M[0,1]/D;
  N[1,0] := - M[1,0]/D;
  v[0] := N[0,0]*u[0] + N[0,1]*u[1];
  v[1] := N[1,0]*u[0] + N[1,1]*u[1];
  eps.x := abs(v[0]);
  eps.y := abs(v[1]);
  A := A - v[0];
  B := B - v[1];
end;

procedure Leibovici(jaar,licht : lijst; var A,B : double; var eps : double);
{
  Least Squares best fit
  to hyperbola y=B/(x-A)
  non-linear regression
}
var
  F,dF,x,y : double;
  P,Q,R,S : double;
  dP,dQ,dR,dS{,dB} : double;
  L,k : integer;
begin
  L := Length(jaar);
  P := 0; Q := 0;
  R := 0; S := 0;
  dP := 0; dQ := 0;
  dR := 0; dS := 0;
  for k := 0 to L-1 do
  begin
    x := jaar[k];
    y := licht[k];
    P := P + y/(x-A);
    dP := dP + y/sqr(x-A);
    Q := Q + 1/(sqr(x-A)*(x-A));
    dQ := dQ + 3/sqr(sqr(x-A));
    R := R + y/sqr(x-A);
    dR := dR + 2*y/(sqr(x-A)*(x-A));
    S := S + 1/sqr(x-A);
    dS := dS + 2/(sqr(x-A)*(x-A));
  end;
  F := P*Q - R*S;
  dF := (P*dQ+dP*Q) - (R*dS+dR*S);
  eps := abs(F/dF);
  A := A - F/dF;
  B := P/S;
{ dB := abs(B - R/Q);
  Writeln(B,dB); }
end;

procedure Raaklijn(jaar,licht,w : lijst; var A,B,dA : double);
{
  Least Squares best fit
  to tangent straight line
  of hyperbola y=B/(x-A)
  with linear regression
}
var
  tel,k : integer;
  m11,m12,m22,R1,R2,
  gewicht,D,H,C,R : double;
begin
  tel := Length(jaar);
  gewicht := 0;
  for k := 0 to tel-1 do
  begin
    gewicht := gewicht + w[k];
  end;
  m11 := 0;  m12 := 0;
  R1 := 0; R2 := 0;
  for k := 0 to tel-1 do
  begin
    m11 := m11 + w[k]*sqr(jaar[k]);
    m12 := m12 + w[k]*jaar[k];
    R1 := R1 + w[k]*jaar[k]*licht[k];
    R2 := R2 + w[k]*licht[k];
  end;
  m11 := m11/gewicht;
  m12 := m12/gewicht;
  m22 := 1;
  R1 := R1/gewicht;
  R2 := R2/gewicht;
  D := m11*m22 - sqr(m12);
  H := (m22*R1 - m12*R2)/D;
  C := (m11*R2 - m12*R1)/D;
  R := 0; 
  for k := 0 to tel-1 do
  begin
    R := R + w[k]*sqr(licht[k]-H*jaar[k]-C);
  end;
  R := sqrt(R/gewicht); { Residual }
  A := C/H; B := -C*A; dA := abs(A/C)*R;
end;

procedure Jacquelin(jaar,licht : lijst; var A,B : double);
{
  Least Squares best fit
  to hyperbola y=B/(x-A)
  with linear regression
}
var
  midden : punt;
  sigma_xx,sigma_yy,sigma_xy,
  Sp,Det,helling,kort,wortel : double;
{ lang,hoek,dA,delta : double; }
  i,teller : integer;
  verzamel : punten;
begin
  teller := Length(jaar);
  SetLength(verzamel,teller);
{ Straight line data transform }
  for i := 0 to teller-1 do 
  begin
    verzamel[i].x := licht[i];
    verzamel[i].y := jaar[i]*licht[i];
  end;
{ Mean / Center of gravity: }
  midden.x := 0;
  midden.y := 0;
  for i := 0 to teller-1 do
  begin
    midden.x := midden.x + verzamel[i].x;
    midden.y := midden.y + verzamel[i].y;
  end;
  midden.x := midden.x/teller;
  midden.y := midden.y/teller;
  for i := 0 to teller-1 do
  begin
    verzamel[i].x := verzamel[i].x - midden.x;
    verzamel[i].y := verzamel[i].y - midden.y;
  end;
{ Moments of inertia: }
  sigma_xx := 0;
  sigma_yy := 0;
  sigma_xy := 0;
  for i := 0 to teller-1 do
  begin
    sigma_xx := sigma_xx + sqr(verzamel[i].x);
    sigma_yy := sigma_yy + sqr(verzamel[i].y);
    sigma_xy := sigma_xy + verzamel[i].x*verzamel[i].y;
  end;
  sigma_xx := sigma_xx/teller;
  sigma_yy := sigma_yy/teller;
  sigma_xy := sigma_xy/teller;

{ hoek := arctan(2*sigma_xy/(sigma_xx-sigma_yy))/2;
  A := -1/tan(hoek);
  Writeln(A); KLOPT
}
  Sp := (sigma_xx + sigma_yy)/2;
  Det := sigma_xx*sigma_yy-sqr(sigma_xy);
  wortel := sqrt(sqr(Sp)-Det);
{ lang := Sp + wortel;} kort := Sp - wortel;
  helling := -sigma_xy/(sigma_yy-kort);
  A := -1/helling;
  B := midden.y - A*midden.x;
{ Writeln(A)
  dA := sqrt(kort/lang);
  hoek := arctan(A); delta := arctan(dA);
  Amin := tan(hoek-delta); ERROR ESTIMATE
  Amax := tan(hoek+delta); DOES NOT WORK!
 
  helling := -sigma_xy/(sigma_yy-lang);
  A := -1/helling;
  Writeln(A); KLOPT
  hoek := arctan(2*sigma_xy/(sigma_xx-sigma_yy))/2+Pi/2;
  Writeln(-1/tan(hoek));
}
end;

END.