Unit Ogenblik;
{
  PunchCard Reader
  ================
  This program has been designed and is CopyLefted by:

* Han de Bruijn; Systems for Research "A little bit of Physics would be  (===)
* and Education (DTO/SOO), Mekelweg 6  NO Idleness in Mathematics"(HdB) @-O^O-@
* 2628 CD Delft, The Netherlands  http://huizen.dto.tudelft.nl/deBruijn  #/_\#
* E-mail: Han.deBruijn@DTO.TUDelft.NL  Tel: +31 15 27 82751. Fax: 81722   ###

All I ask is to be credited when it is appropriate.

  Moments of a Contour
  --------------------
}
INTERFACE

Uses Algemeen;

type
  moment = record
    Inhoud : double; { Length / Area / Volume }
    Spoor, Determinant : double; { Trace & Determinant }
    Traagheid : tensor; { Moments of Inertia }
    Midden : vektor; { Centre of Gravity }
    Hoofdas : vektor; { Main Axis of Inertia }
    spreiding : vektor; { Main Standard Deviations }
  end;
  momenten = array of moment;

  Contour = class(TObject)
  public
    rond : Jordan; { Contour Vertices }
    vlak : moment; { Contour as a Surface }
    lijn : moment; { Contour as a Curve }
  public
    procedure Definitie(curve : Jordan); { Define Contour }
    function Binnen(r : vektor) : boolean; { Point Inside }
    procedure Eenvoudig; { Calculate Length / Area / Volume }
    procedure Hoofdzaken; { Calculate all Inertial Quantities }
    function AfbeeldenOp(andere : Contour) : Jordan; { Mapping }
  private
    procedure VlakRekenen; { Contour as a Surface }
    procedure LijnRekenen; { Contour as a Curve }
    procedure Belangrijk(var welk : moment); { Common }
  end;

  Ponskaart = class(Contour)
  private
    geheel : boolean; { Whole, not a Part }
    dikte : double; { Curve Thickness }
  public
    deltaOpp, deltaSP, deltaDet : double; { Error Estimates }
    Resolutie : double; { Resolution in pixels / meter }
    procedure Omtrek(eps : double); { Whole, not a Part }
    procedure EenGat(eps : double); { Part, not the Whole }
    procedure Grootheden; { Essential Quantities }
    procedure StartWaarden(andere: moment;
                           var A,p,q : vektor); { Starting Values }
    function Georienteerd(andere : Contour; blah : boolean) : double;
  private
    function LinksBoven(andere : Contour) : double;
  end;

procedure Alstublieft(oplossend : double); { Just a Test }

{ Routines for finding the Precise Positions of Punched Holes }
function Vergelijkingen(ponsing : momenten; A,p,q : vektor;
                        var stelsel : matrix; var rechts : kolom) : double;

IMPLEMENTATION

procedure Contour.Definitie;
{
  Definition
}
begin
  rond := curve;
end;

function stukje(a,b : vektor) : double;
begin
  stukje := sqrt(sqr(b.x-a.x) + sqr(b.y-a.y));
end;

procedure Contour.Eenvoudig;
{
  Volume inside a Contour
}
var
  k,L : integer;
  Det,Opp,S,ds : double;
  p,q : vektor;
begin
  Opp := 0;
  L := Length(rond);
  
  q := rond[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := rond[k];
    Det := (p.x*q.y - q.x*p.y);
    Opp := Opp + Det;
  end;
  vlak.Inhoud := Opp/2;
{
  Length of Contour
}
  L := Length(rond);
  S := 0;
  q := rond[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := rond[k];
    ds := stukje(p,q);
    S := S + ds;
  end;
  lijn.Inhoud := S;
end;

procedure Methode1(sigma : tensor; waarde : double;
                   var hoofd : vektor);
{ First Method
  for obtaining Eigenvectors
}
var
  lengte : double;
begin
  if sigma.yy < sigma.xx then
  begin
    hoofd.x := waarde - sigma.yy;
    hoofd.y := sigma.xy;
  end;
  if sigma.xx < sigma.yy then
  begin
    hoofd.x := sigma.xy;
    hoofd.y := waarde - sigma.xx;
  end;
  if sigma.xx = sigma.yy then
  begin
    hoofd.x := 1;
    hoofd.y := 1;
  end;
  if sigma.xy = 0 then
  begin
    hoofd.x := 1;
    hoofd.y := 0;
  end;
  lengte := sqrt(sqr(hoofd.x) + sqr(hoofd.y));
  hoofd.x := hoofd.x/lengte;
  hoofd.y := hoofd.y/lengte;
end;

procedure Methode2(sigma : tensor; waarde : double;
                   var hoofd : vektor);
{ Second Method
  for obtaining Eigenvectors
}
var
  breuk : double;
  theta,c,s,one,two : double;
  speciaal : boolean;

  function uitgewerkt(u,v : double) : double;
  var
    p,q : double;
  begin
    p := (sigma.xx - waarde)*u + sigma.xy*v;
    q := sigma.xy*u + (sigma.yy - waarde)*v;
    uitgewerkt := sqr(p) + sqr(q);
  end;

begin
  theta := 0;
  if sigma.xx = sigma.yy then
    theta := pi/4;
  if sigma.xy = 0 then
    theta := 0;
  speciaal := (sigma.xx = sigma.yy) or (sigma.xy = 0);
  if not speciaal then
  begin
    breuk := 2*sigma.xy / (sigma.xx - sigma.yy);
    theta := ArcTan(breuk)/2;
  end;
  c := cos(theta);
  s := sin(theta);
  one := uitgewerkt(c,s);
  two := uitgewerkt(-s,c);
  if one > two then
  begin
    hoofd.x := c;
    hoofd.y := s;
  end else
  begin
    hoofd.x := - s;
    hoofd.y := c;
  end;
end;

procedure Contour.Belangrijk;
{
  Important Things
}
var
  hoofd : vektor;
  sigma : tensor;
  wortel,waarde : double;
begin
{ Correllations / Moments of inertia: }
  sigma := welk.Traagheid;

  welk.Spoor := sigma.xx + sigma.yy;
  welk.Determinant := sqrt(sigma.xx*sigma.yy - sqr(sigma.xy));

{ Greatest eigenvalue: }
  wortel := sqrt(sqr(sigma.xx - sigma.yy)/4 + sqr(sigma.xy));
  waarde := (sigma.xx + sigma.yy)/2 + wortel;
  welk.spreiding.x := sqrt(waarde);
  welk.spreiding.y := welk.Determinant/welk.spreiding.x;

  Methode2(sigma,waarde,hoofd);
{ Eigenvectors: }
  welk.Hoofdas := hoofd;
end;

procedure Contour.VlakRekenen;
{
  Contour as a Surface
}
var
  k,L : integer;
  O,w,f : double;
  rand : Jordan;
  t,xx,yy,xy : double;
  h : tensor;
  p,q,c : vektor;
begin
  Eenvoudig;
  L := Length(rond);
  O := vlak.Inhoud;
{
  Center of Gravity
}
  f := 1/(O*6);
  c.x := 0;
  c.y := 0;
  q := rond[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := rond[k];
    w := f*(p.x*q.y - q.x*p.y);
    c.x := c.x + w*(p.x+q.x);
    c.y := c.y + w*(p.y+q.y);
  end;
  vlak.Midden := c;

  SetLength(rand,L);
  for k := 0 to L-1 do
  begin
    rand[k].x := rond[k].x - c.x;
    rand[k].y := rond[k].y - c.y;
  end;
{
  Moments of Inertia
}
  f := 1/(O*24);
  xx := 0;
  yy := 0;
  xy := 0;
  q := rand[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := rand[k];
    w := f*(p.x*q.y - q.x*p.y);
    t := 2*(p.x*p.y + q.x*q.y) + p.x*q.y + q.x*p.y;
    xy := xy + w*t;
    t := 2*(p.x*p.x + q.x*q.x + p.x*q.x);
    xx := xx + w*t;
    t := 2*(p.y*p.y + q.y*q.y + p.y*q.y);
    yy := yy + w*t;
  end;
  SetLength(rand,0);
  h.xx := xx;
  h.yy := yy;
  h.xy := xy;
  vlak.Traagheid := h;
  Belangrijk(vlak);
end;

procedure Contour.LijnRekenen;
{
  Contour as a Curve
}
var
  k,L : integer;
  ds,S : double;
  t,xx,yy,xy : double;
  rand : Jordan;
  p,q,c : vektor;
begin
  S := lijn.Inhoud;
  L := Length(rond);
{
  Central Moments
}
  q := rond[L-1];
  c.x := 0;
  c.y := 0;
  for k := 0 to L-1 do
  begin
    p := q;
    q := rond[k];
    ds := stukje(p,q)/(S*2);
    c.x := c.x + ds*(p.x + q.x);
    c.y := c.y + ds*(p.y + q.y);
  end;
  lijn.Midden := c;

  SetLength(rand,L);
  for k := 0 to L-1 do
  begin
    rand[k].x := rond[k].x - c.x;
    rand[k].y := rond[k].y - c.y;
  end;
  xx := 0;
  yy := 0;
  xy := 0;
  q := rand[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := rand[k];
    ds := stukje(p,q);
    t := (p.x*p.x + p.x*q.x + q.x*q.x)/3;
    xx := xx + t*ds/S;
    t := (p.y*p.y + p.y*q.y + q.y*q.y)/3;
    yy := yy + t*ds/S;
    t := (p.x*p.y + q.x*q.y)/3 + (p.x*q.y + p.y*q.x)/6;
    xy := xy + t*ds/S;
  end;
  SetLength(rand,0);
  lijn.Traagheid.xx := xx;
  lijn.Traagheid.yy := yy;
  lijn.Traagheid.xy := xy;
  Belangrijk(lijn);
end;

procedure Contour.Hoofdzaken;
begin
  VlakRekenen;
  LijnRekenen;
end;

procedure Ponskaart.Omtrek;
begin
  geheel := true;
  dikte := eps;
end;

procedure Ponskaart.EenGat;
begin
  geheel := false;
  dikte := eps;
end;

procedure Ponskaart.Grootheden;
{
  Punched Card Characteristics
}
var
  rand : Jordan;
  k,L : integer;
  faktor : double;
begin
  if geheel then
  begin
    L := 5;
    SetLength(rand,L);
    rand[4].x := 0;
    rand[4].y := 0;
    rand[3].x := 187;
    rand[3].y := 0;
    rand[2].x := 187;
    rand[2].y := -72;
    rand[1].x := 181;
    rand[1].y := -83;
    rand[0].x := 0;
    rand[0].y := -83;
  end else
  begin
    L := 4;
    SetLength(rand,L);
    rand[0].x := 0;
    rand[0].y := 0;
    rand[1].x := 1.425;
    rand[1].y := 0;
    rand[2].x := 1.425;
    rand[2].y := -3.5;
    rand[3].x := 0;
    rand[3].y := -3.5;
  end;

{ Assuming 100 DPI gives
  0.254 mm per pixel ==>
  Resolutie = 1/0.254 }

  for k := 0 to L-1 do
  begin
    rand[k].x := rand[k].x * Resolutie;
    rand[k].y := rand[k].y * Resolutie;
  end;
  Definitie(rand);
  SetLength(rand,0);
  Hoofdzaken;

  deltaOpp := lijn.Inhoud * dikte;
  faktor := dikte * lijn.Inhoud / abs(vlak.Inhoud);
  deltaSp := lijn.Spoor * faktor;
  deltaDet := lijn.Determinant * faktor;
end;

function Ponskaart.LinksBoven;
{
  Upper Left Corner of PunchCard
}
var
  u,v,r,p,q,m : vektor;
  k,L : integer;
  sigma,vgl,w,som : double;
begin
  u.x := 187 * Resolutie;
  u.y := -72 * Resolutie;
  v.x := 181 * Resolutie;
  v.y := -83 * Resolutie;
{ Mean }
  r.x := (u.x+v.x)/2; 
  r.y := (u.y+v.y)/2; 
  r := Verplaats(r,vlak.Midden,-1);
  r := Roteer(r,vlak.Hoofdas,-1);
  r := Roteer(r,andere.vlak.Hoofdas,+1);
  r := Verplaats(r,andere.vlak.Midden,+1);
{ Squared Spread }
  sigma := sqr(u.x-v.x)+sqr(u.y-v.y);
  sigma := sigma/sqr(2*3);
  L := Length(andere.rond);
  som := 0;
  q := andere.rond[L-1];
  for k := 0 to L-1 do
  begin
    p := q;
    q := andere.rond[k];
    m.x := (p.x+q.x)/2;
    m.y := (p.y+q.y)/2;
    vgl := sqr(r.x-m.x) + sqr(r.y-m.y);
    if vgl > sqr(2*pi)*sigma then Continue;
    w := sqrt(sqr(p.x-q.x)+sqr(p.y-q.y));
    som := som + exp(-vgl/(2*sigma)) * w;
  end;
  LinksBoven := som;
end;

function Ponskaart.Georienteerd;
{
  Orientation of Model Coordinates
}
var
  OK : boolean;
  een,twee,deze : double;

  procedure Omkeren(var v : vektor);
  {
    Reverse Direction
  }
  begin
    v.x := - v.x;
    v.y := - v.y;
  end;

begin
  OK := true;
  een := LinksBoven(andere);
  Omkeren(vlak.Hoofdas);
  twee := LinksBoven(andere);
  if twee > een then OK := false;
  if OK then
  begin
    Omkeren(vlak.Hoofdas);
    if blah then
      Writeln('Axes Allright');
    deze := een;
  end else
  begin
    if blah then
      Writeln('Axes Reversed');
    deze := twee;
  end;
  if blah then
    Writeln('Best Fit = ',deze:10:5);
  Georienteerd := deze;
end;

function Contour.Binnen;
{
  Decide whether a Point is Inside
}
var
  k,L : integer;
  p,q : vektor;
  rondom : double;
begin
  L := Length(rond);
  q := rond[L-1];
  q.x := q.x - r.x;
  q.y := q.y - r.y;
  rondom := 0;
  for k := 0 to L-1 do
  begin
    p := q;
    q := rond[k];
    q.x := q.x - r.x;
    q.y := q.y - r.y;
    rondom := rondom + Hoek(p,q);
  end;
  Binnen := (rondom/Pi > 1);
end;

function Contour.AfbeeldenOp;
{
  Map This contour upon the Other one
}
var
  k,L : integer;
  rand : Jordan;
  p : vektor;
begin
  L := Length(rond);
  SetLength(rand,L);
  for k := 0 to L-1 do
  begin
    p := rond[k];
    p := Verplaats(p,vlak.Midden,-1);
    p := Roteer(p,vlak.Hoofdas,-1);
    p := Roteer(p,andere.vlak.Hoofdas,+1);
    p := Verplaats(p,andere.vlak.Midden,+1);
    rand[k] := p;
  end;
  AfbeeldenOp := rand;
end;

function Vergelijkingen;
{
  System of Equations
}
var
  sigma,hulp : tensor;
  mu,r : vektor;
  weegt,term,gewicht : double;
  k,L,i,j : integer;
  verdeling : Schuine;
begin
  SetLength(stelsel,6+1,6+1);
  for i := 0 to 6 do
  for j := 0 to 6 do
    stelsel[i,j] := 0;
  SetLength(rechts,6+1);
  for i := 0 to 6 do
    rechts[i] := 0;

  verdeling := Schuine.Create;
  L := Length(ponsing);
  gewicht := 0;
  for k := 1 to L-1 do
  begin
    sigma := ponsing[k].Traagheid;
    mu := ponsing[k].Midden;
    verdeling.Definitie(sigma,mu);
    sigma := verdeling.Gewogen;

    for i := 0 to 80-1 do
    for j := 0 to 12-1 do
    begin
      r.x := a.x + i*p.x + j*q.x;
      r.y := a.y + i*p.y + j*q.y;
      weegt := verdeling.Gauss(r);
      gewicht := gewicht + weegt;

    { By partial differentiation of 
      this exponential function to
      a.x, a.y, p.x, p.y, q.x, q.y
      & putting the result to zero
      linear equations are obtained }

      hulp.xx := weegt*sigma.xx;
      hulp.yy := weegt*sigma.yy;
      hulp.xy := - weegt*sigma.xy;

      stelsel[1,1] := stelsel[1,1] + 1*hulp.yy;
      stelsel[1,2] := stelsel[1,2] + i*hulp.yy;
      stelsel[1,3] := stelsel[1,3] + j*hulp.yy;
      stelsel[1,4] := stelsel[1,4] + 1*hulp.xy;
      stelsel[1,5] := stelsel[1,5] + i*hulp.xy;
      stelsel[1,6] := stelsel[1,6] + j*hulp.xy;

      stelsel[2,1] := stelsel[2,1] + i*1*hulp.yy;
      stelsel[2,2] := stelsel[2,2] + i*i*hulp.yy;
      stelsel[2,3] := stelsel[2,3] + i*j*hulp.yy;
      stelsel[2,4] := stelsel[2,4] + i*1*hulp.xy;
      stelsel[2,5] := stelsel[2,5] + i*i*hulp.xy;
      stelsel[2,6] := stelsel[2,6] + i*j*hulp.xy;

      stelsel[3,1] := stelsel[3,1] + j*1*hulp.yy;
      stelsel[3,2] := stelsel[3,2] + j*i*hulp.yy;
      stelsel[3,3] := stelsel[3,3] + j*j*hulp.yy;
      stelsel[3,4] := stelsel[3,4] + j*1*hulp.xy;
      stelsel[3,5] := stelsel[3,5] + j*i*hulp.xy;
      stelsel[3,6] := stelsel[3,6] + j*j*hulp.xy;

      stelsel[4,1] := stelsel[4,1] + 1*hulp.xy;
      stelsel[4,2] := stelsel[4,2] + i*hulp.xy;
      stelsel[4,3] := stelsel[4,3] + j*hulp.xy;
      stelsel[4,4] := stelsel[4,4] + 1*hulp.xx;
      stelsel[4,5] := stelsel[4,5] + i*hulp.xx;
      stelsel[4,6] := stelsel[4,6] + j*hulp.xx;

      stelsel[5,1] := stelsel[5,1] + i*1*hulp.xy;
      stelsel[5,2] := stelsel[5,2] + i*i*hulp.xy;
      stelsel[5,3] := stelsel[5,3] + i*j*hulp.xy;
      stelsel[5,4] := stelsel[5,4] + i*1*hulp.xx;
      stelsel[5,5] := stelsel[5,5] + i*i*hulp.xx;
      stelsel[5,6] := stelsel[5,6] + i*j*hulp.xx;

      stelsel[6,1] := stelsel[6,1] + j*1*hulp.xy;
      stelsel[6,2] := stelsel[6,2] + j*i*hulp.xy;
      stelsel[6,3] := stelsel[6,3] + j*j*hulp.xy;
      stelsel[6,4] := stelsel[6,4] + j*1*hulp.xx;
      stelsel[6,5] := stelsel[6,5] + j*i*hulp.xx;
      stelsel[6,6] := stelsel[6,6] + j*j*hulp.xx;

      term := hulp.yy*mu.x + hulp.xy*mu.y;
      rechts[1] := rechts[1] +   term;
      rechts[2] := rechts[2] + i*term;
      rechts[3] := rechts[3] + j*term;
      term := hulp.xy*mu.x + hulp.xx*mu.y;
      rechts[4] := rechts[4] +   term;
      rechts[5] := rechts[5] + i*term;
      rechts[6] := rechts[6] + j*term;
    end;
  end;
{ Residual as a measure
  for overall Accuracy }
  Vergelijkingen := gewicht;
end;

procedure Alstublieft;
{
  Just a Test
}
var
  proef : Ponskaart;
  uit : moment;
begin
  proef := Ponskaart.Create;
  proef.Omtrek(oplossend);
  proef.Grootheden;
  uit := proef.vlak;
  Writeln('vlak.Inhoud = ',uit.Inhoud);
  Writeln('vlak.Spoor = ',uit.Spoor);
  Writeln('vlak.Determinant = ',uit.Determinant);
  Writeln('vlak.Midden = (',uit.Midden.x,',',uit.Midden.y,')');
  Writeln('vlak.Hoofdas = (',uit.Hoofdas.x,',',uit.Hoofdas.y,')');
  Writeln('vlak.spreiding = (',uit.spreiding.x,',',uit.spreiding.y,')');
  uit := proef.lijn;
  Writeln('lijn.Inhoud = ',uit.Inhoud);
  Writeln('lijn.Spoor = ',uit.Spoor);
  Writeln('lijn.Determinant = ',uit.Determinant);
  Writeln('lijn.Midden = (',uit.Midden.x,',',uit.Midden.y,')');
  Writeln('lijn.Hoofdas = (',uit.Hoofdas.x,',',uit.Hoofdas.y,')');
  Writeln('lijn.spreiding = (',uit.spreiding.x,',',uit.spreiding.y,')');
end;

procedure Ponskaart.StartWaarden;
{
  Initial Values
}
var
  h : vektor;
begin
  h.x := (5 + 1.425/2) * Resolutie;
  h.y := (5 + 3.500/2) * Resolutie;
  h.y := - h.y;
  h := Verplaats(h,vlak.Midden,-1);
  h := Roteer(h,vlak.Hoofdas,-1);
  h := Roteer(h,andere.Hoofdas,+1);
  h := Verplaats(h,andere.Midden,+1);
  A := h;

  h.x := 2.225 * Resolutie;
  h.y := 0;
  h := Roteer(h,vlak.Hoofdas,-1);
  h := Roteer(h,andere.Hoofdas,+1);
  p := h;

  h.x := 0;
  h.y := 6.315 * Resolutie;
  h.y := - h.y;
  h := Roteer(h,vlak.Hoofdas,-1);
  h := Roteer(h,andere.Hoofdas,+1);
  q := h;
end;

END.