Unit Grafisch;

INTERFACE

Uses Extctrls, Types, Graphics;

var
  Wijd,Hoog : integer;
  xmin,xmax,ymin,ymax : double;

procedure TV(beeld : TImage);
procedure ClearDevice;
function Grijs(getal : double) : TColor;

function x2i(x : double) : integer;
function i2x(i : integer) : double;
function y2j(y : double) : integer;
function j2y(j : integer) : double;

IMPLEMENTATION

var
  scherm : TImage;

procedure TV(beeld : TImage);
begin
  scherm := beeld;
end;

procedure ClearDevice;
{
  Clear Screen
}
var
  rechthoek : TRect;
begin
  Wijd := scherm.Width-1;
  Hoog := scherm.Height-1;
  rechthoek := Rect(0,0,Wijd+1,Hoog+1);
  with scherm.Canvas do
  begin
    Brush.Color := clWhite;
    FillRect(rechthoek);
  end;
end;

function Grijs(getal : double) : TColor;
{
  Translates Real into Grey values
}
var
  even,k : byte;
  effe : Tcolor;
begin
  if getal < 0 then getal := 0;
  if getal > 1 then getal := 1;
  effe := 0;
  for k := 1 to 3 do
  begin
    even := Round(255*(1-getal));
    effe := (effe shl 8) or even;
  end;
  Grijs := $00000000 or effe;
end;

function x2i(x : double) : integer;
begin
  x2i := Round((x-xmin)/(xmax-xmin)*Wijd);
end;

function i2x(i : integer) : double;
begin
  i2x := xmin + i/Wijd*(xmax-xmin);
end;

function y2j(y : double) : integer;
begin
  y2j := Round(Hoog*(1-(y-ymin)/(ymax-ymin)));
end;

function j2y(j : integer) : double;
begin
  j2y := ymin + (Hoog-j)/Hoog*(ymax-ymin);
end;

END.

Unit Gemeen;

INTERFACE

var
  M,N : integer;

function macht(x : double; n : integer) : double;
procedure integraal(Q : double; var F1,F2 : double);

IMPLEMENTATION

function macht(x : double; n : integer) : double;
{
  Power = x^N
}
var
  q : double;
  k : integer;
begin
  q := 1;
  if n > 0 then
  for k := 1 to n do
    q := q * x;
  macht := q;
end;

procedure integraal(Q : double; var F1,F2 : double);
{
  RingWorld Integral
}
var
  k,kk : integer;
  dphi, term : double;

  function kracht(phi : double) : double;
  {
    Gravitational Force
  }
  var
    c,noemer : double;
  begin
    c := cos(phi);
    noemer := sqr(Q) - 2*Q*c + 1;
    kracht := (Q - c)/(noemer*sqrt(noemer));
  end;

begin
  dphi := 2*Pi/M;
  kk := M div 2;
  F1 := 0;
{ Midpoint rule }
  for k := 0 to kk-1 do
  begin
    term := kracht((k+0.5)*dphi);
    F1 := F1 + term*dphi;
  end;
  F1 := 2*F1;
{ Trapezoidal rule }
  F2 := 0;
  for k := 0 to kk-1 do
  begin
    term := kracht(k*dphi)+kracht((k+1)*dphi);
    F2 := F2 + term*dphi;
  end;
{ Writeln('F:',F1,' ',F2); }
end;

END.

unit Unit5;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, Grafisch, Gemeen;

type
  TForm1 = class(TForm)
    Image1: TImage;
    procedure Scheppen(Sender: TObject);
    procedure Toetsdruk(Sender: TObject; var Key: Char);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

const
{ Width of Galaxy Core }
  Sigma : double = 0.1;
var
  R0 : double;
{Sort of grid / density }
  grid,dicht : integer;
{ Normed density constants }
  cA,cB,cC,cD : double;

function straal(k : integer) : double;
{
  Radius of RingWorld
}
var
  D : double;
begin
  D := R0*exp(-ln(R0)*k/N);
  if grid = 1 then D := (k+1)*R0; 
  straal := D;
end;

function rho(R : double) : double;
{
  Mass density as a function
  of the RingWorld's radius
}
var
  D,h : double;
begin
  D := 1/cA;
  h := sqr(R)+sqr(Sigma);
  case dicht of
  1 : D := 1/(cB*h);
{ https://www.astro.umd.edu/~richard/ASTRO421/
  A421_Elliptical%20Galaxies_Lec18.pdf page 19 }
  2 : D := 1/(cC*h*sqrt(h));
  3 : D := exp(-sqr(R/Sigma)/2)/cD;
  end;
  rho := D;
end;

function Ring(i : integer) : double;
{
  Galaxy Integral
  for RingWorld (i)
}
var
  k : integer;
  Ri,Rk,F1,F2,extra : double;
  a1,a2,C,vorig,dR : double;
begin
  a1 := 0; a2 := 0;
  Ri := straal(i);
  Rk := 0;
  extra := 0;
  if grid = 0 then
  begin
    Rk := straal(-1);
  { Inner circle contribution }
    extra := rho(0)*Pi*sqr(Rk)/Ri;
  end;
  for k := 0 to N do
  begin
    if k = i then Continue;
  { Singular force skipped }
    vorig := Rk;
    Rk := straal(k);
    dR := Rk-vorig;
    integraal(Ri/Rk,F1,F2);
    C := Ri/Rk*rho(Rk)*dR;
    a1 := a1 + C*F1;
    a2 := a2 + C*F2;
  end;
{ Writeln('a:',a1,a2); }
  Ring := (a1+a2)/2 + extra;
end;

procedure profiel(kleur : TColor);
{
  Orbital speed profile
}
const
  schaal : double = 0.1;
var
  i,j,k,kk : integer;
  x,y : double;
begin
  Form1.Image1.Canvas.Pen.Color := clRed;
  y := rho(0)*schaal;
  Form1.Image1.Canvas.MoveTo(x2i(0),y2j(y));
  for k := 1 to N do
  begin
    x := k/N ; y := rho(x)*schaal;
    Form1.Image1.Canvas.LineTo(x2i(x),y2j(y));
  end;
  Form1.Image1.Canvas.Pen.Width := 2;
  Form1.Image1.Canvas.Pen.Color := kleur;
  kk := 0;
  for k := 0 to N do
  begin
    y := Ring(k); kk := k;
    if y > 0 then Break;
  end;
{ y is assumed positive }
  i := x2i(straal(kk));
  j := y2j(sqrt(Ring(kk)));
  Form1.Image1.Canvas.MoveTo(i,j);
  for k := kk+1 to N do
  begin
    x := straal(k);
    y := Ring(k);
  { Writeln(y); y := abs(y); }
    y := sqrt(y); { speed }
    i := x2i(x); j := y2j(y);
    Form1.Image1.Canvas.LineTo(i,j);
  end;
end;

procedure test(maal : integer);
{
  Tests
}
var
  keer : integer;
begin
  if grid = 0 then
    Write('Geometric grid ; ');
  if grid = 1 then
    Write('Arithmetic grid ; ');
  if dicht = 0 then
    Writeln('Uniform density');
  if dicht = 1 then
    Writeln('Density ~ Cauchy');
  if dicht = 2 then
    Writeln('Density ~ reference');
  if dicht = 3 then
    Writeln('Density ~ Gaussian');

  for keer := 1 to maal do
  begin
    R0 := macht((1-sin(pi/M))/(1+sin(pi/M)),N);
    if grid = 1 then R0 := 1/(N+1);
    profiel(Grijs(keer/maal));
    N := N*2; M := M*2;
  end;
end;

procedure TForm1.Scheppen(Sender: TObject);
{
  On Create
}
begin
  cA := Pi;
  cB := Pi*ln(1+1/sqr(Sigma));
  cC := 2*Pi*(1/Sigma-1/sqrt(1+sqr(Sigma)));
  cD := Sigma*sqrt(2*Pi);

  TV(Form1.Image1);
  ClearDevice;
  xmin := 0; xmax := 1;
  ymin := 0; ymax := 2;

  N := 20; M := 50;
  grid := 1;
  dicht := 0;
  test(4);

  N := 20; M := 50;
  grid := 0;
  test(4);
end;

procedure TForm1.Toetsdruk(Sender: TObject; var Key: Char);
{
  On Key Press
}
begin
  ClearDevice;
  N := 20; M := 50;
  grid := 1;
  dicht := dicht + 1;
  if dicht > 3 then dicht := 0;
  test(4);

  N := 20; M := 50;
  grid := 0;
  test(4);
end;

END.