program Trouton_Noble_3; ........ Capaciter experiment

program darksim3; ...................... Curvature of space with difference equations

program darksim5; ..................... Curvature of space by ODE in spherical cases only

program haefele_und_keating; ....Clock experiment in the absolute space

program perihel_2; ......................Orbit of Mercury with additional acceleration term

program berechng; ......................Calculates straight forward different constantsprogram ro_d3; .......................... .Path lengthening within the Sun

program energy .......................... Enegy of a moved sphere filled with photons

ioncryc.cpp ................................. Potential of a central ion in a NaCl-type lattice

program Pioneer.cpp ................... Pioneer anomaly by additional acceleration

darksim9.pas ................................ Dark Matter ODE with C*w'^2 term

============================= *** =============================

program Trouton_Noble_3;

uses crt,dos;

(* Zentriertes Ellipsoidmodell vom 1.2.2000:

Berechnet die Energie des elektromagnetischen Feldes U=Integral E²±B² dV

zweier Ladungen im Abstand d=10 in Abhängigkeit von =v/c;

Centered ellipsoid model:

Calculates the energy of electromagnetical field U=integral E²-B² dV

of two charges with distance d=10 depending on b =v/c

Charge "p" positioned at (0',0',0'), charge "m" at (10',0',0')

Slashed coordinates are from the moved frame sigma'

*)

var i,j,k,shell:integer;

beta,xcms,ycms,eps,embq,wembq:real;

xs,ys,zs, xp,yp,zp, xm,ym,zm :real;

rp,rp2,rp3, rm,rm2,rm3 :real;

Exp,Eyp,Ezp, Exm,Eym,Ezm:real;

betasinap,Byp,Bzp, betasinam,Bym,Bzm:real;

n2Ep,nEp, n2Em,nEm:real;

Exps,Eyps,Ezps, Exms,Eyms,Ezms:real;

rmin:real;

Ue,Ueo,Ub,Ubo,E2,B2:real;

Uei,Ubi:array[0..500]of real;

begin

clrscr;

writeln('Trouton-Noble_3: Centered ellipsoid model Integral E²±B² dV');

writeln('===========================================================');

writeln;

writeln(' b Ue(b) Ub(b) ');

writeln(' __________________________________');

writeln;

beta:=0.0; (* =v/c *)

(* charge "p" in (0,0,0) *)

xcms:=0; (* x-coordinate of charge "m* *)

ycms:=5; (* y-coordinate of charge "m" *)

eps:=0.1; (* Epsilon for singularities *)

repeat (* cycle b =0, 0.05, ... *)

embq :=1-beta*beta; (* 1- b ² *)

wembq:=sqrt(embq); (* root of 1- b ² *)

(* Since the variability of the summands of the integral is too large,

it should be integrated in shells *)

for i:=0 to 500 do Uei[i]:=0; (* set shell sums to zero *)

for i:=0 to 500 do Ubi[i]:=0; (* set shell sums to zero *)

for i:=-50 to 60 do

begin

(* gotoxy(50,5);write('i=',i:4); *)

for j:=-50 to 60 do

for k:=-50 to 50 do

begin

xs:=i; (* sigma' coordinates *)

ys:=j;

zs:=k;

(* components of r from "p" (or "q") to (x,y,z) in sigma *)

xp:=xs/embq; yp:=ys/wembq; zp:=zs/wembq;

xm:=(xs-xcms)/embq; ym:=(ys-ycms)/wembq; zm:=zs/wembq;

(* r² in sigma from charge "p" or "m" to point (x,y,z) *)

rp2:= sqr(xp) + sqr(yp) + sqr(zp);

rm2:= sqr(xm) + sqr(ym) + sqr(zm);

rp := sqrt(rp2);

rm := sqrt(rm2);

if (rp>eps) and (rm>eps) then (* excluding singularities *)

begin

rp3:=rp2*rp;

rm3:=rm2*rm;

(* electrical field in point (x,y,z) in components *)

Exp:=xp/rp3; Eyp:=yp/rp3; Ezp:=zp/rp3;

Exm:=xm/rm3; Eym:=ym/rm3; Ezm:=zm/rm3;

(* Das magnetische Feld ist radialsymmetrisch zur x-Achse. Es hat

den Betrag b *sin(a)/r². a ist der Winkel zwischen r und x. Keine

x-Komponente. In x-z-Ebene senkrecht auf r.

Magnetic field is of radial symmetry to x-axis. Absolute value is

b *sin(a)/r². a is the angle between r and x. No x-component.

Lies in y-z-plane perpendicular to r.

*)

betasinap:=beta*sqrt((yp*yp+zp*zp)/rp2);

Byp:=-betasinap*Ezp;

Bzp:= betasinap*Eyp;

betasinam:=beta*sqrt((ym*ym+zm*zm)/rm2);

Bym:=-betasinam*Ezm;

Bzm:= betasinam*Eym;

(* squares and norms of the electrical and

magnetical fields in sigma *)

n2Ep:=sqr(Exp)+sqr(Eyp)+sqr(Ezp); nEp:=sqrt(n2Ep);

n2Em:=sqr(Exm)+sqr(Eym)+sqr(Ezm); nEm:=sqrt(n2Em);

(* Backward transformation of directions without change of norms *)

Exps:=nEp*xp/rp; Eyps:=nEp*yp/rp; Ezps:=nEp*zp/rp;

Exms:=nEm*xm/rm; Eyms:=nEm*ym/rm; Ezms:=nEm*zm/rm;

(* Radial symmetry of B is conserved under transformation *)

(* Calculate E² and B² *)

E2:= sqr(Exps-Exms)+sqr(Eyps-Eyms)+sqr(Ezps-Ezms);

B2:= sqr(Byp -Bym )+sqr(Bzp -Bzm );

rmin:=rp;

if rm<rp then rmin:=rm;

shell:=round(rmin);

Uei[shell]:=Uei[shell]+E2; (* Integration of E² *)

Ubi[shell]:=Ubi[shell]+B2; (* Integration of B² *)

end;

end;

end;

(* Sum the shells *)

Ue:=0;

for shell:=500 downto 0 do Ue:=Ue+Uei[shell];

Ub:=0;

for shell:=500 downto 0 do Ub:=Ub+Ubi[shell];

(* The integral is calculated *)

writeln(beta:9:2,Ue:15:6,Ub:15:6);

beta:=beta+0.1;

until beta>0.51;

writeln;

write('ENTER');

readln;

end.

============================== *** =============================

program darksim3;

(* Simulation of different types of galaxies in a lattice with difference

method

Simulation eines Raumgitters mit unterschiedlichen Belastungen

als Modelle von Galaxien im Weltraum

Especially is:

- matrix w is divided in 3 parts to enlarge rhe capacity. Turbo Pascal 6.0

allows only 64 kByte in one part

- Integration goes from centre to corner. Herefore is needed the index

matrix ijk consisting of 3 index vectors

Besonderheiten des Programms:

- Die w-Matrix wird aus Gruenden der eingeschraenkten Speicherverwaltung

von TurboPascal 6.0 in 3 getrennte Speicherbloecke unterteilt

- Die Integration erfolgt von "innen nach au b en". Dazu ist eine Indexmatrix

"ijk", bestehend aus drei Indexvektoren, erforderlich.

*)

uses crt,graph;

const KDR=0; (* Kontrolldruckparameter for debugging *)

GALAXIETYP=2; (* bar,Balken=2, disc,Scheibe=1, sphere,Kugel=0 *)

ZENTRALLAST=15000.0; (* Nur fuer Balken *)

LAST=1500.0; (* sphere,Kugel =11038.396

disc, Scheibe=11537.388

bar, Balken = 1500.000 *)

LASTFAKTOR=10.0; (* real enlargment of load*)

GWFAKTOR=1; (* enlarging only the numbers of load *)

SIGMA=0.07; (* sigma of density distribution of load *)

BKOEFF=0.0; (* w'ý-coeffizient 0.0002 *)

CKOEFF=0.0; (* w'-coeffizient 0.003 *)

DKOEFF=0.043; (* w-coeffizient 0.011 *)

JDISK=10; (* radius of disc in x-y-plane *)

vth_r10=0.1183;

(* sphere,Kugel =0.1177

disc, Scheibe=0.1180

bar, Balken =0.1183 bei sigma=0.07

=0.135 bei sigma=0.25

Theoretical w'-value at r=10 *)

MAXSTEP=500; (* maximal number of integration steps *)

GAMMA=0.02777729; (* Pseudo Gravitation constant of simulation *)

Imax=22; (* Number of rows in the lattice (y-direction) *)

Jmax=19; (* columns (x-direction) *)

Kmax=23; (* layers (z-direction) *)

NZIJK=31737; (* length of index vectors *)

NGW=5916; (* Number of lattice points until radius JDISK=10 *)

(* Die Koordinationen eines Punktes ijk cordinations of a point

1.Block: gerade Zeile, gerade Schicht 1. block even row even layer

2. unger. gerade 2. odd even

3. gerade unger. 3. even odd

4. unger. unger. 4. odd odd

coordination: how to find the 12 neighbours of a point

*)

kor:array[1..4,1..3,1..12]of integer=

(((0, 1, 1, 0,-1,-1, 0, 0,-1, 0, 0,-1),

(1, 0,-1, -1,-1, 0, 1, 0, 0, 1, 0, 0),

(0, 0, 0, 0, 0, 0, 1, 1, 1, -1,-1,-1)),

((0, 1, 1, 0,-1,-1, 0, 0,-1, 0, 0,-1),

(1, 1, 0, -1, 0, 1, 1, 0, 1, 1, 0, 1),

(0, 0, 0, 0, 0, 0, 1, 1, 1, -1,-1,-1)),

((0, 1, 1, 0,-1,-1, 1, 0, 0, 1, 0, 0),

(1, 0,-1, -1,-1, 0, -1,-1, 0, -1,-1, 0),

(0, 0, 0, 0, 0, 0, 1, 1, 1, -1,-1,-1)),

((0, 1, 1, 0,-1,-1, 1, 0, 0, 1, 0, 0),

(1, 1, 0, -1, 0, 1, 0,-1, 0, 0,-1, 0),

(0, 0, 0, 0, 0, 0, 1, 1, 1, -1,-1,-1)));

(* Gerades oder ungerades k (Schicht)

x-coordinate of neighbours for a point in a even or odd layer *)

xvek:array[1..2,1..12]of real=

((1.0, 0.5, -0.5, -1.0, -0.5, 0.5,

0.5, -0.5, 0.0, 0.5, -0.5, 0.0),

(1.0, 0.5, -0.5, -1.0, -0.5, 0.5,

0.0, -0.5, 0.5, 0.0, -0.5, 0.5));

(* y-coordinates of neighboures for a point in a even or odd layer *)

yvek:array[1..2,1..12]of real=

(( 0.0, 0.866, 0.866, 0.0, -0.866, -0.866,

0.289, 0.289, -0.577, 0.289, 0.289, -0.577),

( 0.0, 0.866, 0.866, 0.0, -0.866, 0.866,

0.577, -0.289, -0.289, 0.577, -0.289, -0.289));

(* z-coordinates dont depend on layer position *)

zvek:array[1..12]of real=

(0.0, 0.0, 0.0, 0.0, 0.0, 0.0,

0.816, 0.816, 0.816, -0.816, -0.816, -0.816);

type s9=string[10];

(* W-Speicherblock partitioned matrix of the w-coordinates of the lattice *)

wmat1=array[ -22..-8, -JMAX..Jmax, -KMAX..Kmax] of integer;

wmat2=array[ -7.. 7, -JMAX..Jmax, -KMAX..Kmax] of integer;

wmat3=array[ 8..22, -JMAX..Jmax, -KMAX..Kmax] of integer;

(* Index-vektor *)

indexvektor=array[1..NZIJK]of integer;

var sz,zeile,zeilnum:integer;

i,ii,j,k,l,m,n,kbl,kxy:integer;

ww:integer;

wmin,wijk:integer;

rmax,r,minr,maxr,ra,rb,rjd,rgrenz:real;

w3h,hdw3,w2d3:real;

gwsumme:real;

x,y,z:real;

f:real;

wert:real;

iterat:integer;

wkoord:integer;

swx,swy,swz,nablaw:real;

aa,bb,a,b,sx,sy,sxx,syy,syx,saqxx,saqyy,sapyx,xq,yq:real;

w1:^wmat1;

w2:^wmat2;

w3:^wmat3;

iiv,jiv,kiv:^indexvektor;

start,stopp:array[-IMAX..Imax,-KMAX..Kmax]of integer;

gw:array[1..NGW]of integer;

wr,vkurv,vkurvth:array[0..100]of real;

nr,nvkurv:array[0..100]of integer;

vkmax,dmk:real;

jd:integer;

vth,wsth,norm,skapro,r2:real;

ojdv,eojdv,oijv,rv,erv:array[1..3]of real;

wvek:array[1..12]of integer;

wrmin:real;

ket1:string[1];

ket6:string[6];

graphdriver,graphmode :integer;

antwort:char;

datei :text;

Function f9(x:real):s9;

(* converts number for output *)

var ket11:string[11];

begin

str(x:11,ket11);

if x<0 then f9:='-'+copy(ket11,2,6)+copy(ket11,10,2)

else f9:='+'+copy(ket11,2,6)+copy(ket11,10,2);

end;

Function lesew(i,j,k:integer):integer; (* read from partioned matrix *)

begin

case i of

-22..-8: lesew:=w1^[i,j,k];

-7.. 7: lesew:=w2^[i,j,k];

8..22: lesew:=w3^[i,j,k];

end;

end;

Procedure schreibew(i,j,k,w:integer); (* write to partioned matrix *)

begin

case i of

-22..-8: w1^[i,j,k]:=w;

-7.. 7: w2^[i,j,k]:=w;

8..22: w3^[i,j,k]:=w;

end;

end;

begin

clrscr;

writeln('Programm DARKSIM3');

getmem(w1,54990);

getmem(w2,54990);

getmem(w3,54990);

getmem(iiv,NZIJK*2);

getmem(jiv,NZIJK*2);

getmem(kiv,NZIJK*2);

writeln('Speicher ist allokiert storage allocated');

(* Berechne die (i,k)-Matrix der Reihenstopps

Jmax gibt die Entfernung der Randpunkte links rechts

A sphere inside the lattice is taken. Matrices START and STOPP

give the first and the last index of i,k to form a sphere

*)

rmax := Jmax-0.6; (* Maximal radius *)

w3h := sqrt(3.0)/2; (* distances in the cell *)

hdw3 := 0.5/sqrt(3.0);

w2d3 :=sqrt(2.0/3.0);

maxr :=-1e25;

minr := 1e25;

for i:=-Imax to Imax do

for k:=-Kmax to Kmax do

begin

start[i,k]:=99; (* starting values *)

stopp[i,k]:=99;

(* Berechne die Abstaende r vom Punkt 0,0,0 calculate distance r *)

for j:=-Jmax to Jmax do

begin

(* r depends on layer number (even or odd) *)

if (i mod 2)=0 then x:=j else x:=j+0.5;

if (k mod 2)=0 then y:=i*w3h else begin x:=x-0.5; y:=i*w3h+hdw3; end;

z:=k*w2d3;

r:=sqrt(sqr(x)+sqr(y)+sqr(z));

if (j<1) and (r<rmax) and (start[i,k]=99) then

begin

start[i,k]:=j;

if maxr<r then maxr:=r;

if minr>r then minr:=r;

end;

end;

for j:=Jmax downto -Jmax do

begin

if (i mod 2)=0 then x:=j else x:=j+0.5;

if (k mod 2)=0 then y:=i*w3h else begin x:=x-0.5; y:=i*w3h+hdw3; end;

z:=k*w2d3;

r:=sqrt(sqr(x)+sqr(y)+sqr(z));

if (j>-1) and (r<rmax) and (stopp[i,k]=99) then

begin

stopp[i,k]:=j;

if maxr<r then maxr:=r;

if minr>r then minr:=r;

end;

end;

end;

writeln('Maxr=',maxr:10:3);

writeln('Minr=',minr:10:3);

(* ---------------Kontrolle for control ------------------------------

for i:=0 to Imax do

begin

for k:=0 to kmax do write(stopp[i,k]:2);

writeln;

if (i mod 10)=0 then readln;

end;

readln;

---------------------------------------------------------- *)

(* Berechne die Matrix "ijk", d.h. die drei Indexvektoren iiv,jiv,kiv

mit den Gitterkoordinaten der hoehenveraenderlichen Punkte nach auf-

steigendem r vom Punkt (0.0314159265; 0.027182818; 0.014131711) aus

gerechnet. Dazu wird r in Schrittweite 0.25 von 0 bis maxr variiert.

Bei r=34.5 ist NZIJK=32098 erreicht.

Calculate matrix ijk, i.e. three vectors iiv,jiv,kiv. They contain the

point-indexes of all points sorted by ascending r. Distance taken not from

0,0,0, but from (0.0314159265; 0.027182818; 0.014131711) to bring some

randomness in the serie. r is varied from 0 to rmax in steps of 0.25.

At r=34.5 the maximum number of NZIJK=32098 is reached.

*)

ra:=0.0; (* Untere Intervallgrenze lower limit *)

zeilnum:=0; (* Indiziert die Zeile der Matrix ijk - the rows of "ijk" *)

repeat

rb:=ra+0.25;

for i:=-IMAX+1 to Imax-1 do

for k:=-KMAX+1 to Kmax-1 do

if (stopp[i,k]<>99)and(start[i,k]<>99) then

begin

for j:=start[i,k] to stopp[i,k] do

begin

(* radius r depends on layer number (even-odd) *)

if (i mod 2)=0 then x:=j else x:=j+0.5;

if (k mod 2)=0 then y:=i*w3h else begin x:=x-0.5; y:=i*w3h+hdw3; end;

z:=k*w2d3;

r:=sqrt(sqr(x-0.0314159265)+sqr(y-0.027182818)+sqr(z-0.014131711));

if (r>ra) and (r<=rb) then

begin

zeilnum:=zeilnum+1;

iiv^[zeilnum]:=i; (* Fill matrix *)

jiv^[zeilnum]:=j;

kiv^[zeilnum]:=k;

end;

end;

end;

writeln('Berechnung Matrix "ijk": r=',rb:6:3,' zeilnum=',zeilnum);

ra:=rb;

until zeilnum=NZIJK;

(* Kontrolle ---------------------for control-------------------- *)

writeln;writeln('Matrix ijk mit ',zeilnum,' Zeilen');

for i:=1 to 20 do writeln(iiv^[i]:5,jiv^[i]:5,kiv^[i]:5);

writeln('ENTER');

readln;

(* ------------------------------------------------------- *)

(* Berechne eine kugelfoermige Gewichtsmatrix GW mit NGW Punkten.

LAST verteilt sich zum Rand abnehmend auf alle Punkte.

ZENTRALLAST>0 wird zusaetzlich auf Punkt 0,0,0 gepackt.

Calculate mtrix of weights with spherical shape and NGW points.

Load (or weight) is so distributed, that it fades to the edge.

Central Load added to point 0,0,0

*)

gwsumme:=0; (* Sum of weights *)

for zeile:=1 to NGW do

begin

i:=iiv^[zeile];

j:=jiv^[zeile];

k:=kiv^[zeile];

(* Wenn Kugel, dann alle Punkte all points in case of sphere

Wenn Scheibe, dann nur in x-y-Ebene points in x-y-plane if disc

Wenn Balken, dann nur auf x-Achse points on x-axis if bar *)

gw[zeile]:=0;

if ((GALAXIETYP=2)and(k=0)and(i=0)) or

((GALAXIETYP=1)and(k=0)) or

(GALAXIETYP=0) then

begin

(* Abstand r des Punktes ijk vom Zentrum berechnen

radius r from centre to the point *)

if (i mod 2)=0 then x:=j else x:=j+0.5;

if (k mod 2)=0 then y:=i*w3h else begin x:=x-0.5; y:=i*w3h+hdw3; end;

z:=k*w2d3;

r:=sqrt(sqr(x)+sqr(y)+sqr(z));

(* =======Lastverteilung Load distribution ==================== *)

if (r<=JDISK) then

case GALAXIETYP of

0: gw[zeile]:=round(LAST*GWFAKTOR*exp(-sqr(r/SIGMA/JDISK)));

1: gw[zeile]:=round(LAST*GWFAKTOR*exp(-sqr(r/SIGMA/JDISK))*(1+r));

(* Der Balken besteht aus einer Zentrallast und einer linear

abfallenden Last, beide gleich gross

The bar has a central load and linear descendent weights.

Both parts of equal sum

*)

2: if zeile>1 then gw[zeile]:=round(LAST*(1-r/JDISK))

else gw[zeile]:=round(LAST*(1-r/JDISK)+ZENTRALLAST);

end;

gwsumme:=gwsumme+gw[zeile];

(* ====Ende Lastverteilung end load distribution ============ *)

end;

end;

writeln('GW-Summe =',gwsumme/GWFAKTOR*LASTFAKTOR/60000);

(* ------------------------------------Kontrolle------control------

writeln('Gewichte gw');

for i:=1 to 20 do writeln('gw[i]=',i:2,' ',gw[i]:6);

write('ENTER');readln;

-------------------------------------------------------------------- *)

(* Setze alle w-Auslenkungen 0 (kodiert als 30.000

w-matrix is of type integer to save storage. Starting value is

30.000 (means zero)

*)

for i:=-IMAX to Imax do

for j:=-JMAX to Jmax do

for k:=-KMAX to Kmax do schreibew(i,j,k,30000);

repeat

writeln;

write('Continue with old case? (j/n) : ');readln(antwort);

until antwort in ['j','n'];

if antwort='j' then

begin

assign(datei,'darksim3.txt');reset(datei);

for i:=-IMAX to Imax do

for j:=-JMAX to Jmax do

for k:=-KMAX to Kmax do

begin

readln(datei,wijk); (* get data from external file *)

schreibew(i,j,k,wijk);

end;

close(datei);

end;

sz:=0; (* Schrittzaehler step counter *)

rgrenz:=1.5*JDISK; (* r-limit *)

writeln('RGrenz=',rgrenz:10:5);

repeat

wmin:=32000; (* search minimum w *)

for i:=0 to 100 do begin vkurv[i]:=0; nvkurv[i]:=0; end; (* starting values *)

(* Berechne die Kraftresultierenden bei aufsteigender Indizierung

calculate the vectorial sums of forces in each point

during ascending indexing *)

for zeile:=1 to zeilnum do (* ascending indexing *)

begin

i:=iiv^[zeile];

j:=jiv^[zeile];

k:=kiv^[zeile];

(* Welcher der Koordinationsbloecke 1-4 ist zustaendig?

coordination depends on layer number (even-odd) *)

if (k mod 2) =0 then kbl:=1 else kbl:=3;

if (i mod 2)<>0 then kbl:=kbl+1;

(* Gerade Schichtnummer (kxy=1) oder ungerade (kxy=2)

even layer or odd *)

if (k mod 2) =0 then kxy:=1 else kxy:=2;

f:=0; (* Kraftresultierende vectorial sum of forces *)

wijk:=lesew(i,j,k); (* Auslenkung des Punktes ijk displacement of point *)

(* Die 12 Koordinationen zu Punkt ijk the 12 neighbours of point ijk *)

for m:=1 to 12 do

begin

wkoord:=lesew(i+kor[kbl,1,m], j+kor[kbl,2,m], k+kor[kbl,3,m]);

f:=f+wkoord-wijk;

wvek[m]:=wkoord-wijk;

end;

(* Berechnung des Gradienten Nabla-W gradient of w(x,y,z)

Zuerst die 3 Produktsummen first 3 sums of products *)

swx:=0; swy:=0; swz:=0;

for m:=1 to 12 do

begin

swx:=swx+wvek[m]*xvek[kxy,m];

swy:=swy+wvek[m]*yvek[kxy,m];

swz:=swz+wvek[m]*zvek[m];

end;

nablaw:=sqrt(sqr(swx)+sqr(swy)+sqr(swz))/4.0;

(* Belastung des Punktes mit nablaw*BKoeff

load point with Gradient(w)*Bcoefficient *)

f:=f-sqr(nablaw)*BKOEFF-nablaw*CKOEFF+DKOEFF*(wijk-30000);

(* In der Diskusebene bzw der Kugel: within disc or sphere *)

(* Ist dieser Punkt ijk mit einerm Gewicht belastet?

is this point with load? *)

if ((GALAXIETYP=2)and(k=0)and(i=0)) or

((GALAXIETYP=1)and(k=0)) or

(GALAXIETYP=0) then

begin

if zeile<=NGW then f:=f-gw[zeile]*LASTFAKTOR/GWFAKTOR;

end;

(* Merke den Gradienten und das r fuer die Kurvendarstellung

Hier reichen die Punkte auf der k=0-Ebene

Store the gradient and radius r for drawing curve

The points of k=0 plane are sufficiently *)

if k=0 then

if ((GALAXIETYP=2)and(i=0))or(GALAXIETYP<2) then

begin

if (i mod 2)=0 then x:=j else x:=j+0.5;

if (k mod 2)=0 then y:=i*w3h else begin x:=x-0.5; y:=i*w3h+hdw3; end;

z:=k*w2d3;

r:=sqrt(sqr(x)+sqr(y)+sqr(z));

if (r/rgrenz*100)<=100.49 then

begin

m:=round(r/rgrenz*100);

vkurv[m]:=vkurv[m]+sqrt(nablaw*r/60000.0);

nvkurv[m]:=nvkurv[m]+1;

end;

end;

(* Festlegung der Verschiebung in W-Richtung

calculate displacement in w-direction *)

wijk:=round(wijk+f/12);

schreibew(i,j,k,wijk);

if wmin>wijk then wmin:=wijk;

end;

dmk:=sqr(vkurv[67]/nvkurv[67]/vth_r10); (* dark matter coefficient *)

sz:=sz+1;

writeln('sz=',sz:4,' von ',MAXSTEP,' wmin=',wmin:7, ' dmk=',dmk:6:3,

' Taste q');

(* Abbruch auf Tastendruck q oder bei Erreichen der Iterationszahl

Stopp with key "q" or if reaching maximum step number *)

ii:=0;

if keypressed then begin ii:=ord(readkey); end;

until (sz=MAXSTEP) or (ii=113);

(* Ausgabe der w-Daten fuer eventuelle Nachiterationen

store w-data on file for later additional iterations *)

assign(datei,'darksim3.txt');rewrite(datei);

for i:=-IMAX to Imax do

for j:=-JMAX to Jmax do

for k:=-KMAX to Kmax do

writeln(datei,lesew(i,j,k):6);

close(datei);

(* Berechne die w(r)-Kurve mit 100 Stuetzstellen nur in der xy-Ebene

Calculate the w(r)-curve within 100 points only in x-y plane *)

for i:=0 to 100 do wr[i]:=0;

for i:=0 to 100 do nr[i]:=0;

wrmin:=0;

k:=0;

for i:=-IMAX+1 to Imax-1 do

if (start[i,k]<>99)and(stopp[i,k]<>99) then

begin

(* Berechne die Abstaende r vom Punkt 0,0,0

distance r from 0,0,0 *)

for j:=start[i,k] to stopp[i,k] do

begin

if (i mod 2)=0 then x:=j else x:=j+0.5;

y:=i*w3h;

r:=sqrt(sqr(x)+sqr(y));

if r<=rgrenz then

begin

wert:=lesew(i,j,k)/60000.0-0.5;

m:=round(r/rgrenz*100.0);

wr[m]:=wr[m]+wert; (* sum w-values *)

nr[m]:=nr[m]+1; (* count *)

end;

end;

end;

writeln('Summation of w(r) ready');

for i:=0 to 100 do if nr[i]>0 then

begin

wr[i]:=wr[i]/nr[i]; (* mean of this radius *)

if wrmin>wr[i] then wrmin:=wr[i];

end;

vkmax:=-1e25; (* Maximum for scaling plot *)

for i:=0 to 100 do if nvkurv[i]>0 then

begin

vkurv[i]:=vkurv[i]/nvkurv[i];

if vkmax<vkurv[i] then vkmax:=vkurv[i];

end;

writeln('vkmax=',vkmax:20:10);

write('ENTER');readln;

(* Berechnung der theoretischen w'(r)-Werte. Das Problem ist die

Vermeidung eines belasteten Punktes selbst.

Calculation of theoretical w'(r)-curve. Avoid to touch loaded points

since there the value is infinitely

*)

vkurvth[0]:=0.0;

for jd:=1 to 100 do

begin

wsth:=0;

(* Berechne Vektor ojdv vom Zentrum zum Kurvenpunkt auf der x-Achse

vector ojdv from centre in direction point on x-axis *)

rjd:=jd*rgrenz/100;

ojdv[1]:=rjd;

ojdv[2]:=0.43;

ojdv[3]:=0.43; (* 0.43 if bar, 0.0 if sphere or disc *)

norm:=sqrt(sqr(ojdv[1])+sqr(ojdv[2])+sqr(ojdv[3]));

for l:=1 to 3 do eojdv[l]:=ojdv[l]/norm;

(* Berechne die Summe der Einzelgravitationen zum Punkt (r,0.43,0)

Sum of gravitational influence of each point to point (r,0,0.43) *)

for zeile:=1 to NGW do

begin

i:=iiv^[zeile];

j:=jiv^[zeile];

k:=kiv^[zeile];

(* Berechne Vektor oijv vector from origin to (x,y,z) *)

if (i mod 2)=0 then x:=j else x:=j+0.5;

if (k mod 2)=0 then y:=i*w3h else begin x:=x-0.5; y:=i*w3h+hdw3; end;

z:=k*w2d3;

oijv[1]:=x;

oijv[2]:=y;

oijv[3]:=z;

(* Berechne r-Vektor rv=oijv-ojdv und rý

differnce vector *)

for l:=1 to 3 do rv[l]:=oijv[l]-ojdv[l];

r2:=sqr(rv[1])+sqr(rv[2])+sqr(rv[3]);

(* Bilde Einheitsvektor erv unit vector *)

r:=sqrt(r2);

for l:=1 to 3 do erv[l]:=rv[l]/r;

(* Berechne das Skalarprodukt eojdv*erv inner product *)

skapro:=0;

for l:=1 to 3 do skapro:=skapro+eojdv[l]*erv[l];

(* Addiere die Gravitationssumme auf sum gravitations *)

wsth:=wsth+skapro*gw[zeile]*LASTFAKTOR/GWFAKTOR/r2;

end;

(* scaling *)

wsth:=wsth*GAMMA/60000;

(* Jetzt Berechnung der Newton-v(r)-Kurve

Theoretical Newtonian v(r)-curve *)

vth:=sqrt(abs(-wsth*rjd));

vkurvth[jd]:=vth;

writeln('r,v(r)-th.',rjd:10:3,vth:10:5);

end;

write('ENTER');readln;

graphdriver:=0;graphmode:=0; (* Start graphics *)

initgraph(graphdriver,graphmode,'c:\tp\lib');

(* L"sche Schirm clear screen *)

setcolor(0);

for i:=1 to 479 do

begin

moveto(0,i);

lineto(639,i);

end;

(* Achsen axes *)

setcolor(2);

outtextxy(200,10,'Auslenkung w(r) der Membran');

moveto(100,40);lineto(100,460);

moveto(90,50);lineto(620,50);

for i:=1 to 10 do begin moveto(95,40*i+50);lineto(105,40*i+50);end;

for i:=1 to 10 do begin moveto(40*i+100,45);lineto(40*i+100,455);end;

outtextxy(505,55,'Rmax');

outtextxy(5,130,'W');

for i:=1 to 9 do

begin

str(i,ket1);

outtextxy(i*40+90,25,ket1+'0%');

end;

for i:=0 to 10 do

begin

x:=wrmin*i/10;

str(x:6:3,ket6);

outtextxy(40,i*40+45,ket6);

end;

(* Zeichne beobachtete Kurve draw observed curve *)

setcolor(2); outtextxy(500,400,'beob.w(r)');

j:=0;

for i:=0 to 100 do

if nr[i]>0 then

begin

if j=0 then begin moveto(100+4*i,50+round(wr[i]/wrmin*400.0));j:=1;end;

if j>0 then lineto(100+4*i,50+round(wr[i]/wrmin*400.0));

end;

(* Zeichne v-Kurve draw v-curve *)

setcolor(14); outtextxy(500,420,'beob.v(r)');

j:=0;

for i:=0 to 100 do

if nvkurv[i]>0 then

begin

if j=0 then begin moveto(100+i*4,450-round(vkurv[i]/vkmax*400.0));j:=1;end;

if j>0 then lineto(100+i*4,450-round(vkurv[i]/vkmax*400.0));

end;

(* Zeichne theoretische v-Kurve draw theoretical v(r)-curve *)

setcolor(12); outtextxy(500,440,'theo.v(r)');

j:=0;

for i:=5 to 100 do

begin

if j=0 then

begin

moveto(100+i*4,450-round(vkurvth[i]/vkmax*400.0));

j:=1;

end;

if j>0 then lineto(100+i*4,450-round(vkurvth[i]/vkmax*400.0));

end;

repeat until keypressed;

i:=ord(readkey);

closegraph;

(* Ausgabe des Radius und der 4 Kurven auf Datei und vmax=300 normiert

writing on file 4 curves normed to vmax=300 *)

assign(datei,'darksim3.dat');rewrite(datei);

for i:=2 to 100 do

begin

r:=i*rgrenz/100;

if vkurv[i]<>0.0 then

writeln(datei,f9(r),' ',f9(vkurv[i]*300/vkmax),' '

,f9(vkurvth[i]*300/vkmax),' '

,f9(wr[i]));

end;

close(datei);

for i:=60 to 70 do

begin

r:=i*rgrenz/100;

if vkurv[i]<>0.0 then

writeln('r,v,vth: ',i,' ',r,' ',vkurv[i],' '

,vkurvth[i]);

end;

end.

========================= *** ===============================

program darksim5;

uses crt,graph;

(* Berechnet die Biegung der Membran fuer folgendes Modell:

Computes the curvature of the membrane with the following model:

w"=-2w'/r + A*ro(r) - C*(w')ý - D*w*abs(w');

Factor "-C" means a force in w-direction (opposite of gravity)

Factor "-D" means (by negative sign of w) a force in direction of gravity

Der Unterschied zum Programm DARKSIM4.PAS sind die originalen Konstanten

in Kg, m, s und die dadurch notwendige logarithmische Rechnung.

The difference to program darksim4.pas is the use of the original constants

with Kg,m,s and the therefore needed logarithmic calculation

*)

const

lnC=-40.8; (* -40.8 w'-Koeffizient coefficient of w' *)

lnD=-59.4; (* -59.4 w-Koeffizient coefficient of w *)

W1o=-0.95e6; (* -0.95e+6 Startwert w1(r) initial value *)

W2o=-0.593e6; (* -0.593e6 Startwert w2(r) initial value *)

FRI=1.5; (* 1.5 Integrationsgrenzenfaktor *)

(* factor for integration boundary *)

Wo=1.432e+6; (* Raumtiefe am Sonnenrand space depths sun edge *)

Rs=6.9635e+8; (* Sonnenradius sun radius *)

Gs=280.1; (* Beschleunigung am Sonnenrand acceleration *)

gamma=6.6726e-11;

ha=1473.1; (* Halber Schwarzschildradius 1/2 Schw-Schildr *)

Ms=1.991e+30; (* Sonnenmasse sun mass *)

LJ=0.946728e+16; (* Lichtjahr Meter light year *)

FAKWs=1e15; (* Faktor fuer ws factor for w' *)

FAKWss=1e30; (* Faktor fuer wss factor for w" *)

MAXSTEP=25000.0; (* Zahl Integrationsschritte number of steps *)

SIGMA=0.2; (* 0.2 SIGMA Dichteverteilung sigma density distr. *)

RGLJ=50000.0; (* Galaxieradius Lichtjahre radius of galaxy *)

MG=5e10; (* Galaxiemasse Sonnenmassen mass of galaxy fac *)

type s9=string[10];

var

lnWo:real;

lnRs:real;

lnGs:real;

lnMs:real;

lnha:real;

lnRo,Ro:real; (* Dichtewert im Zentrum density at centre *)

a,lna:real; (* A-Koeffizient der DGL *)

Fo,lnFo:real; (* Membranspannung membrane tension *)

Ae,lnAe:real; (* Aetherbeschleunigung ether acceleration *)

r,lnr,r_in_Rg:real;(* Radius *)

RG,lnRG:real;

RI,lnRI:real; (* Integrationsradius radius of integration *)

dr,lndr:real; (* Schrittweite step width *)

lnFAKws:real;

lnFAKWss:real;

lnAterm:real;

lnBterm:real;

lnCterm:real;

lnDterm:real;

i,j:integer; (* Hilfsgroessen auxillary integers *)

w1,w1s,w1ss:real; (* w1(r), w1'(r), w1"(r) mit dunkler Materie *)

w1km,w2km:real; (* Raumtiefe in Km space depth in Km *)

vzw1,lnW1:real;

lnw1s,lnw1ss:real; (* resolving curves with dark matter *)

vzw1s:real;

w2,w2s,w2ss:real; (* Kurven ohne dunkle Materie without dark matter *)

lnW2:real;

lnw2s,lnw2ss:real;

vzw2s:real;

v1,v2:real; (* Bahngeschwindigkeit mit und ohne dunkler Mat. *)

(* trajectory velocity with/without dark matter *)

v1rmax,v2rmax:real;(* Maximalwert von v1 oder v2 *)

(* maximum value of v1 or v2 *)

dmk:real; (* dmk=v1/v2 Dunkle-Materie-Koeffizient *)

(* dark-matter coeffizent *)

graphdriver,graphmode:integer; (* Hilfsgroessen auxillaries *)

ket1:string[1];

vr1,vr2:array[1..500]of real; (* Kurvenpunkte curve points *)

datei,v_curves,w_curves,d_curve:text; (* ASCII-Files *)

Function expmk(lnx:real):real;

(* e hoch x mit Kontrolle exponential function with check *)

begin

if lnx<(-70) then expmk:=0

else expmk:=exp(lnx);

end;

Function f9(x:real):s9;

(* Konvertiert Zahl in Zeichenkette converts number to string *)

var ket11:string[11];

begin

str(x:11,ket11);

if x<0 then f9:='-'+copy(ket11,2,6)+copy(ket11,10,2)

else f9:='+'+copy(ket11,2,6)+copy(ket11,10,2);

end;

Function lnroh(x:real):real;

var y:real;

(* Logarithmus der Massenverteilung

Logarithm of mass distribution *)

begin

y:=sqr(x/(SIGMA*RG));

if y<60

then lnroh:=lnRo-y

else lnroh:=-170.0;

end;

begin

clrscr;

writeln('Der "Dunkle-Materie-Koeffizient d als Verh"ltnis von');

writeln('Membranneigung mit und ohne C,D-Glieder der DGL');

writeln;

writeln('The dark-matter coeffizient d as quotient of the slope');

writeln('of the membrane with and without C-D-members of the equation');

writeln;

(* Dateien oeffenen open files *)

assign(datei, 'darksim5.dat');rewrite(datei);

assign(v_curves,'dsim5-vc.dat');rewrite(v_curves);

assign(w_curves,'dsim5-wc.dat');rewrite(w_curves);

assign(d_curve ,'dsim5-dc.dat');rewrite(d_curve);

writeln('Wo =',Wo);

lnWo:=ln(Wo); (* Raumtiefe am Sonnenrand space depths on sun edge *)

lnRs:=ln(Rs); (* Sonnenradius sun radius *)

lnGs:=ln(Gs); (* Beschleunigung am Sonnenrand acceleration sun edge *)

lnha:=ln(ha); (* 1/2 Schwarzschildradius *)

lnMs:=ln(Ms); (* Sonnenmasse mass off sun *)

lnAe:=lnGs+lnRs-lnWo; (* Aetherbeschleunigung ether acceleration *)

Ae:=expmk(lnAe);

writeln('Ae =',Ae);

lnFo:=lnMs+lnGs-ln(4.0*pi)-2*lnWo; (* Membranspannung memrane tension *)

Fo:=expmk(lnFo);

writeln('Fo =',Fo);

(* A=Ae/Fo A-Koeffizient *)

A:=4*pi*Rs*Wo/Ms;

lnA:=ln(A);

writeln('A =',A);

(* Um die Zahlen im Bereich 1e-30 bis 1e+30 zu halten

To get the numbers between 1e-30 and 1e+30 *)

lnFAKWs:=ln(FAKWs); (* Multiplikator von w' *)

lnFAKWss:=ln(FAKWss); (* Multiplikator von w" *)

RG:=RGLJ*LJ; lnRG:=ln(RG); (* Galaxieradius radius of galaxy *)

RI:=FRI*RG; lnRI:=ln(RI); (* Integrationsradius radius of integr.*)

writeln('RG =',RG);

writeln('RI =',RI);

(* Materiedichte density of matter Galaxiemasse=10^12 Sonnenmassen *)

lnRo:=lnMs+ln(MG)-ln(pi)*1.5-ln(SIGMA)*3-lnRG*3;

Ro:=exp(lnRo);

writeln('Ro =',Ro);

(* Startwerte der DGL initial values of the ODE *)

dr:=RG/MAXSTEP; lndr:=ln(dr); (* Schrittweite stepwidth *)

w1:=W1o;

if w1<0 then begin vzw1:=-1.0; lnw1:=ln(-w1o);end;

if w1=0 then begin vzw1:=0; lnw1:=-70.0;end;

if w1>0 then begin vzw1:=+1.0; lnw1:=ln( w1o);end;

w2:=W2o;

w1s:=0; lnW1s:=-170.0; vzw1s:=1;

w2s:=0; lnW2s:=-170.0; vzw2s:=1;

r:=0;

v1rmax:=0; (* Startwert Maximum Initial value of maximum *)

v2rmax:=0;

i:=0; (* Vektorindex index of vector element *)

readln; (* Pause stopp *)

r_in_Rg:=0;

w1km:=w1/1000;

w2km:=w2/1000;

v1:=0;

v2:=0;

dmk:=1.0;

writeln(' r w1 w2 w1s w2s v1',

' v2 dmk');

writeln(datei,f9(r),' ',f9(w1),' ',f9(w2),' ',f9(w1s),' ',f9(w2s),

' ',f9(v1),' ',f9(v2),' ',f9(dmk));

writeln(v_curves,f9(r_in_Rg),' ',f9(v1),' ',f9(v2));

writeln(w_curves,f9(r_in_Rg),' ',f9(w1km),' ',f9(w2km));

writeln(d_curve ,f9(r_in_Rg),' ',f9(dmk));

repeat

r:=r+dr; lnr:=ln(r);

(* Differentialgleichung differential equation *)

(* Terme von w1ss terms of w1" *)

lnAterm:=lnA+lnroh(r)+lnFAKWss;

(* writeln('lnAterm=',lnAterm); *)

lnCterm :=lnC+lnw1s-lnFAKWs+lnFAKWss;

(* writeln('lnCterm=',lncterm); *)

lnDterm :=lnD+lnW1+lnFAKWss +lnw1s-lnFAKWs;

(* writeln('lnDterm=',lnDterm); *)

w1ss:=-2.0*vzw1s*w1s/r*FAKWss/FAKWs+expmk(lnAterm)

-vzw1s* expmk(2*lnCterm) - vzw1*expmk(lnDterm);

(* writeln('w1ss=',w1ss); *)

w2ss:=-2.0*w2s/r*FAKWss/FAKWs+expmk(lnAterm);

(* writeln('w2ss=',w2ss); *)

(*Formel: w2ss:=-2.0*w2s/r+rohzl(r)+roh(r); *)

(* Euler Integration Eulerian integration *)

w1s:=w1s+w1ss*dr/FAKWss*FAKWs;

if w1s<0 then begin vzw1s:=-1.0; lnw1s:=ln(-w1s);end;

if w1s=0 then begin vzw1s:=0; lnw1s:=-70.0;end;

if w1s>0 then begin vzw1s:=+1.0; lnw1s:=ln( w1s);end;

(* writeln('w1s=',w1s); *)

w2s:=w2s+w2ss*dr/FAKWss*FAKWs;

(* writeln('w2s=',w2s); *)

w1:=w1+w1s*dr/FAKWs;

if w1<0 then begin vzw1:=-1.0; lnw1:=ln(-w1);end;

if w1=0 then begin vzw1:=0; lnw1:=-70.0;end;

if w1>0 then begin vzw1:=+1.0; lnw1:=ln( w1);end;

(* writeln('w1=',w1); *)

w2:=w2+w2s*dr/FAKWs;

(* writeln('w2=',w2); *)

(* Bahngeschwindigkeiten velocities on the trajectory *)

v1:=sqrt(Ae*abs(w1s)*r/FAKWs)/1000;

v2:=sqrt(Ae*abs(w2s)*r/FAKWs)/1000;

(* Dunkler-Materie-Koeffiziet dark-matter coeffizient *)

dmk:=sqr(v1)/sqr(v2);

r_in_Rg:=r/Rg;

w1km:=w1/1000;

w2km:=w2/1000;

(* Immer, wenn r einen Punkt RI/500 nahekommt, Werte speichern

Store value if near a point RI/500

*)

if abs((i+1)*RI/500-r)<(0.6*dr) then

begin

writeln(datei,f9(r),' ',f9(w1),' ',f9(w2),' ',f9(w1s),' ',f9(w2s),

' ',f9(v1),' ',f9(v2),' ',f9(dmk));

writeln(v_curves,f9(r_in_Rg),' ',f9(v1),' ',f9(v2));

writeln(w_curves,f9(r_in_Rg),' ',f9(w1km),' ',f9(w2km));

writeln(d_curve ,f9(r_in_Rg),' ',f9(dmk));

if (i mod 25)=0 then

writeln(f9(r),' ',f9(w1),' ',f9(w2),' ',f9(w1s),' ',f9(w2s),

' ',f9(v1),' ',f9(v2),' ',f9(dmk));

i:=i+1;

vr1[i]:=v1;

vr2[i]:=v2;

if v1>v1rmax then v1rmax:=v1;

if v2>v2rmax then v2rmax:=v2;

end;

until i=500;

readln;

writeln('v1rmax=',v1rmax);

writeln('v2rmax=',v2rmax);

writeln('w1 =',w1);

writeln('w2 =',w2);

readln;

graphdriver:=0;graphmode:=0;

initgraph(graphdriver,graphmode,'c:\tp\lib');

(* L"sche Schirm clean screen *)

setcolor(0);

for i:=1 to 479 do

begin

moveto(0,i);

lineto(639,i);

end;

(* Achsen axis *)

setcolor(2);

outtextxy(200,10,'Geschwindigkeit velocity v1(r) v2(r)');

moveto(100,40);lineto(100,460);

moveto(90,50);lineto(620,50);

moveto(90,450);lineto(620,450);

for i:=1 to 10 do begin moveto(95,40*i+50);lineto(105,40*i+50);end;

for i:=1 to 10 do begin moveto(50*i+100,45);lineto(50*i+100,455);end;

outtextxy(595,35,'RI');

outtextxy(50,130,'v(r)');

outtextxy(520,35,'r');

(* Zeichne Kurve mit Dark-Matter-Effekt gruen

Plott curve with dark-matter effect with green

*)

setcolor(10); outtextxy(10,0,'v1(r) with dark matter');

j:=0;

for i:=1 to 500 do

begin

if j=0 then begin moveto(100+i,450-round(vr1[i]/v1rmax*400.0));j:=1;end;

if j>0 then lineto(100+i,450-round(vr1[i]/v1rmax*400.0));

end;

(* Zeichne Kurve ohne Effekt rot curve without effect with red *)

setcolor(12); outtextxy(10,20,'v2(r) without');

j:=0;

for i:=1 to 500 do

begin

if j=0 then begin moveto(100+i,450-round(vr2[i]/v1rmax*400.0));j:=1;end;

if j>0 then lineto(100+i,450-round(vr2[i]/v1rmax*400.0));

end;

readln;

closegraph;

close(datei);

close(v_curves);

close(w_curves);

close(d_curve);

end.

====================== *** =========================

program haefele_und_keating;

(* geschrieben am written at 23.02.2000 *)

uses crt;

(* Uhr 1 fliegt in 10 km Höhe in Ostrichtung um den Globus. Nach 40 Stunden

ist sie wieder daheim. Diese Uhr bewegt sich mit der Drehrichtung.

Uhr 2 fliegt in Westrichtung, Rest dito.

Clock 1 is flying in altitude 10 Km in direction east round the globe.

After 40 h it returns. This clock moves together with earth rotation.

Clock 2 is flying in direction west, remaining is the same

*)

var n,c,c2,g,Bl,R,h,It,u,v,we,w1,w2,dt1,dt2,dt3,dt,t,Gk:real;

begin

clrscr;

writeln('Uhrenexperiment von Haefele und Keating');

writeln('Clock Experiment of Haefele and Keating');

writeln('=======================================');

writeln;

c:=3e+8; (* Lichtgeschwindigkeit speed of light *)

c2:=sqr(c);

g:=9.81; (* Erdbeschleunigung Acceleration *)

Bl:=36000.0; (* Breitenkreislänge Km Perimeter *)

R:=1000.0*Bl/2.0/pi; (* Breitengradradius m radius *)

h:=10000.0; (* Flughöhe altitude *)

u:=900; (* Fluggeschwindigkeit Kmh speed *)

It:=3600.0*(Bl/u); (* Flugzeit in Sekunden bei u Kmh time *)

v:=500000; (* Erdgeschwindigkeit im Raum speed of earth *)

we:=2*pi/24/3600; (* Winkelgeschwindigkeit der Erde angle speed *)

w1:=we+2*pi/It; (* Winkelgeschwindigkeit Uhr 1 *)

w2:=we-2*pi/It; (* Winkelgeschwindigkeit Uhr 2 *)

(* Berechnung der 3 Integrale integrals

dt1=integral 0.5 w²/c² dt mit w²= (v+R*w1*cos(w1*t)) ² + (R*w1*sin(w1*t)) ²

dt2 dito

dt3=integral 0.5 w²/c² dt mit w²= (v+R*we*cos(we*t)) ² + (R*we*sin(we*t)) ²

*)

n:=1000000.0; (* Zahl der Integrationsschritte number of steps *)

dt:=It/n; (* Zeitdifferential time differential *)

dt1:=0; (* Integralwerte integral values *)

dt2:=0;

dt3:=0;

t:=0; (* Zeit time *)

writeln('c =',c ,' c2=',c2);

writeln('g =',g ,' Bl=',Bl);

writeln('h =',h ,' It=',It);

writeln('u =',u ,' v =',v );

writeln('we=',we,' w1=',w1);

writeln('w2=',w2,' n =',n );

writeln('dt=',dt,' It=',It);

writeln('wait , please');

writeln;

repeat

t:=t+dt;

dt1:=dt1+0.5*(sqr(v+R*w1*cos(w1*t)) + sqr(R*w1*sin(w1*t)))/c2*dt;

dt2:=dt2+0.5*(sqr(v+R*w2*cos(w2*t)) + sqr(R*w2*sin(w2*t)))/c2*dt;

dt3:=dt3+0.5*(sqr(v+R*we*cos(we*t)) + sqr(R*we*sin(we*t)))/c2*dt;

until abs(It-t) < (0.5*dt);

(* Gravitationskorrektur für Uhr1 und Uhr 2: Diese Uhren laufen schneller,

da die Flugzeuge etwa in Hoehe h fliegen

Correction caused by gravitation for clock 1 and 2: This clocks are

running faster, since the aircrafts are flying in altitude h

*)

Gk:=It*h*g/c2; (* Gravitationskorrektur gravitational correction *)

writeln('dt1=',(dt1*1e9):10:5,' dt1-dt3=',((dt1-dt3)*1e9):10:5);

writeln('dt2=',(dt2*1e9):10:5,' dt2-dt3=',((dt2-dt3)*1e9):10:5);

writeln('dt3=',(dt3*1e9):10:5);

writeln('Gk =',(Gk *1e9):10:5);

Writeln;

(* dt1 und dt2 werden negatv erwartet, da die Uhren langsamer gehen

Gk wird positiv erwartet, da die Uhren schneller gehen

dt1 and dt2 are expected to be negative, since the clocks are

runnung more slowly

Gk is expected to be positive

*)

writeln('Uhr1: dt=',((Gk-dt1+dt3)*1e9):10:5);

writeln('Uhr2: dt=',((Gk-dt2+dt3)*1e9):10:5);

writeln('Uhr3: dt=',((dt3)*1e9):10:5);

writeln;

writeln('ENTER');

readln;

end.

============================ *** =========================

program perihel_2;

(* Berechnet Orbit des Merkur mit Zusatzbeschleunigung -K*g*M*a/r^3

durch Integration der originalen Bewegungsgleichungen. Die Genauigkeit

der Integration wird geschätzt.

Große Halbachse A=5.80e10, Kleine B=5.6745e10, g*M=13.245e19

Startwerte v_minus=5.9309e4 r_minus=4.6e10 r_plus=7.0e10

calculates the orbit of Mercury with add. accel. -K*g*M*a/r^3

by integration of the original equations of movement. The

accaracy of integration will be estimated.

Starting values

Great half axis A=5.80e10,

small half axis B=5.6745e10,

Gamma*Mass_of_Sun=13.213e19

speed v_minus=5.9309e4 in the perihelion,

radius r_minus=4.6e10 in the perihelion,

radius r_plus=7.0e10 in the aphelion

Letzte Änderung last change 05/2000

*)

uses crt;

const faktor=0.0; (* Factor of additional accel. -HA*gm/r^3 *)

v_minus=5.8949766e4; (* [m/s] speed in point nearest to sun *)

N=1e7; (* Number of integration steps *)

r_plus=7.0e10; (* [m] largest radius of orbit *)

r_minus=46.0e9; (* [m] smallest radius of orbit *)

gm=13.245e19; (* [m^3/s²] Gravitat.Const times mass of sun *)

p=55.517e9; (* [m] Semilatus rectum of Ellipsa *)

eps=0.2068966; (* Numerical Excentricity Mercury *)

c=2.99792e8; (* [m/s] speed of light *)

Mm=0.33021e24; (* [Kg] Mass of Mercury *)

Ms=1.9851e30; (* [Kg] Mass of Sun *)

Tr=87.9; (* [d] Time of revolution Mercury *)

var i :integer;

x,y :real; (* coordinates from focal point *)

vx,vy :real; (* speeds *)

vr,vphi :real;

r,rq :real; (* radii *)

sinphi,cosphi :real; (* trigonometrics *)

dt,t,Ts :real; (* time *)

a,ax,ay :real; (* accelerations *)

phi,phi0 :real; (* angles *)

phimi,phipl :real;

phi_print :real;

xmi,xpl,ymi,ypl :real; (* zero positions *)

vxmi,vxpl,vymi,vypl:real; (* zero speeds *)

vrmi,vrpl :real;

phivrmi,phivrpl :real; (* zero angles *)

stop,first,nearend :boolean; (* for control *)

dphi,x1,y1,vx1,vy1 :real; (* values after 1 revolution *)

phi1,x180 :real;

error :real;

HA :real; (* Half Schwarzschild radius *)

r1,v1 :real; (* radius and speed after 1 revolution *)

begin

clrscr;

writeln('Calculation orbit of Mercury');

writeln;

x:=r_minus;

y:=0;

vx:=0;

vy:=v_minus;

HA:=gm/sqr(c);

Ts:=87.9*24.0*3600.0;

dt:=Ts/N;

t:=0;

phi_print:=0;

first:=true;

nearend:=false;

x180:=0;

xpl:=0;

vrpl:=0;

repeat

rq:=sqr(x)+sqr(y);

r:=sqrt(rq);

a:=-gm/rq-faktor*HA/r*gm/rq;

cosphi:=x/r;

sinphi:=y/r;

ax:=a*cosphi;

ay:=a*sinphi;

vx:=vx+ax*dt;

vy:=vy+ay*dt;

x:=x+vx*dt;

y:=y+vy*dt;

phi:=arctan(y/x);

if x<0 then phi:=phi+pi;

if (x>0) and (y<0) then phi:=phi+2*pi;

vr := cosphi*vx+sinphi*vy;

vphi:=-sinphi*vx+cosphi*vy;

t:=t+dt;

if phi>=phi_print then

begin

if first then

begin

writeln(' Phi x y r vphi vr');

writeln('-----------------------------------------------------------');

first:=false;

end;

writeln(phi:10:5,

(x*1e-9):10:5,

(y*1e-9):10:5,

(r*1e-9):10:5,

vphi:10:1,

vr:10:1);

phi_print:=phi_print+2*pi/20;

end;

if (phi>=pi) and (x180=0) then x180:=x;

if (phi>6.2) and not nearend then nearend:=true;

if nearend then

begin

phi0:=sinphi;

if vr<0 then

begin

phivrmi:=phi0; vrmi:=vr; xmi:=x; ymi:=y; vxmi:=vx; vymi:=vy;

end;

if (vr>0) and (vrpl=0) then

begin

phivrpl:=phi0; vrpl:=vr; xpl:=x; ypl:=y; vxpl:=vx; vypl:=vy;

end;

end;

until vrpl>0;

write('ENTER');readln;

(* Estimation of r, v, phi (after 1 revolution) *)

dphi:=phipl-phimi; (* distance of 2 phi-values *)

x1:=xmi-(xpl-xmi)/dphi*phimi; (* since phimi<0 *)

y1:=ymi-(ypl-ymi)/dphi*phimi;

r1:=sqrt(sqr(x1)+sqr(y1));

vx1:=vxmi-(vxpl-vxmi)/dphi*phimi;

vy1:=vymi-(vypl-vymi)/dphi*phimi;

v1:=sqrt(sqr(vx1)+sqr(vy1));

phi1:=phivrmi-(phivrpl-phivrmi)/(vrpl-vrmi)*vrmi;

writeln('r :',(r_minus*1e-9):20:10, (r1*1e-9):20:10);

writeln('v :',v_minus:20:10,v1:20:10);

writeln('x180 :',(x180*1e-9):20:10);

writeln('PHI :',(0.0):20:10, phi1:20:10);

writeln('Faktor :',faktor:20:10);

writeln('Tu :',(t/24/3600):20:10);

error:=sqrt(sqr((r_minus-r1)/r_minus)+

sqr((v_minus-v1)/v_minus));

writeln('Error :',error:10);

write('ENTER');readln;

end.

============================ *** ====================

program berechng;

(* Berechnet verschiedene Werte aus Konstanten *)

(* Calculates different values from constants *)

uses crt;

const

(* Allgemeine Konstanten *)

c=2.99792458e8; (* m/s Lichtgeschwindigkeit speed of light*)

Qel=1.60217646263e-19; (* As Elektronenladung charge of electron*)

eps0=8.8541878e-12; (* F/m Dielektrizitaetskonstante *)

Gamma=6.67310e-11; (* m^3/(Kg*s*s) Gravitationskonstante *)

h=6.6260687652e-34; (* Js Planck-Konstante *)

lambda_C=2.426310215e-12; (* m Comptonwellenl"nge Elektron *)

Melek=9.1093818872e-31; (* Kg Elektronenmasse mass of electron*)

Mprot=1.6726215813e-27; (* Kg Protonenmasse mass of proton *)

Ho=75; (* (km/s)/Mpc Hubble-Konstante *)

OmegaZ=1.5527e21; (* 1/s Kreisfrequenz Zitterbewegung *)

magneton=9.28477e-24; (* Am^2 Bohr's Magneton *)

(* Zeit *)

dsun=86400.0; (* s Sonnentag mean sun-day*)

dstar=86164.2; (* s Sterntag siderischer Tag *)

Jtrop=365.2421897; (* dsun Tropisches Jahr mean tropical year *)

AstrEinh=1.4959787066e11; (* m Erdbahnradius Astronomical Unit *)

LJ=9.460730472e15; (* m Lichtjahr light year *)

pc=3.0856776e16; (* m Parsec *)

(* Erde *)

Raequ=6.3781e6; (* m Erdradius Aequator *)

Rpole=6.3568e6; (* m Erdradius Pole *)

Rerde=6.3710e6; (* m Erdkugelradius radius of sphere *)

Merde=5.9736e24; (* Kg Erdmasse *)

Gaequ=9.78031; (* m/s*s Fallbeschleunigung Aequator *)

Gpole=9.83217; (* m/s*s Fallbeschleunigung Pole *)

Rperi=1.4709e11; (* Perihel Erdbahn um COM (Center of Mass) *)

Raphe=1.5210e11; (* Aphel Erdbahn um COM *)

Tside=365.256; (* dsun Siderische Umlaufszeit orbital time *)

Ttrop=365.242; (* dsun Tropische Umlaufszeit *)

Verde=29.78e3; (* m/s Mittlere Umlaufgeschwindigkeit Erde *)

(* Sonne *)

Msonne=1.989e30; (* Kg Sonnenmasse mass of sun *)

Rsonne=6.9550826e8; (* m Sonnenradius radius of sun *)

var Gsonne, rc, rohc, Fo, Wo :real;

a, a3, Rohsonne, WsR, Ae: real;

RCOM, Rb, DAQ, DWs :real;

begin

clrscr;

writeln('Zuerst aus Comptonwellenl"nge Lambda_C die Raumtiefe Wo berechnen');

writeln('First from Compton wave-length depth of space Wo at edge of sun');

writeln;

Gsonne:=Gamma*Msonne/Rsonne/Rsonne;

writeln(Gsonne:15, ' m/s*s Beschleunigung am Sonnenrand acceleration');

(* Comptonradius *)

rc:=lambda_c/pi;

writeln(rc:15, ' m Comptons Elektronenradius Lambda/pi');

writeln;

(* Membrandichte aus Compton-Radius*)

rohc:=(Qel/(16*pi*pi*c*rc*rc))*(Qel/(rc*rc*c*eps0));

writeln(rohc:15, ' Kg/m^3 Density RohC from Comptons Rc calculated');

(* Membranspannung *)

Fo:=rohc*sqr(c);

writeln(Fo:15, ' N/m*m Membrane tension from Compton calculated');

(* Raumtiefe *)

Wo:=sqrt( Msonne*Gsonne/(4*pi*rohc*sqr(c)) );

writeln(Wo:15, ' m Depth of space from Compton calculated');

writeln;

(* 1/2 Schwarzschildradius=G*M/(c*c) *)

a:=Gamma*Msonne/c/c; a3:=a/3;

writeln(a:15, ' m 1/2 Schwarzschildradius der Sonne');

writeln(a3:15, ' m Feynmans Excessradius der Sonne');

(* Mittlere Sonnendichte *)

Rohsonne:=3*Msonne/4/pi/Rsonne/Rsonne/Rsonne;

writeln(Rohsonne:15, ' Kg/m^3 Mean density of sun');

writeln;

wsR:=Wo/Rsonne;

writeln(WsR:15, ' [ ] Slope of membrane at edge of sun');

Ae:=Gamma*Msonne/Wo/Rsonne;

writeln(Ae:15, ' m/s*s Aether acceleration');

writeln;

writeln(' Einige Formeln zu Gravitationswellen');

writeln('Some formulas to gravitational waves');

writeln;

RCOM:=Merde*AstrEinh/Msonne;

writeln(RCOM:15, ' m Distance sun to COM (Center of Mass)');

Rb:=c*Jtrop*dstar/4/pi;

writeln(Rb:15, ' m Distance Rb of point of balance Pb to COM');

writeln;

DAQ:=Gamma*Merde/sqr(Rb)*sqr(AstrEinh)*3/sqr(Rb);

writeln(DAQ:15, ' m/s^2 Acceleration difference at Pb');

DWs:=DAQ/Ae;

writeln(DWs:15, ' [ ] Change of slope of membrane at Pb');

write('ENTER');readln;

readln;

end.

============================ *** ====================

program ro_d3;

(* Berechnet fuer d=3 und frei waehlbare Dichtefunktion ro(r) die

Wegverlaengerung dS

Calculates for d=3 and free choosable density function ro(r) the

path lengthening *)

uses crt;

const N=1000;

Wo=1.4317e6;

M=1.989e30; (* Kg Sonnenmasse mass of sun *)

G=6.67310e-11; (* m^3/(Kg*s*s) Gravitationskonstante *)

Rs=6.9550826e8; (* m Sonnenradius radius of sun *)

c=2.99792458e8; (* m/s Lichtgeschwindigkeit speed of light*)

Eps1=1e-5;

gamma=0.001;

var i:integer;

r,rq,dr,Fo,Fe,wsR,Srorq,roo:real;

tanalfa,tanalfa0,tanbeta,dFo:real;

dSr,dSu:real;

ro,w,ws:array[0..N]of real;

begin

clrscr;

writeln('Berechnet fuer d=3 und frei waehlbare Dichtefunktion');

writeln('ro(r) die Wegverlaengerung dS');

writeln;

(* Eingangskonstante *)

dr:=Rs/(N+0.5);

wsR:=Wo/Rs;

Fe:=G*M/Wo/Rs;

Fo:=M/4/pi/Wo/Rs*Fe;

writeln('Schrittweite dr ',dr:12);

writeln('Startneigung wsR ',wsR:12);

writeln('Aetherbeschleunigung Fe ',Fe:12);

writeln('Membranspannung Fo ',Fo:12);

write('ENTER');readln;

(* Berechne Dichtekurve und das Integral density curve and integral *)

Srorq:=0;

for i:=0 to N do

begin

r:=0.5*dr+i*dr;

ro[i]:= 1407.0; (* Konstante mittlere Sonnendichte 1407.667 mean density*)

(*

( 3.517e3 * (1.0-sqr(r/Rs))

+2.462e3 * (1.0-sqr(sqr(r/Rs)))

+3.889e3 * cos(pi*r/Rs/2)

+2.537e3 * sqrt(abs(cos(pi*r/Rs/2))) )/4;

*)

Srorq:=Srorq+ro[i]*sqr(r);

end;

roo:=M/4/pi/dr/Srorq;

writeln('roo ',roo:12);

write('ENTER');readln;

for i:=0 to N do

begin

r:=0.5*dr+i*dr;

ro[i]:=1407; (* roo*ro[i]; *)

end;

(* Zur Kontrolle die Masse *)

Srorq:=0;

for i:=0 to N do

begin

r:=0.5*dr+i*dr;

Srorq:=Srorq+4.0*pi*ro[i]*sqr(r)*dr;

end;

writeln('Masse ',Srorq:12);

write('ENTER');readln;

writeln;

writeln('r,ro(r)');

for i:=0 to N do

begin

r:=0.5*dr+i*dr;

if (i mod 10)=5 then writeln('r=',r:12,' ro(r)=',ro[i]:12);

end;

write('ENTER');readln;

(* Berechne Raumtiefenkurve. Fo iterativ ver"ndert, bis Winkel bei

r=0 zu Null wird. *)

repeat

w[N] :=-Wo;

ws[N]:=wsR;

for i:=N-1 downto 1 do

begin

r:=0.5*dr+i*dr;

rq:=sqr(r);

(* Berechne die notwendigen Aenderungen von w fuer Kraefteausgleich.

Im Punkt N ist z.B. vorgegeben: ws[N]=wsR. Daraus folgt

ws[N-1]. *)

if i=(N-1) then tanbeta:=wsR else tanbeta:=tanalfa;

tanalfa:=ws[i+1]*sqr(r+dr)/rq - Fe/Fo*dr*ro[i]*sqr((r+dr/2)/r);

if i=1 then tanalfa0:=tanalfa;

w[i] :=w[i+1]-tanalfa*dr;

ws[i]:=tanalfa;

if (i mod 10)=5 then writeln(

'r=',r:12,' w(r)=',w[i]:12,' ws(r)=',ws[i]:12);

end;

dFo:=tanalfa0*Fo;

Fo:=Fo-dFo*gamma;

writeln;

writeln('tanalfa0=',tanalfa:12,' dFo=',dFo:12,' Fo=',Fo:12);

until abs(tanalfa0)<(Eps1*abs(wsR));

(* Berechnung der beiden Wegverl"ngerungen *)

dSr:=0;

for i:=1 to N do dSr:=dSr+sqr(ws[i]);

dSr:=dSr*0.5*dr;

dSu:=0;

for i:=30000 downto 0 do

begin

r:=Rs+i*dr;

wsR:=Wo*Rs/r/r;

dSu:=dSu+sqr(wsr);

end;

dSu:=dSu*0.5*dr;

writeln;

writeln('Wo=',Wo:15,' innen dSr=',dSr:15,' aussen dSu=',dSu:15);

writeln('W1=',w[1]:15);

readln;

end.

============================ *** =============================

============================= *** =============================

============================= *** =============================

============================= *** =============================

============================= *** =============================