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.