{
This software has been designed and it is
CopyLefted by Han de Bruijn:

                          (===)
                         @-O^O-@
                          #/_\#
                           ###

UNIFIED ALTERNATIVE COSMOLOGY GALAXY PROJECT
============================================
}
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 Unit9;

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.13;
var
  R0 : double;
{Sort of grid / density }
  grid : integer;

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,C : double;
begin
  C := 2*Pi/(sqrt(1+1/sqr(sigma)) + 1+1/sqr(sigma));
  h := 1+sqr(R/sigma);
  D := 1/(C*h*sqrt(h));
{ https://www.astro.umd.edu/~richard/ASTRO421/
  A421_Elliptical%20Galaxies_Lec18.pdf page 19 }
  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;

function echt : double;
{
  Dimensionless G.W/(V^2.L)
}
const
{ Gravity constant (G) }
  G : double = 6.67408E-11;
{ Speed of light }
  light : double = 299792458;
{ Number of seconds in a year }
  year : double = 31556926;
{ Mass of our sun }
  sun : double = 1.989E30;
{ Speed of sun (V) }
  speed : double = 250000;
var
  radius, mass : double;
begin
{ Radius of milky way (L) }
  radius := 52850*light*year;
{ Mass of milky way (W) }
  mass := 1.5E12*sun;
{ Dimensionless number }
  echt := G*mass/(radius*sqr(speed));
end;

procedure profiel(kleur : TColor);
{
  Orbital speed profile
}
const
  schaal : double = 0.25;
  gok : double = 3.837911894;
var
  i,j,k,kk : integer;
  x,y,C : double;
begin
  Form1.Image1.Canvas.Pen.Width := 2;
  Form1.Image1.Canvas.Pen.Color := clRed;
  y := rho(0)*0.25*schaal;
  Form1.Image1.Canvas.MoveTo(x2i(0),y2j(y));
  for k := 1 to N do
  begin
    x := k/N ; y := rho(x)*0.25*schaal;
    Form1.Image1.Canvas.LineTo(x2i(x),y2j(y));
  end;
  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;
  Writeln('First ring ',k,': radius = ',straal(k));

  C := echt;
{ y is assumed positive }
  i := x2i(straal(kk));
  y := C*Ring(kk); { speed }
  j := y2j(schaal*sqrt(y));
  Form1.Image1.Canvas.MoveTo(i,j);
  for k := kk+1 to N do
  begin
    x := straal(k);
    y := C*Ring(k); { speed squared }
    i := x2i(x); j := y2j(schaal*sqrt(y));
    Form1.Image1.Canvas.LineTo(i,j);
  end;
{ Unified Alternative Cosmology }
  x := straal(kk);
  y := C*Ring(kk); { speed squared }
  i := x2i(x);
  j := y2j(schaal*sqrt(gok*x*y));
  Form1.Image1.Canvas.MoveTo(i,j);
  for k := kk+1 to N do
  begin
    x := straal(k);
    y := C*Ring(k); { speed squared }
    i := x2i(x);
    j := y2j(schaal*sqrt(gok*x*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 ; ');

  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;
  Form1.Image1.Picture.SaveToFile('2022jun.bmp');
end;

procedure TForm1.Scheppen(Sender: TObject);
{
  On Create
}
begin
  TV(Form1.Image1);
  ClearDevice;
  xmin := 0; xmax := 1;
  ymin := 0; ymax := 2;

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

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

procedure TForm1.Toetsdruk(Sender: TObject; var Key: Char);
{
  On Key Press
}
begin
  Exit;
end;

END.
=============================================================================
program Project9;

uses
  Forms,
  Unit9 in 'Unit9.pas' {Form1};

{$R *.res}

begin
  Application.Initialize;
  Application.CreateForm(TForm1, Form1);
  Application.Run;
end.