The Turbo
Pascal- or C- source codes of most of the calculations and simulations in
relation with the membrane theory of gravity.
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
Membran4DView ……………… Calculation of the
curved 3D membrane in the 4D bulk space
=============================
*** =============================
program Trouton_Noble_3;
uses crt,dos;
(* Zentriertes Ellipsoidmodell vom 1.2.2000:
Berechnet die Energie des elektromagnetischen Feldes U=Integral E2±B2 dV
zweier Ladungen im Abstand d=10 in Abhängigkeit von b=v/c;
Centered ellipsoid model:
Calculates the energy of
electromagnetical field U=integral E2-B2 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 E2±B2 dV');
writeln('===========================================================');
writeln;
writeln(' b Ue(b) Ub(b) ');
writeln(' __________________________________');
writeln;
beta:=0.0; (*b =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 2 *)
wembq:=sqrt(embq); (* root
of 1- b 2 *)
(* 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;
(* r2 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)/r2. 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)/r2. 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 E2
and B2 *)
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 E2 *)
Ubi[shell]:=Ubi[shell]+B2;
(* Integration of B2 *)
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 w2/c2 dt mit w2=
(v+R*w1*cos(w1*t)) 2 + (R*w1*sin(w1*t)) 2
dt2 dito
dt3=integral 0.5 w2/c2
dt mit w2= (v+R*we*cos(we*t)) 2 + (R*we*sin(we*t)) 2
*)
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/s2]
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
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
(* Raumtiefe *)
Wo:=sqrt( Msonne*Gsonne/(4*pi*rohc*sqr(c)) );
writeln(Wo:15, ' m Depth of
space from
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.
============================
*** =============================
program energy;
(* Calculates the energy integral surrounding a moved
source
considering the Doppler effect.
We consider an uniform distribution of photons of the
same
wavelenght over the spatial angle in a resting
mirroring sphere.
No we consider the same objekt moving with speed v.
The number of
photons travelling in x direction is greater than the
number of
those travelling in -x direction because of the speed
v:
n = n+ + n- or n=np+nm with np=n*b/(c-v), nm=n*b/(c+v)
mean energy is E(alfa=0) = b/(c-v) * C/lambdap +
b(c+v) * c/lambdam
oder E = b*nue0(
1/(c-v)^2 + 1/(c+v)^2 )
*)
uses crt;
const N=10000.0; (* Number of integration steps *)
var i, jbeta:integer; (* beta = v / c *)
beta, b, a, integral :real; (* description see below
*)
d_alfa_s, alfa_s, sinas :real;
y_s, x_s, rsq, r_s, bx, r_ss : real;
abq, gamma :real;
bb, alfa_ss :real;
DE, fak, dilat, q : real;
outfile:text;
begin
clrscr;
assign(outfile,'dilation.prt');
rewrite(outfile);
d_alfa_s:=pi/N; (* Step for angle *)
for jbeta:=0 to 19 do
begin
beta:=0+jbeta*0.05;
b := 1-sqr(beta); (* Short axis in x direction *)
a := sqrt(b); (* longer axis in y direction *)
abq:=a/sqr(b);
integral:=0; (* Starting value *)
alfa_s:=d_alfa_s/2; (* Starting angle between x and
r_s *)
repeat (* Integrate for alfa from 0 to pi *)
sinas:=sin(alfa_s);
y_s:=a*sinas; (* y' or altitude of point P *)
x_s:=b*cos(alfa_s); (* x'-coordinate of P *)
rsq:= sqr(y_s)+ sqr(x_s); (* r' squared *)
r_s:=sqrt(rsq); (* Distance A-P or r' *)
bx:=beta*x_s / b; (* dummy for r'' calculation *)
r_ss:=bx + sqrt( sqr(bx) + rsq/b ); (* Distance S-P or
r" *)
bb:=beta*r_ss + x_s; (* Distance S-Px *)
alfa_ss:=arctan(y_s/bb); (* angle between r" and
x *)
gamma:=arctan(-x_s*abq/sqrt(1-sqr(x_s/b)))+pi/2;
if bb<0 then alfa_ss:=alfa_ss+pi; (* lack of this
arctan function *)
(* DE:=1.0/(1.0-cos(alfa_ss)*beta); Doppler effect *)
(* Integration considering radius y' of the ring,
Doppler effect
and visible width of the ring
writeln('beta =',beta);
writeln('alfa_s =',alfa_s);
writeln('alfa_ss
=',alfa_ss);
writeln('gamma
=',gamma);
writeln('r_s =',r_s);
writeln('r_ss =',r_ss);
readln;
*)
fak:= sin(alfa_s)/sqr(r_ss);
integral:=integral + fak * d_alfa_s;
alfa_s:=alfa_s+d_alfa_s; (* increase angle alfa' *)
until (alfa_s>pi); (* Test end of integration *)
integral:=integral/2-1; (* Norm to 1 and subtract 1 *)
dilat:=1/sqrt(1-sqr(beta))-1; (* SR time dilation
minus 1 *)
if jbeta=0 then
begin
writeln('beta=',beta:10:6,' Int=',integral:15:10,
' Dil=',dilat:15:10);
writeln(outfile,
'beta=',beta:10:6,' Int=',integral:15:10,
' Dil=',dilat:15:10);
end
else
begin
q:=integral/dilat; (* ratio *)
writeln('beta=',beta:10:6,'
Int=',integral:15:10,
' Dil=',dilat:15:10,' q=',q:10:7);
writeln(outfile,
'beta=',beta:10:6,' Int=',integral:15:10,
' Dil=',dilat:15:10,' q=',q:10:7);
end;
end;
close(outfile);
readln;
end.
=============================
*** =============================
// Program Ioncryc.cpp calculates the potential of a
central
// ion in a lattice of positive and negative ions
(NaCl-type
// crystal). Lattice constant is 2, i.e. we find the
same kind
// of ions in the nearest distance d=2. The other sort
of ions
// is placed in the center of a cube built from 8 ions
of the
// first kind.
// Programmer: Stefan von Weber, May 2004
#include "Stdlib.h"
#include "Math.h"
#include "Stdio.h"
// N = Greatest radius from central ion
// A defines thickness of shells by dr=1/A
// AN = Number of shells
const int N=500, A=150, AN=A*N, AN1=AN+1;
int i,j,k,sph;
double dx,dyz,beta,one_minus_beta_squared,pot,
// sphere[i] is a shell surrounding the central ion
double sphere[AN1];
void add_pot()
// Adds the contribution of one ion to the potential
// of the central ion. The contributions are collected
// in different shells
{int spherenum; // Shellnumber
int odds; // Even positions or odd positions
double r,rs,x,y,z,xx,yy,zz;
odds=0;
if ((abs(i) % 2)==1){ odds++;}
if ((abs(j) % 2)==1){ odds++;}
if ((abs(k) % 2)==1){ odds++;}
if (odds==0) // an ion of the same kind as the central
{
x=dx*i; xx=i; // Contracted / not contracted
coordinates
y=dyz*j; yy=j;
z=dyz*k; zz=k;
// Euclidean and Retarded distance
r=sqrt(xx*xx + yy*yy + zz*zz);
rs=sqrt(x*x + one_minus_beta_squared*(y*y + z*z));
if (r>0)
{ spherenum=floor(r*A); // Number of shell
pot=one_minus_beta_squared/rs; // According equ. (4)
if (spherenum<=AN)
{ sphere[spherenum]=sphere[spherenum]+pot;}
}
}; // odds==0
if (odds==3) // the other kind of ions
{
x=dx*i; xx=i;
y=dyz*j; yy=j;
z=dyz*k; zz=k;
// Euclidean and Retarded distance
r=sqrt(xx*xx + yy*yy + zz*zz);
rs=sqrt(x*x + one_minus_beta_squared*(y*y + z*z));
if (r>0)
{ spherenum=floor(r*A);
pot=one_minus_beta_squared/rs;
if (spherenum<=AN)
{ sphere[spherenum]=sphere[spherenum]-pot;} //
subtract
}
}; // odds==3
return;
}
void main()
{
//program ioncryst;
// Calculates potential of NaCl-type ion lattice in a
moved frame
// using the cross and length contraction dependent on
speed v
// at the origin.
// For comparation additional calculations are done
with other kinds of
// contractions: No contraction, Lorentz-contraction,
Cross- and Length
// i runs in x-axis (parallel to speed v)
// j in y-axis
// k in z-axis
// positive ions we find at even integer positions
// negative ions we find at odd integer positions
int jbeta;
double P1,P2,P3;
FILE* aus;
aus=fopen("D:\\p\\waves\\ioncrya.dat","w");
fprintf(aus,"Potential of an ionic crystall of
NaCl type");
for (jbeta=0;jbeta<7;jbeta++)
{
beta=0.001*jbeta; // beta=v/c
one_minus_beta_squared=1.0-beta*beta; // 1-beta^2
// No contraction
dx =1.0; // Lattice distances
dyz=1.0;
// Starting values of the sums
for(sph=0;sph<=AN;sph++) {sphere[sph]=0.0;}
// The potential in the single spheres
for( i=-N;i<=N;i++)
{
printf("\nP1 beta=%6.3f i=%i von
%i",beta,i,N);
for (j=-N;j<=N;j++)
{ for( k=-N;k<=N;k++){ add_pot();}
}
};
// Add potential of the spheres
// Damping with exponential factor similarly to
Laplace transform
// Because of error minimization go from right to left
P1=0.0;
for(sph=AN;sph>=0;sph--)
{
P1=P1+ sphere[sph]*exp(-sph*11.5/AN);
printf("\nbeta=%6.4f i=%i p=%10.5f",
beta, sph, P1);
};
if(jbeta==0) Po=P1; // Store the zero case
P2=P3=0.0;
if(jbeta!=0)
{
// only Lorentz contraction
dx=sqrt(one_minus_beta_squared);
dyz=1.0;
for(sph=0;sph<=AN;sph++)
{sphere[sph]=0.0;}
// The potential in the single spheres
for( i=-N;i<=N;i++)
{
printf("\nP2 beta=%6.3f %5i von
%i",beta,i,N);
for (j=-N;j<=N;j++)
{ for( k=-N;k<=N;k++){ add_pot();}
}
};
// Add potential of the spheres
// Because of error minimization go from right to left
P2=0.0;
for(sph=AN;sph>=0;sph--)
{
P2=P2+ sphere[sph]*exp(-sph*11.5/AN);
printf("\nbeta=%6.4f i=%i p=%10.5f",
beta, sph, P2);
};
// Cross and lengt contraction
dx=one_minus_beta_squared;
dyz=sqrt(one_minus_beta_squared);
for(sph=0;sph<=AN;sph++)
{sphere[sph]=0.0;}
// The potential in the single spheres
for( i=-N;i<=N;i++)
{
printf("\n P3 beta=%6.3f %5i von
%i",beta,i,N);
for (j=-N;j<=N;j++)
{ for( k=-N;k<=N;k++){ add_pot();}
}
};
// Add potential of the spheres
// Because of error minimization go from right to left
P3=0.0;
for(sph=AN;sph>=0;sph--)
{
P3=P3+ sphere[sph]*exp(-sph*11.5/AN);
printf("\nbeta=%6.4f i=%i p=%10.5f",
beta, sph, P3);
}
}; // Case jbeta!=0
// Output to file
fprintf(aus,
"\n beta=%6.3f P1=%15.10f P2=%15.10f P3=%15.10f
",
beta,P1,P2,P3);
fprintf(aus,
"\n beta=%6.3f p1=%15.10f p2=%15.10f p3=%15.10f
\n",
beta,P1/Po,P2/Po,P3/Po);
} //Schleife jbeta
printf("\n\n");
fclose(aus);
}
=============================
*** =============================
// Program Pioneer : Accelerations by different
potentials
//
// started 12/2005
// Programmer
Stefan von Weber
//
// The normal acceleration is -Gamma*Ms/r^2
//
// Gerber-Einstein acceleration is -Gamma*Ms*6*a/r^3
//
// The correction term is automatically calculated
// (by the NASA) if the least square algorithmus fits
the planet
// data to the general model
//
// If the spacecraft is leaving the planet belt the
fit does not
// further hold
#include <Iostream.h>
#include <Math.h>
#include <Stdio.h>
// the correction_term controls the minimum
// of the sum of squares of the errors
#define correction_term -1.9e-10 // -10.5e-10 for
summation until 20 AU
// -3.2e-10e-10
until 100 AU
// -1.9e-10 until
200 AU
#define Rs
6.9635e+8 // 6.9635e8 Sonnenradius sun radius
#define gamma
6.6726e-11 // gravitational constant
#define ha 1473.1
// Halber Schwarzschildradius 1/2 Schw-Schildr
#define Ms 1.991e+30 // Sonnenmasse Mass of Sun sun
mass
#define SUMRAD 200.0 // Radius of summation in AU
void main()
{
// The radii of the planets in AU (1.5*10^11 m), and
then
// other error measuring positions freely defined
// Merc Ven Earth Mars Jup Sat Ur
double radii[40]={0.387, 0.723, 1.0, 1.523, 5.202,
9.555, 19.278,
25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0,
65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 99.5,
125.0, 130.0, 135.0, 140.0, 145.0, 150.0, 155.0,
160.0,
165.0, 170.0, 175.0, 180.0, 185.0, 190.0, 195.0, 199.5,
500.0};
double error, errorqsum;
double r,r_in_AU;
// Radius m AU
int j; // counter
double a1; // acceleration with additional term
double a2; // acceleration without add. term
char zeile[200];
j=0; // planet position counter
errorqsum=0.0; //sum of squares of the errors
do
{ r_in_AU=radii[j];
r=r_in_AU*1.5e11; // r from astronomical units to
meter
a2=gamma*Ms/r/r;
a1=a2*(1.0+6.0*ha/r);
error= a2*6.0*ha/r + correction_term;
errorqsum+=error*error;
if(j==0)
{ cout<<"\n r [AU] a1(r) a2(r) error";
}
sprintf(zeile,"\n%8.4f
%15e %15e %15e",
r_in_AU,a1,a2,error);
cout<<zeile;
j++; // next planet or position
}while (radii[j]<SUMRAD);
cout<<"\n
Errorqsum="<<errorqsum<<" ";
}
=============================
*** =============================
program darksim9;
uses crt,graph;
(* Berechnet die
Biegung der Membran fr folgendes Modell:
Computes the curvature of the membrane with the
following model:
w"=-2w'/r + A*ro(r) + C*(w'^2);
Factor "C" means a force in opposite of
w-direction
*)
const
lnC=-13.002; (* -13.002 ln of coefficient of w'w' *)
W1o=-0.95e6; (*
-0.95e+6 Startwert w1(r) initial value *)
W2o=-0.593e6; (*
-0.593e6 Startwert w2(r) initial value *)
FRI=2.0; (* 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 fr ws factor for w' *)
FAKWss=1e30; (* Faktor fr wss factor for w" *)
MAXSTEP=250000.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; (*
Žtherbeschleunigung 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;
lnEterm:real;
i,j:integer; (* Hilfsgr"áen 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; (* Hilfsgr"áen
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
"ffenen 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;
(* Žtherbeschleunigung 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+2*(lnw1s-lnFAKWs)+lnFAKWss;
(* writeln('lnCterm=',lncterm); *)
w1ss:=-2.0*vzw1s*w1s/r*FAKWss/FAKWs+expmk(lnAterm)
+expmk(lnCterm);
(* 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 grn
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.
=============================
*** =============================
Calculation of the curved 3D
membrane in the 4D bulk space
// Membrane4DView.cpp : Implementierung der Klasse CMembrane4DView
//
#include
"stdafx.h"
#include
"Membrane4D.h"
#include
"Membrane4DDoc.h"
#include
"Membrane4DView.h"
#ifdef
_DEBUG
#define
new DEBUG_NEW
#undef
THIS_FILE
static
char THIS_FILE[] = __FILE__;
#endif
/////////////////////////////////////////////////////////////////////////////
//
CMembrane4DView
IMPLEMENT_DYNCREATE(CMembrane4DView,
CView)
BEGIN_MESSAGE_MAP(CMembrane4DView,
CView)
//{{AFX_MSG_MAP(CMembrane4DView)
// HINWEIS - Hier werden Mapping-Makros vom Klassen-Assistenten eingefügt und entfernt.
// Innerhalb dieser generierten Quelltextabschnitte NICHTS VERÄNDERN!
//}}AFX_MSG_MAP
// Standard-Druckbefehle
ON_COMMAND(ID_FILE_PRINT,
CView::OnFilePrint)
ON_COMMAND(ID_FILE_PRINT_DIRECT,
CView::OnFilePrint)
ON_COMMAND(ID_FILE_PRINT_PREVIEW,
CView::OnFilePrintPreview)
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CMembrane4DView Konstruktion/Destruktion
CMembrane4DView::CMembrane4DView()
{
// ZU ERLEDIGEN: Hier Code zur Konstruktion einfügen,
}
CMembrane4DView::~CMembrane4DView()
{
}
BOOL
CMembrane4DView::PreCreateWindow(CREATESTRUCT& cs)
{
// ZU ERLEDIGEN: Ändern Sie hier die Fensterklasse oder das Erscheinungsbild, indem Sie
//
CREATESTRUCT cs modifizieren.
return CView::PreCreateWindow(cs);
}
#include
"Math.h"
void
GaussJordan(int n, double A[4][4], double b[4], double x[4], double *PP)
// Loesung eines linearen Gleichungssystem Ax=b bis n=4
// Matrizen muessen als [4][4]-Matrizen und Vektoren als [4]-Vektoren deklariert werden
{ int i, ii, j, jj, ip, jp, k, rx[4], idummy;
double B[4][4], xx[4], F, Piv, Pivv[4],
dummy;
// Umspeichern A nach B und b nach xx um A und b unverändert zu erhalten
for(i=0; i<n; i++)
{ for(j=0;
j<n; j++) B[i][j] = A[i][j]; // Koeffizientenmatrix
xx[i] = b[i]; // Rechte Seite
rx[i]=i; // Reihenfolge der x-Elemente
}//for i
for(i=0; i<n; i++)
{ // Suche groessten Wert ip, jp im Rest der Matrix
Piv=-1e25;
for(ii=i; ii<n; ii++)
{ for(jj=i;
jj<n; jj++)
{ if(Piv<fabs(B[ii][jj])) { Piv=fabs(B[ii][jj]); ip=ii; jp=jj;
} }
} //for ii
// Vertausche Zeilen i und ip einschliesslich Rechte Seite
for(jj=0; jj<n; jj++) { dummy=B[ip][jj]; B[ip][jj]=B[i][jj]; B[i][jj]=dummy; }
dummy = xx[ip]; xx[ip]=xx[i]; xx[i]=dummy;
// Vertausche Spalten i und jp einschliesslich x-Reihenfolge
for(ii=0; ii<n; ii++) { dummy=B[ii][i]; B[ii][i]=B[ii][jp]; B[ii][jp]=dummy; }
idummy = rx[i]; rx[i]=rx[jp]; rx[jp]=idummy;
Piv = B[i][i]; // Pivotwert im i-ten Schritt
Pivv[i]=Piv; // Merke alle Pivotwerte
// Stelle Null her, indem von einer Zeile ii die Pivotzeile i subtrahiert wird
for(ii=0; ii<n; ii++)
{ if(ii!=i) { F
= B[ii][i]/Piv; // Faktor
for(k=i+1;
k<n; k++) B[ii][k] -= F*B[i][k];
xx[ii] -= F*xx[i]; // Dasselbe für die Rechte Seite
}//if ii
}//for ii
// Division der Pivotzeile
for(j=i+1; j<n; j++) B[i][j] /=
Piv;
xx[i] /= Piv;
}// for i
// Umspeichern der Lösung xx nach x, dabei Reihenfolge wieder herstellen
for(i=0; i<n; i++) x[ rx[i] ] = xx[i];
// n-ter Pivot durch 1. Pivot ist das Gütemaß PP der Koeffizientenmatrix
*PP = fabs(Pivv[n-1] / Pivv[0]);
return;
}
/////////////////////////////////////////////////////////////////////////////
//
CMembrane4DView Zeichnen
#include
"Stdio.h"
#include
"Malloc.h"
#include
"string.h"
#define N 101
#define MAXSTEP1 1000 // Iterationsschritte für die anfängliche Membrankrümmung
#define MAXSTEP2 10000 // Iterationsschritte für die Planetenbewegung (Zeitschritte)
#define MAXSTEP3 5 // Iterationsschritte für die Umgebung des Planeten innerhalb eines Zeitschrittes
#define
PI 3.14159265
void
CMembrane4DView::OnDraw(CDC* pDC)
{
CMembrane4DDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
// ZU ERLEDIGEN: Hier Code zum Zeichnen der ursprünglichen Daten hinzufügen
char zeile[200];
FILE *ein, *aus, *aus1;
int i,j,k,m,n, Nh, N2, N3, N3max, N3max1, N3h;
int ii, jj, kk, ix,iix,jy,jjy, ind, ind0, ind1, ind2, ind3, tape;
int nj, nk, sortend, sortstep, num, mm;
int step, step0, step0file, step1, kbz;
int i00, i11, i22, i33, j00, j11, j22,
j33, k00,k11,k22,k33, gefunden;
int ktp0, ktp1, ktp2, ktp3, ktpv[4];
int Rox, Roy,Rux, Ruy; //Graphikgrenzen
int ibz, jbz, anz, iv[500], jv[500], kv[500];
int step3;
int pnr, ii00, jj00, kk00;
double dx, dxx,dyy, dzz, dww, Lo, Lp, dt,
Fo, Mp;
double dh, dh3, ht, dxh, xx, yy, zz, ww,
radius, Rmax, Kp, Fx, Fy, Fz, Fw;
double oldrad, radsi, rad0, rad1;
double W1, W2, W2h, Rmaxh, Rmax95;
double Ee, AltEe, Ediff, epsilon, Fo2dx,
rq, wdx, mindist, dist;
double x0, y0, z0, v0x,v0y,v0z;
double dtMp, eps1, vxx, vyy, vzz, xpl,
ypl, zpl;
double rqmin, rqv[150], r02, r12, r03,
r13, r23, Lpv[4], A[4][4], B[4];
double dx95, dx105, xmax, xmin, PP;
double sx, sy, sz, sw, sxx, syy, szz, sxy,
sxz, sxw, syz, syw, szw;
double SAQxy[4][4], SAQw[4], bv[4], b1,
b2, b3;
double xreg[13], yreg[13], zreg[13], wreg[13];
double reggew[4], pgew;
dx = 0.02; //
dx = Punktabstand
Lo = 10.0; //
Lo = Zentrale Last
Lp = 0.5; //
Lp = Last Planet
dt = 0.001; // dt = Zeitschritt (bei v=0.1 mit etwa 50 Schritten ein dx durchmessen)
Fo = 10.0; // Fo = Membranspannung
Mp = 1.0; // Mp = Masse des Planeten
W1 = 0.008; // W1 = Zentrale Trichtertiefe 1
W2 = 0.002; // W2 = Zentrale Trichtertiefe 2
epsilon=1e-10; // epsilon zum Abbruch der statischen Krümmungsberechnung
step0=1; // Bei Fortsetzung der Anfangskrümmung mittels backup file hier eine 1 einsetzen
step1=1; // Bei Fortsetzung der Planetenbahn mittels backup file hier eine 1 einsetzen
Nh = (N-1)/2;
x0 = (Nh/2 - Nh)*dx;// Startposition des Planeten
y0 = 0.0; // Startposition des Planeten
z0 = 0.0; // Startposition des Planeten
v0x = 0.0; // Startgeschwindigkeit des Planeten
v0y = 0.07; // Startgeschwindigkeit des Planeten
v0z = 0.0; // Startgeschwindigkeit des Planeten
sprintf(zeile,"dx=%10.5f",dx); pDC->TextOut(200, 10,zeile);
sprintf(zeile,"Lo=%10.5f",Lo); pDC->TextOut(350, 10,zeile);
sprintf(zeile,"Lp=%10.5f",Lp); pDC->TextOut(500, 10,zeile);
sprintf(zeile,"dt=%10.5f",dt); pDC->TextOut(650, 10,zeile);
sprintf(zeile,"Fo=%10.5f",Fo); pDC->TextOut(800, 10,zeile);
N2 = N*N;
N3=N2*N;
Rmax=Nh*dx*sqrt(3.0/4.0)*sqrt(8.0/9.0)*1.0001; // In z-Richtung ist der Kubus am kleinsten (Nh*ht)
Kp=5*dx; //Der Verbesserungsfaktor wird heuristisch bestimmt
N3h=N3/2;
dx95=dx*0.95;
dx105=dx*1.05;
// Das Gitter besteht aus einer tetraedrischen kubischen Struktur.
// Aus dem Würfel wird eine Kugel ausgeschnitten mit Radius Rmax.
// Jeder Massenpunkt ist durch 12 elastische Fäden mit den Nachbarn verbunden
// Das Gitter hat die Form eines Würfels
// Zu jedem Punkt die 4 räumlichen Koordinaten x,y,z,w
// Auf Matrix elems sind Zeilennummern i und Spaltennummern j in
// radialer Reihenfolge (Abstand zum Punkt x=0.0017392847351,
// y=0.001352974376, z=0,0011314159265)
// Auf Matrix conn sind in den 13 Spalten die Anzahl und Nummern der Nachbarn
// Auf Matrix reih steht die Reihenfolge der Gitterpunkte in radialer
// Reihenfolge
double *x, *y, *z, *w, *rads, *tape0, *tape1; // Vektoren rads, tape0, tap1 sind Hilfsspeicher zum Sortieren
double *ux, *uy, *uz, *uw; //Koordinaten der vom Planeten ungestörten Membran
int *ielems, *jelems, *kelems,
*conn, *reih, *itape0, *itape1;
rads
= (double *) malloc( N3h*sizeof(double));
reih
= (int *) malloc( N3*sizeof(int));
// Herstellung der tetraedrischen Struktur. Grundbaustein ist der gleichseitige
// Tetraeder mit Seitenlänge dx, Höhe der Grundfläche dh=sqrt(3/4)*dx;
// Höhe Tetraeder ht=sqrt( dh^2 - ( (1/3)dh )^2 ) = sqrt(8/9)dh
// Speicherabbildungsfunktion ist elementnummer = i*N*N + j*N + k
dh = sqrt(3.0/4.0)*dx; // Höhe Grunddreieck
ht = sqrt(8.0/9.0)*dh; // Höhe Tetraeder
dxh= dx/2.0; // Seitliche Verschiebung jeder 2. Reihe
dh3=dh/3.0; // 1/3 Höhe Grunddreieck
xmin=-Nh*dx; xmax=Nh*dx;
// rads wird als Hilfsspeicher (Radien) benutzt. reih enthält die laufende Nummer
// Es werden nur Gitterpunkte innerhalb des Kugelradius Rmax genommen
n=0;
for(i=0; i<N; i++)
{ for(j=0; j<N; j++)
{
for(k=0; k<N; k++)
{ // Jede 2. Reihe ist um dx/2 versetzt
// Zusätzlich ist jede ungerade Schicht um dx/2 verschoben
xx = ((i%2) - (k%2)) * dxh;
// Jede ungerade Schicht ist um 1/3 dh in y-Richtung versetzt gegenüber der 0-ten Schicht
yy= (k%2)*dh/3.0;
xx= xx + (j-Nh)*dx; // x ist mit j gekoppelt
yy= yy + (Nh-i)*dh; // y ist mit -i gekoppelt
zz= (k-Nh)*ht; // z ist mit k gekoppelt
radius=sqrt(
(xx-0.0017392847351)*(xx-0.0017392847351) +
(yy-0.001352974376)*(yy-0.001352974376)
+
(zz-0.0011314159265)*(zz-0.0011314159265) );
//sprintf(zeile,"radius=%10.5f",radius); pDC->TextOut(150, 90,zeile);
if(radius<=Rmax) { rads[n]=radius; reih[n]=n;
n++; }
}// for k
}// for j
} // for i
N3max=n; //Anzahl der Gitterpunkte in der Kugel bis Rmax
N3max1=N3max +1;
tape0 = (double *) malloc(
N3max1*sizeof(double));
tape1 = (double *) malloc(
N3max1*sizeof(double));
itape0= (int *) malloc(
N3max1*sizeof(int));
itape1= (int *) malloc(
N3max1*sizeof(int));
// Sortieren der Radien aufsteigend. Bubble-Sort ist zu langsam (Ordnung n*n)
// Dreibandsortierung geht schneller (Ordnung n*log2(n) )
// Da die Radien selbst nicht interessieren, sondern die Reihenfolge, wird diese
// ebenfalls umsortiert.
sortstep=0;
do{ //Trennschritt. Die Werte auf rads werden auf die beiden Bänder verteilt.
// Immer wenn sich die Sortierung auf rads ändert, wird das Band gewechselt
tape=0; j=k=0; oldrad=-1e25;
for(i=0; i<N3max; i++)
{
if(oldrad>(radsi=rads[i]))
{ if(tape==0) tape=1; else tape=0; } //Wechsle
das Band
if(tape==0) {tape0[j]=radsi; itape0[j]=reih[i]; j++; }
else
{tape1[k]=radsi; itape1[k]=reih[i]; k++; }
oldrad=radsi;
}//for i
// Zusammenführung der Bänder auf rads.
nj=j; nk=k; //Füllstand der Bänder
tape0[nj]=-1e26; tape1[nk]=-1e26; //Markiere das Ende der Bänder
i=j=k=0; sortend=1; oldrad=-1e25;
if ( (nj>0) && (nk>0)
)
{
rad0=tape0[j]; rad1=tape1[k];
do{ // Übernimm immer den kleineren Wert von Band 0 oder Band 1,
// solange die Werte nicht absteigen. Dann beginnt eine neue Teilfolge
switch( (rad0<oldrad)*2 + (rad1<oldrad) ) // verzweigt in die Werte 0, 1, 2, 3
{ case
0: { //Fall (rad0>=oldrad) &&
(rad1>=oldrad)
// Nimm den kleineren der beiden Werte
if (rad0<=rad1) {oldrad=rads[i]=rad0; reih[i]=itape0[j]; j++;
rad0=tape0[j]; }
else
{oldrad=rads[i]=rad1; reih[i]=itape1[k]; k++; rad1=tape1[k]; }
break;}
case 1: { //Fall (rad0>=oldrad) && (rad1<oldrad) Nimm rad0
oldrad=rads[i]=rad0;
reih[i]=itape0[j]; j++; rad0=tape0[j];
break;}
case 2: { //Fall (rad0<oldrad) &&
(rad1>=oldrad) Nimm rad1
oldrad=rads[i]=rad1;
reih[i]=itape1[k]; k++; rad1=tape1[k];
break;}
case 3: { // Fall (rad0<oldrad) && (rad1<oldrad) Setze oldrad zurück
// und nimm dann den kleineren Wert wie in case 0
sortend=0; oldrad = -1e25;
if
(rad0<=rad1) {oldrad=rads[i]=rad0;
reih[i]=itape0[j]; j++; rad0=tape0[j]; }
else
{oldrad=rads[i]=rad1; reih[i]=itape1[k]; k++; rad1=tape1[k]; }
break;}
}//end switch
i++;
}while( (j<nj) &&
(k<nk) );
}//if(nj
// Übernimm den Rest des Bandes, das nicht am Ende ist
if(j<nj)
{ do{ radsi=rads[i]=tape0[j]; reih[i]=itape0[j];
j++;
if(oldrad>radsi)
sortend=0;
oldrad=radsi;
i++;
}while(j<nj);
}
if(k<nk)
{ do{ radsi=rads[i]=tape1[k]; reih[i]=itape1[k];
k++;
if(oldrad>radsi)
sortend=0;
oldrad=radsi;
i++;
}while(k<nk);
}
sortstep++;
sprintf(zeile,"sortstep=%2i",sortstep); pDC->TextOut(200, 30,zeile);
}while( sortend==0 );
sprintf(zeile,"N3max=%7i",N3max); pDC->TextOut(350, 30,zeile);
// Auf Vektor reih stehen die Punktnummern, sortiert nach aufsteigenden Radien.
// Auf ielems, jelems, kelems sollen die Gitterpunktindizes in derselben sortierten Reihenfolge stehen
// Da sie aber in der Rehenfolge der Laufanweisungen for(i=0,... for(j=0;... for(k=0;... entstehen,
// wird auf itape1 die inverse Zuordnung zur Reihenfolge auf reih erzeugt: itape1[0] sagt, wo die
// Gitterpunktindizes des 0-ten Punkts stehen, itape1[1] die des 1., itape1[2] des 2. usw.
// Z.B. 3 Punkte P0, P1, P2 sind nach aufsteigenden Radien geordnet: reih[0]=P1, reih[1]=P2, reih[2]=P0
// Dann ist die inverse Zuordnung: itape1[P0]=2, itape1[P1]=0, itape1[P2]=1, d.h. Punkt P0 hatte den grössten
// Radius und damit itape1[P0]=2, Punkt P1 den kleinsten Radius Radius und damit itape1[P1]=0, und
// Punkt P2 den mittleren Radius, und damit itape1[P2]=1.
for(m=0; m<N3max; m++)
itape1[reih[m]]=m;
free(rads);
free(tape0);
free(tape1);
free(itape0);
ielems = (int *) malloc(
N3max*sizeof(int));
jelems = (int *) malloc(
N3max*sizeof(int));
kelems = (int *) malloc(
N3max*sizeof(int));
x = (double *) malloc( N3*sizeof(double));
y = (double *) malloc( N3*sizeof(double));
z = (double *) malloc( N3*sizeof(double));
w = (double *) malloc( N3*sizeof(double));
//Für die Planetenbewegung werden zur Berechnung der Flächennormalen auch die Werte der
//vom Planeten ungestörten Membran benötigt.
ux = (double *) malloc( N3*sizeof(double));
uy = (double *) malloc( N3*sizeof(double));
uz = (double *) malloc( N3*sizeof(double));
uw = (double *) malloc( N3*sizeof(double));
n=0;
W2h=W2/2;
Rmaxh=Rmax/2;
Rmax95=Rmax*0.95;
for(i=0; i<N; i++)
{ for(j=0;
j<N; j++)
{ for(k=0;
k<N; k++)
{ // Jede 2. Reihe ist um dx/2 versetzt
// Zusätzlich ist jede ungerade Schicht um dx/2 verschoben
xx = ((i%2) - (k%2)) * dxh;
// Jede ungerade Schicht ist um 1/3 dh in y-Richtung versetzt gegenüber der 0-ten Schicht
yy= (k%2)*dh/3.0;
ind=i*N2+j*N+k;
x[ind]= xx= xx + (j-Nh)*dx; // x ist mit j gekoppelt
y[ind]= yy= yy + (Nh-i)*dh; // y ist mit -i gekoppelt
z[ind]= zz= (k-Nh)*ht; // z ist mit k gekoppelt
w[ind]= 0.0;
radius=sqrt(
(xx-0.0017392847351)*(xx-0.0017392847351) +
(yy-0.001352974376)*(yy-0.001352974376) +
(zz-0.0011314159265)*(zz-0.0011314159265) );
//sprintf(zeile,"radius=%10.5f",radius); pDC->TextOut(150, 90,zeile);
if(radius<=Rmax) //Nimm nur Gitterpunkte in der Kugel
{
ind1=itape1[n]; // Wo steht dieser Gitterpunkt in der sortierten Folge
ielems[ind1]=i; jelems[ind1]=j; kelems[ind1]=k;
// Für schnellere Anfangsiteration der Krümmung Trichter anlegen
if(radius<Rmaxh) { w[ind]= -W2h -
(Rmaxh - radius)*W1/Rmax; }
else
{if(radius<Rmax95) { w[ind]= -W2h + (radius-Rmaxh)*W2/Rmax; } }
n++;
}//if(radius
}//for(k
}//for(j
}//for(i
free(itape1);
conn
= (int *) malloc( 13*N3max*sizeof(int));
// Auf Matrix reih zu jedem Gitterpunkt seine Rangnummer speichern (seine Position in elems)
// (-1 markiert unbelegte Gitterpunkte, d.h. die, die nicht in der Kugel liegen)
for(i=0; i<N3; i++) reih[i]= -1;
for(m=0; m<N3max; m++) { i=ielems[m];
j=jelems[m]; k=kelems[m]; reih[i*N2+j*N+k]=m;
}
// Auf Matrix conn sind in den 13 Spalten die Anzahl und Nummern der Nachbarn einzutragen
for(num=0; num<N3max; num++)
{ m=0; //Anzahl der Nachbarpunkte
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
// Suche in der Umgebung Nachbarpunkte im Abstand dx
for(ii=i-2; ii<(i+3); ii++)
{ for(jj=j-2;
jj<(j+3); jj++)
{ for(kk=k-2; kk<(k+3); kk++)
{ if( (ii>=0) && (ii<N) && (jj>=0)
&& (jj<N) && (kk>=0) && (kk<N) )
{ ind2=ii*N2+jj*N+kk;
dxx=x[ind]-x[ind2];
dyy=y[ind]-y[ind2];
dzz=z[ind]-z[ind2];
radius=sqrt( dxx*dxx + dyy*dyy + dzz*dzz);
if((radius>(0.99*dx)) &&
(radius<(1.01*dx)))
{
if(reih[ind2]>=0) {m++;
conn[13*num+m]=reih[ind2];}
}// if/radius
}//if ((ii>=0)
}//for kk
}// for jj
}// for ii
conn[13*num+0]=m;
if((num%100)==0) { sprintf(zeile,"num=%7i",num); pDC->TextOut(500, 30,zeile); }
}// for num
// Belastung des zentralen Punktes (Nummer 0) mit Lo, um eine Gitterkrümmung zu erhalten
// Membranspannung ist Fo
// Start eines Rechenzyklus, bei dem die Gitterpunkte abwechseln von innen nach außen
// und dann von außen nach innen angesprochen werden.
// Die statische Kraft auf einen Gitterpunkt ist die w-Komponente der Strings zu den
// Nachbarn. Die Kraft der Strings nimmt linear mit der Länge zu (ideale Elastizität)
// Ist ein Punkt nicht im statischen Kräftegleichgewicht, dann wird er in Richtung der
// resultierenden Kraft um eine Distanz dxx, dyy, dzz, dww bewegt, die proportional zur Kraft ist.
// Der Proportionalitätsfaktor Kp wird heuristisch bestimmt.
// Elastische Energie ist hier die Summe der elastischen Spannungen der einzelnen Verbindungen
// Bei Streckung dx (dx=Gitterkonstante) hat ein String die Grundspannung Fo
// und Energie E = Integral von 0 bis dx über uFo/dx du = [ u^2 (Fo/(2dx)) ] = Fodx/2
// Bei Streckung auf Länge r hat der String die Spannung Fo r/dx und Energie E = r^2 Fo / 2dx
// Die x-Komponente der Stringkraft ist Fx =r Fo/dx cos( alfa_x) = r Fo/dx (dxx/r) = dxx Fo /dx
// Um die lange Iterationszeit zu verkürzen, bis die Anfangskrümmung der Membran steht, werden zwei
// lineare Trichterfunktionen als Startauslenkungen w eingespeichert. Trichter 1 geht von Radius
// 0 bis Rmax/2 und hat eine zentrale Tiefe W1. Trichter 2 geht bis 0.95*Rmax mit Tiefe W2.
AltEe=1e25; //Alter Wert der elastischen Energie
step = 0; //Schrittzähler Anfangsiteration
step0file=0; // Auf dem Backupfile gespeicherte Schritte
pDC->TextOut(650,
30, "no backup");
char
ket[20], *pEnd;
if(step0>0)
{
//Eingabe der Membranauslenkung von einer Backup-datei
ein=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MemXyzw.txt","r"); //AT-Labor
fscanf(ein,"%6i\n",
&step0file);
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
//Warum so umständlich? Weil die Scanfunktion eine Krankheit ist!!!
fscanf(ein,"%s\n", &ket); x[ind]=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); y[ind]=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); z[ind]=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); w[ind]=strtod(ket,&pEnd);
}
fclose(ein); pDC->TextOut(650, 30, "mit backup" );
}
do{
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt in radialer Reihenfolge
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
// Nimm keine Randpunkte (solche mit weniger als n=12 Nachbarn)
n=conn[13*num+0];
Fx=Fy=Fz=Fw=0.0; if(num==0) Fw= -Lo; //Zentrale Belastung (Punkt Nummer 0)
if(n==12)
{ for(m=1;
m<=n; m++)
{ mm=conn[13*num+m]; ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm];
ind2=ii*N2+jj*N+kk;
dxx = (x[ind2] - x[ind]);
dyy = (y[ind2] -
y[ind]);
dzz = (z[ind2] - z[ind]);
dww = (w[ind2] - w[ind]);
//radius=sqrt( dxx*dxx +
dyy*dyy +dzz*dzz +dww*dww);
Fx
= Fx+dxx; Fy=Fy+dyy; Fz=Fz+dzz; Fw=Fw+dww;
}// for m
x[ind]+=Kp*Fx;
y[ind]+=Kp*Fy;
z[ind]+=Kp*Fz;
w[ind]+=Kp*Fw;
dxx=0;
}//if n==12
}//for num
// Nimm Punkte von außen nach innen
for(num=N3max-1; num>=0; num--)
{ // Nimm einen Gitterpunkt in radialer Reihenfolge
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
// Nimm keine Randpunkte (solche mit weniger als n=12 Nachbarn)
n=conn[13*num+0];
Fx=Fy=Fz=Fw=0.0; if(num==0) Fw= -Lo; //Zentrale Belastung (Punkt Nummer 0)
if(n==12)
{ for(m=1;
m<=n; m++)
{ mm=conn[13*num+m]; ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm];
ind2=ii*N2+jj*N+kk;
dxx = (x[ind2] - x[ind]);
dyy = (y[ind2] -
y[ind]);
dzz = (z[ind2] - z[ind]);
dww = (w[ind2] - w[ind]);
//radius=sqrt( dxx*dxx +
dyy*dyy +dzz*dzz +dww*dww);
Fx = Fx+dxx;
Fy=Fy+dyy; Fz=Fz+dzz; Fw=Fw+dww;
}// for m
x[ind]+=Kp*Fx;
y[ind]+=Kp*Fy;
z[ind]+=Kp*Fz;
w[ind]+=Kp*Fw;
}//if n==12
}//for num
//Berechnung der Elastischen Energie Ee
// Elastische Energie ist hier die Summe der elastischen Spannungen der einzelnen Verbindungen
// Bei Streckung dx (dx=Gitterkonstante) hat ein String die Grundspannung Fo
// und Energie E = Integral von 0 bis dx über uFo/dx du = [ u^2 (Fo/(2dx)) ] = Fodx/2
// Bei Streckung auf Länge r hat der String die Spannung Fo r/dx und Energie E = r^2 Fo / 2dx
// Die x-Komponente der Stringkraft ist Fx =r Fo/dx cos( alfa_x) = r Fo/dx (dxx/r) = dxx Fo /dx
Ee=0.0;
Fo2dx=Fo/(2*dx);
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt in radialer Reihenfolge
i=ielems[num]; j=jelems[num]; k=kelems[num];
ind=i*N2+j*N+k;
n=conn[13*num+0];
for(m=1; m<=n; m++)
{ mm=conn[13*num+m];
//Nimm nur die Verbindung zu Punkten mit größerer Nummer
if(mm>num)
{ ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm];
ind2=ii*N2+jj*N+kk;
dxx = x[ind2] - x[ind];
dyy = y[ind2] - y[ind];
dzz = z[ind2] - z[ind];
dww = w[ind2] - w[ind];
rq= dxx*dxx + dyy*dyy
+dzz*dzz + dww*dww;
Ee += rq*Fo2dx;
}//if mm>num
}// for m
}//for num
Ediff = fabs(AltEe-Ee)/Ee;
AltEe = Ee;
step++;
sprintf(zeile,"step=%5i",step); pDC->TextOut(800, 30,zeile);
sprintf(zeile,"Ee=%10.5f",Ee); pDC->TextOut(200, 50,zeile);
// Zeichnen des Gitters, aber nur die x-y-Ebene;
// y-Position 1/2 Pixelzahl nach oben und 1/3 nach rechts
// w-Position direkt nach oben bzw. unten
// x ist mit j gekoppelt, y ist mit -i gekoppelt in der Koordinatenmatrix
// Die Graphik ist unabhängig davon: x zeigt nach rechts, y schräg nach hinten
Rox=0; Roy=100;
Rux=1200; Ruy=700;
pDC->Rectangle(Rox, Roy, Rux, Ruy);
eps1 = dx*0.001;
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
// Nimm nur Punkte mit z=0
if(fabs(z[ind])<eps1)
{
// Nimm alle Nachbarn mit z=0, aber ziehe nur Linien zu Nachbarn mit größerer Nummer
n=conn[13*num];
for(m=1; m<=n; m++)
{ mm=conn[13*num+m];
if(mm>num)
{
ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm]; ind1=ii*N2+jj*N+kk;
if(fabs(z[ind1])<eps1)
{ //Bestimme Pixelkoordinaten der beiden Punkte
ix =floor(600.499 + 14*(x[ind]
+y[ind]/3) /dx);
iix=floor(600.499 +
14*(x[ind1]+y[ind1]/3)/dx);
jy =floor(400.499 - 8*
(y[ind] /2 + 20*w[ind]) /dx);
jjy=floor(400.499 - 8*
(y[ind1]/2 + 20*w[ind1]) /dx);
if((ix>=Rox)
&& (ix<=Rux) && (jy>=Roy) && (jy<=Ruy)
&&
(iix>=Rox)
&& (iix<=Rux) && (jjy>=Roy) && (jjy<=Ruy))
{ pDC->MoveTo(ix,jy);
pDC->LineTo(iix,jjy);
}
}//if(fabs(z
}//if(mm>num
}//for m=1
}//if(fabs(z
}//for num=0
//Ausgabe der Membranauslenkung auf eine Backup-datei nach jedem 19-ten Schritt
if((step%19)==0)
{
aus=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MemXyzw.txt","w"); //AT-Labor
fprintf(aus,"%6i\n",
(step+step0file) );
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
i=ielems[num]; j=jelems[num]; k=kelems[num];
ind=i*N2+j*N+k;
fprintf(aus,"%12.5f\n",x[ind]);
fprintf(aus,"%12.5f\n",y[ind]);
fprintf(aus,"%12.5f\n",z[ind]);
fprintf(aus,"%12.5f\n",w[ind]);
}
fclose(aus);
}
}while((step<=MAXSTEP1)
&& (Ediff>epsilon));
//Ausgabe der Membranauslenkung auf eine Backup-datei
aus=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MemXyzw.txt","w"); //AT-Labor
fprintf(aus,"%6i\n",
(step+step0file) );
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
i=ielems[num]; j=jelems[num]; k=kelems[num];
ind=i*N2+j*N+k;
fprintf(aus,"%12.5f\n",x[ind]);
fprintf(aus,"%12.5f\n",y[ind]);
fprintf(aus,"%12.5f\n",z[ind]);
fprintf(aus,"%12.5f\n",w[ind]);
}
fclose(aus);
//Ausgabe eines Teils der Daten auf eine Datei zur Kontrolle der Krümmung
wdx=sqrt(7.0)*dx;
aus=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MemRW.txt","w"); //AT-Labor
fprintf(aus," x y z w");
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
xx=x[ind]; yy=y[ind]; zz=z[ind]; ww=w[ind];
if( (fabs(xx)<wdx) && (
(fabs(yy)<wdx) || (fabs(zz)<wdx) ) )
fprintf(aus,"\n%12.5f%12.5f%12.5f%12.5f", xx,yy,zz,ww);
}
fclose(aus);
if(Ediff>epsilon) return; // Muss später weiteriteriert werden
// Modell:
// Schwerkrafttrichter eines Planeten mit Last Lp
// Belastung des Startpunkts des Planeten mit Lp um einen Trichter zu erhalten
// >>>>>>>> Beginn der Berechnung der Bewegung <<<<<<<<<<<<
// Im ersten Anlauf hat der Planet schwere und träge Masse, aber nicht die Membran
// Er krümmt somit lediglich die Membran.
// Seine Beschleunigung erfährt der Planet durch durch die Hangabtriebskraft
// Hangabtriebskraft: Der Planet liegt auf einer Ausgleichsfläche, die durch den
// nächstliegenden Punkt und seinen 12 Nachbarn definiert wird. Die Komponenten der negativen
// Flächennormalen ergeben die Kraftkomponenten, die auf den Planeten wirken
// Auf einen Massenpunkt der Membran wirken die statischen Zugkräfte Fx, Fy, Fz, Fw der 12 Nachbarn
// und eventuell eine Last.
// Kräfte am Planeten ergeben die Beschleunigungen ax, ay, az. Diese integriert ergeben die
// Geschwindigkeiten vx, vy, vz, diese integriert die neuen x, y, z
// Startpunkt des kleinen Trichters im Modell 2, ktp=Startpunktnummer in "elems"
// Startgeschwindigkeit v0, Bezugspunkt für die Suche der nächsten Punkte ist ibz, jbz, kbz
if(step1==0) //Gegeben sind x0, y0, z0, v0x, v0y, v0z.
{
xpl=x0; ypl=y0; zpl=z0; vxx=v0x; vyy=v0y; vzz=v0z;
pDC->TextOut(350, 50, "Neustart Orbit" );
//Öffne Datei für Ausgabe: Schritt, xpl, ypl, zpl, vxx, vyy, vzz
//aus=fopen("c:\\Dokumente und Einstellungen\\user\\MembranOrbit.txt","w"); //In der CAD-Hall und zu Hause
aus=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MembranOrbit.txt","w"); //AT-Labor
fprintf(aus," x y z\n");
}
else
{ // Eingabe der Daten xpl, ypl, zpl, vxx, vyy, vzz von der Backup-Datei
ein=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MemXV.txt","r"); //AT-Labor
//Warum so umständlich? Weil die Scanfunktion eine Krankheit ist!!!
fscanf(ein,"%s\n", &ket); xpl=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); ypl=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); zpl=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); vxx=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); vyy=strtod(ket,&pEnd);
fscanf(ein,"%s\n", &ket); vzz=strtod(ket,&pEnd);
fclose(ein); pDC->TextOut(350, 50, "Fortsetzung Orbit" );
//Öffne Datei für Ausgabe: Schritt, xpl, ypl, zpl, vxx, vyy, vzz
//aus=fopen("c:\\Dokumente und Einstellungen\\user\\MembranOrbit.txt","w"); //In der CAD-Hall und zu Hause
aus=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MembranOrbit.txt","a"); //AT-Labor
}
// Suche den nächstgelegenen Gitterpunkt ibz, jbz, kbz.
mindist=1e25;
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
xx=x[ind]-xpl; yy=y[ind]-ypl; zz=z[ind]-zpl;
dist=sqrt(xx*xx + yy*yy + zz*zz);
if(dist<mindist) {mindist=dist; ibz=i; jbz=j; kbz=k;}
}//for num
dtMp=dt/Mp;
// Konsevieren der ungestörten Membrankoordinaten
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
i=ielems[num]; j=jelems[num];
k=kelems[num]; ind=i*N2+j*N+k;
ux[ind]=x[ind]; uy[ind]=y[ind]; uz[ind]=z[ind]; uw[ind]=w[ind];
}
step=0;
// Einrichtung einer Zeichenfläche für die Bahnbewegung im x-y-Koordinatensystem
pDC->Rectangle(0,0,100,100);
do{
// Mache einen Bewegungsschritt des Planeten
xpl=xpl+vxx*dt;
ypl=ypl+vyy*dt;
zpl=zpl+vzz*dt;
//Suche in der Umgebung des Planeten die nächsten 4 Punkte, die einen Tetraeder bilden
// Indizes und Abstand auf den Vektoren iv, jv, kv und rqv speichern
// unbelegte Gitterpunkte erkennt man an reih[i*N2+j*N+k]<0
anz=0;
for(i=ibz-2; i<ibz+3; i++)
{ for(j=jbz-2;
j<jbz+3; j++)
{ for(k=kbz-2; k<kbz+3; k++)
{ ind=i*N2+j*N+k;
if(reih[ind]>=0)
{ dxx=ux[ind]-xpl;
dyy=uy[ind]-ypl;
dzz=uz[ind]-zpl;
rq = dxx*dxx + dyy*dyy +
dzz*dzz; // r-Quadrat
iv[anz]=i; jv[anz]=j; kv[anz]=k, rqv[anz]=rq; anz++;
}//end if
}}//for k//for j
}//for i
// Sortiere aufsteigend
for(m=0; m<(anz-1); m++)
{ rqmin=rqv[m];
for(i=m+1; i<anz; i++)
{ if(rqv[i]<rqmin)
{ rqmin=rqv[i]; // Das neue Minimum
rq=rqv[m]; ii=iv[m]; jj=jv[m]; kk=kv[m];
rqv[m]=rqv[i]; iv[m]=iv[i]; jv[m]=jv[i]; kv[m]=kv[i];
rqv[i]=rq; iv[i]=ii; jv[i]=jj; kv[i]=kk;
}// end
if
}//for i
}//for k
// Suche Punkt P2 so, dass er zu P0 und P1 etwa Abstand dx hat
i=iv[0]; j=jv[0]; k=kv[0];
ind=i*N2+j*N+k;
i=iv[1]; j=jv[1]; k=kv[1];
ind1=i*N2+j*N+k;
m=2; gefunden=0;
do{ i=iv[m]; j=jv[m]; k=kv[m]; ind2=i*N2+j*N+k;
dxx=ux[ind]-ux[ind2];
dyy=uy[ind]-uy[ind2];
dzz=uz[ind]-uz[ind2];
r02 = sqrt(dxx*dxx + dyy*dyy + dzz*dzz); // Abstand P0-P2
dxx=ux[ind1]-ux[ind2];
dyy=uy[ind1]-uy[ind2];
dzz=uz[ind1]-uz[ind2];
r12 = sqrt(dxx*dxx + dyy*dyy + dzz*dzz); // Abstand P1-P2
if( (r02>dx95) &&
(r02<dx105) && (r12>dx95) && (r12<dx105) )
{ //Vertausche Punkt m mit Punkt 2
rq=rqv[m]; ii=iv[m]; jj=jv[m]; kk=kv[m];
rqv[m]=rqv[2]; iv[m]=iv[2]; jv[m]=jv[2]; kv[m]=kv[2];
rqv[2]=rq; iv[2]=ii; jv[2]=jj; kv[2]=kk;
gefunden=1;
}
m++;
}while((m<anz) && (gefunden==0));
// Suche Punkt P3 so, dass er zu P0, p1 und P2 etwa Abstand dx hat
i=iv[2]; j=jv[2]; k=kv[2];
ind2=i*N2+j*N+k;
m=3; gefunden=0;
do{ i=iv[m]; j=jv[m]; k=kv[m]; ind3=i*N2+j*N+k;
dxx=ux[ind]-ux[ind3];
dyy=uy[ind]-uy[ind3];
dzz=uz[ind]-uz[ind3];
r03 = sqrt(dxx*dxx + dyy*dyy + dzz*dzz); // Abstand P0-P3
dxx=ux[ind1]-ux[ind3];
dyy=uy[ind1]-uy[ind3];
dzz=uz[ind1]-uz[ind3];
r13 = sqrt(dxx*dxx + dyy*dyy + dzz*dzz); // Abstand P1-P3
dxx=ux[ind3]-ux[ind3];
dyy=uy[ind3]-uy[ind3];
dzz=uz[ind3]-uz[ind3];
r23 = sqrt(dxx*dxx + dyy*dyy + dzz*dzz); // Abstand P2-P3
if( (r03>dx95) &&
(r03<dx105) && (r13>dx95) && (r13<dx105) &&
(r23>dx95) && (r23<dx105) )
{ //Vertausche Punkt m mit Punkt 3
rq=rqv[m]; ii=iv[m]; jj=jv[m]; kk=kv[m];
rqv[m]=rqv[3]; iv[m]=iv[3]; jv[m]=jv[3]; kv[m]=kv[3];
rqv[3]=rq; iv[3]=ii; jv[3]=jj; kv[3]=kk;
gefunden=1;
}
m++;
}while((m<anz) && (gefunden==0));
sprintf(zeile,"ii=%3i, jj=%3i",ii,jj); pDC->TextOut(500, 50,zeile);
// P0, P1, P2, P3 sind die vier nächsten Gitterpunkte am Planeten
// Sie bilden ein Tetraeder
i00=iv[0]; j00=jv[0]; k00=kv[0]; ind0=i00*N2+j00*N+k00;
i11=iv[1]; j11=jv[1]; k11=kv[1]; ind1=i11*N2+j11*N+k11;
i22=iv[2]; j22=jv[2]; k22=kv[2]; ind2=i22*N2+j22*N+k22;
i33=iv[3]; j33=jv[3]; k33=kv[3]; ind3=i33*N2+j33*N+k33;
x0 = x[ind0]; y0 = y[ind0]; z0 = z[ind0]; // Koordinaten P0
// Die 4 Belastungen der vier Eckpunkte P0, P1, P2, P3 des Tetraeders,
// in dem sich der Planet gerade befindet, ergeben sich aus der Lösung des folgenden
// Gleichungssystems:
// x0 g0 + x1 g1 + x2 g2 + x3 g3 = xpl Lp
// y0
g0 + y1 g1 + y2 g2 + y3 g3 = ypl Lp
// z0
g0 + z1 g1 + z2 g2 + z3 g3 = zpl Lp
//
g0 + g1 + g2 +
g3 = Lp
// Koeffizientenmatrix A
A[0][0] = x0; A[0][1] = x[ind1]; A[0][2]
= x[ind2]; A[0][3] = x[ind3];
A[1][0] = y0; A[1][1] = y[ind1]; A[1][2]
= y[ind2]; A[1][3] = y[ind3];
A[2][0] = z0; A[2][1] = z[ind1]; A[2][2]
= z[ind2]; A[2][3] = z[ind3];
A[3][0] = 1.0; A[3][1] = 1.0; A[3][2]
= 1.0; A[3][3] = 1.0;
//Rechte Seite B
B[0] = xpl*Lp;
B[1] = ypl*Lp;
B[2] = zpl*Lp;
B[3] = Lp;
GaussJordan(4, A, B, Lpv, &PP ); // Löst Gleichungssystem A*Lpv = B
// Der Bewegungsschritt wird gleich am Anfang der Schleife ausgeführt (siehe oben)
// Es folgt ein kleiner Zyklus, in dem die neue Raumkrümmung in unmittelbarer Umgebung
// des Planeten iteriert wird. Danach erfolgt ein Gesamtiterationsschritt mit allen Punkten
// Wir benötigen die ktp-Nummern der 4 Randpunkte P0,P1,P2,P3 des Gittertetraeders
ktpv[0]=ktp0=reih[ind0];
ktpv[1]=ktp1=reih[ind1];
ktpv[2]=ktp2=reih[ind2];
ktpv[3]=ktp3=reih[ind3];
if((step%10)==0)
{
for(step3=0; step3<MAXSTEP3; step3++)
{
for(i=ibz-3; i<ibz+4; i++)
{ for(j=jbz-3;
j<jbz+4; j++)
{ for(k=kbz-3; k<kbz+4; k++)
{
// Die Punktnummer num zu einem Punkt i,j,k finden wir ebenfalls in Matrix "reih"
ind = i*N2+j*N+k; num = reih[ind];
n=conn[13*num];
// Punktbelastungen des Zentralpunktes und der Tetraederecken
Fx=Fy=Fz=Fw=0.0;
if(num==0) Fw=
-Lo;
if(num==ktp0) Fw=
-Lpv[0];
if(num==ktp1) Fw=
-Lpv[1];
if(num==ktp2) Fw=
-Lpv[2];
if(num==ktp3) Fw= -Lpv[3];
if(n==12) // Nimm keine Randpunkte (solche mit weniger als n=12 Nachbarn)
{ for(m=1;
m<=n; m++)
{ mm=conn[13*num+m];
// Indizes des Nachbarpunktes und Matrixposition ind1
ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm]; ind1=ii*N2+jj*N+kk;
// Abstand zum Nachbarpunkt
dxx = (x[ind1] - x[ind]);
dyy = (y[ind1] - y[ind]);
dzz = (z[ind1] - z[ind]);
dww = (w[ind1] - w[ind]);
Fx +=
dxx; Fy += dyy; Fz += dzz; Fw
+= dww;
}// for m
// Iterative Änderung der Gitterpunktkoordinaten
x[ind] += Kp*Fx;
y[ind] +=
Kp*Fy;
z[ind] += Kp*Fz;
w[ind] += Kp*Fw;
}//if n==12
}//for k
}//for j
}//for i
} // for step3
//sprintf(zeile,"z0=%10.5f",z[i0*N+j0]); pDC->TextOut(650, 70,zeile);
// Gesamtiterationsschritt ähnlich wie bei der Berechnung der Anfangskrümmung
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt in radialer Reihenfolge
i=ielems[num]; j=jelems[num]; k=kelems[num];
ind=i*N2+j*N+k;
// Nimm keine Randpunkte (solche mit weniger als n=12 Nachbarn)
n=conn[13*num];
Fx=Fy=Fz=Fw=0.0; if(num==0) Fw= -Lo;
if(num==ktp0) Fw= -Lpv[0]; // Punktbelastungen
if(num==ktp1) Fw= -Lpv[1];
if(num==ktp2)
Fw= -Lpv[2];
if(num==ktp3)
Fw= -Lpv[3];
if(n==12)
{ for(m=1;
m<=n; m++)
{ mm=conn[13*num+m]; ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm]; ind1=ii*N2+jj*N+kk;
dxx = (x[ind1] -
x[ind]);
dyy = (y[ind1] -
y[ind]);
dzz = (z[ind1] - z[ind]);
dww = (w[ind1] - w[ind]);
//radius=sqrt( dxx*dxx +
dyy*dyy +dzz*dzz );
Fx = Fx+dxx;
Fy=Fy+dyy; Fz=Fz+dzz; Fw += dww;
}// for m
x[ind] += Kp*Fx;
y[ind] += Kp*Fy;
z[ind] += Kp*Fz;
w[ind] += Kp*Fw;
}//if n==12
}//for num
// Nimm Punkte von außen nach innen
for(num=N3max-1; num>=0; num--)
{ // Nimm einen Gitterpunkt in radialer Reihenfolge
i=ielems[num]; j=jelems[num]; k=kelems[num]; ind=i*N2+j*N+k;
// Nimm keine Randpunkte (solche mit weniger als n=12 Nachbarn)
n=conn[13*num];
Fx=Fy=Fz=Fw=0.0; if(num==0) Fw= -Lo;
if(num==ktp0) Fw= -Lpv[0]; // Punktbelastungen
if(num==ktp1) Fw= -Lpv[1];
if(num==ktp2)
Fw= -Lpv[2];
if(num==ktp3)
Fw= -Lpv[3];
if(n==12)
{ for(m=1;
m<=n; m++)
{ mm=conn[13*num+m]; ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm]; ind1=ii*N2+jj*N+kk;
dxx = (x[ind1] -
x[ind]);
dyy = (y[ind1] -
y[ind]);
dzz = (z[ind1] - z[ind]);
dww = (w[ind1] - w[ind]);
//radius=sqrt( dxx*dxx +
dyy*dyy +dzz*dzz );
Fx = Fx+dxx;
Fy=Fy+dyy; Fz=Fz+dzz; Fw += dww;
}// for m
x[ind] += Kp*Fx;
y[ind] += Kp*Fy;
z[ind] += Kp*Fz;
w[ind] += Kp*Fw;
}//if n==12
}//for num
// Zeichnen des Gitters, aber nur die gebogene x-y-Ebene auf der z-Höhe des Planeten;
// y-Position 1/2 Pixelzahl nach oben und 1/3 nach rechts
// w-Position direkt nach oben bzw. unten
// x ist mit j gekoppelt, y ist mit -i gekoppelt in der Koordinatenmatrix
// Die Graphik ist unabhängig davon: x zeigt nach rechts, y schräg nach hinten
Rox=0; Roy=100;
Rux=1200; Ruy=700;
pDC->Rectangle(Rox, Roy, Rux, Ruy);
eps1 = dx*0.001;
for(num=0; num<N3max; num++)
{ // Nimm einen Gitterpunkt
k=kelems[num];
// Nimm nur Punkte der z-Ebene, die dem Planeten am nächsten liegt
if(k==k00)
{ i=ielems[num]; j=jelems[num]; ind=i*N2+j*N+k;
// Nimm alle Nachbarn auf derselben z-Ebene, aber ziehe nur Linien zu Nachbarn mit größerer Nummer
n=conn[13*num];
for(m=1; m<=n; m++)
{ mm=conn[13*num+m];
if(mm>num)
{
kk=kelems[mm];
if(kk==k00)
{ ii=ielems[mm]; jj=jelems[mm]; ind1=ii*N2+jj*N+kk;
//Bestimme Pixelkoordinaten der beiden Punkte
ix =floor(600.499 + 14*(x[ind]
+y[ind]/3) /dx);
iix=floor(600.499 +
14*(x[ind1]+y[ind1]/3)/dx);
jy =floor(400.499 - 8*
(y[ind] /2 + 20*w[ind]) /dx);
jjy=floor(400.499 - 8*
(y[ind1]/2 + 20*w[ind1]) /dx);
if((ix>=Rox)
&& (ix<=Rux) && (jy>=Roy) && (jy<=Ruy)
&&
(iix>=Rox)
&& (iix<=Rux) && (jjy>=Roy) && (jjy<=Ruy))
{ pDC->MoveTo(ix,jy);
pDC->LineTo(iix,jjy);
}
}//if(fabs(z
}//if(mm>num
}//for m=1
}//if(fabs(z
}//for num=0
//Ausgabe der Bahnkurve in die kleine Zeichenfläche
ix =floor((xpl-xmin)/(xmax-xmin)*100);
jy =100-floor((ypl-xmin)/(xmax-xmin)*100);
pDC->SetPixelV(ix,jy, RGB(0,0,0));
//Schreibe Datei für Ausgabe: Schritt, xpl, ypl, zpl
fprintf(aus,"%12.5f%12.5f%12.5f\n",xpl, ypl, zpl);
// Ausgabe von Ort und Geschwindigkeit für eine spätere Fortsetzung der Rechnung
aus1=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MemXV.txt","w"); //AT-Labor
fprintf(aus1,"%12.5f\n%12.5f\n%12.5f\n%12.5f\n%12.5f\n%12.5f\n", xpl, ypl, zpl, vxx, vyy, vzz);
fclose(aus1);
} //if((step%10
// Lege durch die 4 zentralen Punkte und ihre jeweils 12 Nachbarpunkte (n=12) eine
// gewichtete Ausgleichsebene. Gewichte sind die Punktbelastung, aber ihre Summe auf 1 normiert.
// Ebenengleichung: w = b0 + b1*x + b2*y + b3*z
// Die Koeffizienten b1, b2, b3 bzw. deren negative Werte geben die Raumnormale ( -b1 , -b2, -b3, 1 )
// Koeffizientenmatrix:
// (Summe x^2) - ((Summe x)^2)/n (Summe x*y) - (Summe x)(Summe y)/n (Summe x*z) - (Summe x)(Summe z)/n
//
(Summe x*y) - (Summe x)(Summe y)/n
(Summe y^2) - ((Summe y)^2)/n (Summe y*z) - (Summe y)(Summe z)/n
//
(Summe x*z) - (Summe x)(Summe z)/n
(Summe y*z) - (Summe y)(Summe z)/n
(Summe z^2) - ((Summe z)^2)/n
// Rechte Seite
//
(Summe x*w) - (Summe x)(Summe w)/n
(Summe y*w) - (Summe y)(Summe w)/n
(Summe z*w) - (Summe z)(Summe w)/n
// Normieren der 4 Gewichte auf Vektor reggew[4]
for(pnr=0; pnr<4; pnr++) reggew[pnr]=Lpv[pnr]/Lp; // Summe der Lpv ist immer Lp
// Löschen der SAQxy-Matrix und des SAQw-Vektors
for(ii=0; ii<3; ii++) { SAQw[ii]=0.0; for(jj=0; jj<3; jj++) SAQxy[ii][jj]=0.0;}
for(pnr=0; pnr<4; pnr++)
{
// Löschen der Summen zu einem Punkt
sx =sy =sz =sw =sxx =syy =szz =sxy =sxz =sxw =syz =syw =szw = 0.0;
// Punktgewicht ist pgew
pgew = reggew[pnr];
if(fabs(pgew)>0.0001)
{
//Speichere die Indizes des aktuellen Punktes pnr nach ii00, jj00, kk00,
ii00=iv[pnr]; jj00=jv[pnr]; kk00=kv[pnr];
// Zusammensuchen der 12 Koordinatentupel
ind=ii00*N2+jj00*N+kk00; xreg[0]=ux[ind]; yreg[0]=uy[ind]; zreg[0]=uz[ind]; wreg[0]=uw[ind];
// Nimm alle 12 Nachbarn des aktuellen Gitterpunktes zum Planeten. Dieser Punkt hat Nummer ktpv[pnr]
for(m=1; m<=12; m++)
{ mm=conn[13*ktpv[pnr]+m];
ii=ielems[mm]; jj=jelems[mm]; kk=kelems[mm]; ind=ii*N2+jj*N+kk;
xreg[m]=ux[ind]; yreg[m]=uy[ind]; zreg[m]=uz[ind]; wreg[m]=uw[ind];
}//form=1
// Berechnung der Summen
for(i=0; i<13; i++)
{ xx=xreg[i]; yy=yreg[i]; zz=zreg[i]; ww=wreg[i]; //Die Werte eines Punktes
sx += xx; sy += yy; sz += zz; sw +=
ww; //Summen
sxx += xx*xx; syy += yy*yy; szz += zz*zz; //Quadratsummen
sxy += xx*yy; sxz += xx*zz; sxw += xx*ww; //Produktsummen
syz += yy*zz; syw += yy*ww;
szw += zz*ww;
}//for(i=0
// Obere Hälfte der Koeffizientenmatrix SAQxy und Rechte Seite SAQw
SAQxy[0][0] += (sxx - sx*sx/13)*pgew;
SAQxy[1][1] += (syy - sy*sy/13)*pgew;
SAQxy[2][2] += (szz - sz*sz/13)*pgew;
SAQxy[0][1] += (sxy - sx*sy/13)*pgew;
SAQxy[0][2] += (sxz - sx*sz/13)*pgew;
SAQxy[1][2] += (syz - sy*sz/13)*pgew;
SAQw[0] += (sxw-sx*sw/13)*pgew;
SAQw[1] += (syw-sy*sw/13)*pgew;
SAQw[2] += (szw-sz*sw/13)*pgew;
}//if(fabs(pgew)
}//for(pnr=0
SAQxy[1][0]= SAQxy[0][1]; // Wegen Symmetrie der Matrix
SAQxy[2][0]= SAQxy[0][2];
SAQxy[2][1]= SAQxy[1][2];
// Lösung des kleinen Gleichungssystems SAQxy * bv = SAQw
GaussJordan(3, SAQxy, SAQw, bv, &PP ); // Löst Gleichungssystem
// Die Regressionskoeffizienten und damit die Komponenten der Flächennormalen sind
b1=bv[0]; b2=bv[1]; b3=bv[2];
// Die seitlich auf den Planeten wirkenden Kräfte sind -Lp*bi
Fx = -Lp*b1; Fy = -Lp*b2; Fz =
-Lp*b3;
// An dieser Stelle liegen die Kräfte Fx, Fy und Fz vor, die den Planeten seitlich beschleunigen
// Beschleunigung ax=Fx/Mp usw.. Die Integration erfolgt mit Euler. v_neu = v_alt + a * dt .
vxx = vxx + Fx*dtMp; vyy = vyy + Fy*dtMp; vzz = vzz + Fz*dtMp;
ibz=i00; jbz=j00; kbz=k00; // Im nächsten Schritt Zentralpunkt der Suche
sprintf(zeile,"x=%8.5f",xpl); pDC->TextOut(200, 70,zeile);
sprintf(zeile,"y=%8.5f",ypl); pDC->TextOut(350, 70,zeile);
sprintf(zeile,"z=%8.5f",zpl); pDC->TextOut(500, 70,zeile);
sprintf(zeile,"vx=%8.5f",vxx); pDC->TextOut(650, 70,zeile);
sprintf(zeile,"vy=%8.5f",vyy); pDC->TextOut(800, 70,zeile);
sprintf(zeile,"vz=%8.5f",vzz); pDC->TextOut(950, 70,zeile);
step++; sprintf(zeile,"step=%i", step); pDC->TextOut(650,50,zeile);
}while(step<=MAXSTEP2);
fclose(aus);
// Ausgabe von Ort und Geschwindigkeit für eine spätere Fortsetzung der Rechnung
aus1=fopen("c:\\Dokumente und Einstellungen\\AT Vorlesung\\Eigene Dateien\\MemXV.txt","w"); //AT-Labor
fprintf(aus1,"%12.5f\n%12.5f\n%12.5f\n%12.5f\n%12.5f\n%12.5f\n", xpl, ypl, zpl, vxx, vyy, vzz);
fclose(aus1);
free(x);
free(y);
free(z);
free(w);
free(ux);
free(uy);
free(uz);
free(uw);
free(ielems); free(jelems); free(kelems);
free(conn);
free(reih);
pDC->TextOut(200,70,"EndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEndeEnde");
}
/////////////////////////////////////////////////////////////////////////////
//
CMembrane4DView Drucken
BOOL
CMembrane4DView::OnPreparePrinting(CPrintInfo* pInfo)
{
// Standardvorbereitung
return DoPreparePrinting(pInfo);
}
void
CMembrane4DView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
// ZU ERLEDIGEN: Zusätzliche Initialisierung vor dem Drucken hier einfügen
}
void
CMembrane4DView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)
{
// ZU ERLEDIGEN: Hier Bereinigungsarbeiten nach dem Drucken einfügen
}
/////////////////////////////////////////////////////////////////////////////
//
CMembrane4DView Diagnose
#ifdef
_DEBUG
void
CMembrane4DView::AssertValid() const
{
CView::AssertValid();
}
void
CMembrane4DView::Dump(CDumpContext& dc) const
{
CView::Dump(dc);
}
CMembrane4DDoc* CMembrane4DView::GetDocument() // Die endgültige (nicht zur Fehlersuche kompilierte) Version ist Inline
{
ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CMembrane4DDoc)));
return (CMembrane4DDoc*)m_pDocument;
}
#endif
//_DEBUG
============================= *** =============================