summaryrefslogtreecommitdiff
path: root/Physikunit.pas
diff options
context:
space:
mode:
Diffstat (limited to 'Physikunit.pas')
-rw-r--r--Physikunit.pas228
1 files changed, 145 insertions, 83 deletions
diff --git a/Physikunit.pas b/Physikunit.pas
index 4861a14..5d8881a 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -27,6 +27,7 @@ type
eqRho,eqJX,eqJY
);
tMaterieFeldGroesze = (mfPY);
+ tMaterieSpeicherGroesze = (msPX,msPXSqr,msPY,msVX,msVY,msN,msPXRipple);
tTeilchenSpeziesGroeszen = (tsgMasse,tsgIQdrMasse,tsgLadung);
tEMQuellen = array[tEMQuellGroesze] of extended;
@@ -46,9 +47,11 @@ type
tAusgabeDatei = class
private
ableitung: boolean;
- emFelder: tEMFeldGroesze;
- emQuellen: tEMQuellGroesze;
- teilchen, // -2: EM-Feld; -1: Gesamt-EM-Quelle; sonst: EM-Quelle der entsprechenden Teilchen
+ emFeld: tEMFeldGroesze;
+ emQuelle: tEMQuellGroesze;
+ matFeld: tMaterieSpeicherGroesze;
+ wasSpeichern, // 0: EM-Feld; 1: EM-Quelle; 2: Mat-Feld
+ teilchen, // -1: Gesamt-Feld; sonst: Feld des entsprechenden Teilchen
nNum,tCnt: longint;
pre,suf: string;
datei: file;
@@ -132,7 +135,7 @@ type
{$UNDEF LiKotRaumPunktHeader}
function nichtnegativieren: extended; inline;
function impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; overload; inline;
- function impulsIntegral(teilchen: longint): extended; overload; inline;
+ function impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended; overload; inline;
procedure initialisiereDichte(teilchen: longint; breite,n: extended); inline;
procedure reflektiereTeilchen(rechts: boolean); inline;
procedure setzeNull; inline;
@@ -211,6 +214,9 @@ const
emQuellNamen: array[tEMQuellGroesze] of string = (
'RHO','JX','JY'
);
+ matSpeicherNamen: array[tMaterieSpeicherGroesze] of string = (
+ 'PX','PXSQR','PY','VX','VY','N','PXRIPPLE'
+ );
minAusgabeBuffer = 1024*1024;
var
@@ -222,6 +228,7 @@ constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokol
var
emF: tEMFeldGroesze;
emQ: tEMQuellGroesze;
+ maF: tMaterieSpeicherGroesze;
num: longint;
abl,gef: boolean;
begin
@@ -238,8 +245,10 @@ begin
else
teilchen:=-1;
- emFelder:=efAX;
- emQuellen:=eqRho;
+ emFeld:=efAX;
+ emQuelle:=eqRho;
+ matFeld:=msN;
+ wasSpeichern:=0;
feldName:=ansiUpperCase(feldName);
nNum:=0; // Header 0 wurde also noch nicht geschrieben
zeitAnz:=-1;
@@ -250,7 +259,8 @@ begin
if not gef then
for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin
- emFelder:=emF;
+ emFeld:=emF;
+ wasSpeichern:=0;
ableitung:=abl;
teilchen:=-2;
gef:=true;
@@ -261,24 +271,45 @@ begin
for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
if emQuellNamen[emQ] = feldName then begin
ableitung:=false;
- emQuellen:=emQ;
+ emQuelle:=emQ;
+ wasSpeichern:=1;
+ gef:=true;
+ break;
+ end;
+ if not gef then
+ for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do
+ if matSpeicherNamen[maF] = feldName then begin
+ ableitung:=false;
+ matFeld:=maF;
+ wasSpeichern:=2;
gef:=true;
break;
end;
if not gef then begin
pro.schreibe('tAusgabeDatei.create kennt Größe '''+feldName+''' nicht!',true);
+ nam:='';
+ for abl:=false to true do
+ for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
+ nam:=nam+', '+copy('D',1,byte(abl))+emFeldNamen[emF];
+ for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
+ nam:=nam+', '+emQuellNamen[emQ];
+ for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do
+ nam:=nam+', '+matSpeicherNamen[maF];
+ delete(nam,1,2);
+ pro.schreibe('Ich kenne nur: '+nam,true);
+
pro.destroyAll;
halt(1);
end;
- if teilchen=-2 then
- pre:=emFeldNamen[emF]
- else begin
- pre:=emQuellNamen[emQ];
- if teilchen>=0 then
- pre:=pre+inttostr(teilchen+1);
+ case wasSpeichern of
+ 0: pre:=emFeldNamen[emFeld];
+ 1: pre:=emQuellNamen[emQuelle];
+ 2: pre:=matSpeicherNamen[matFeld];
end;
+ if teilchen>=0 then
+ pre:=pre+inttostr(teilchen+1);
if ableitung then
pre:='d'+pre;
@@ -320,13 +351,13 @@ end;
function tAusgabeDatei.dump: string;
begin
- if teilchen=-2 then
- result:=emFeldNamen[emFelder]
- else begin
- result:=emQuellNamen[emQuellen];
- if teilchen>=0 then
- result:=result+inttostr(teilchen);
+ case wasSpeichern of
+ 0: result:=emFeldNamen[emFeld];
+ 1: result:=emQuellNamen[emQuelle];
+ 2: result:=matSpeicherNamen[matFeld];
end;
+ if teilchen>=0 then
+ result:=result+'['+inttostr(teilchen+1)+']';
if ableitung then
result:=result+'''';
result:=result+' -> '+pre+'-$i'+suf;
@@ -385,40 +416,40 @@ begin
inc(bufPos,sizeof(integer));
sX:=cX-Min(gitter.dX,sDX)/2;
- if teilchen=-2 then begin
- for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
- while cX>=sX do begin
- dec(cnt);
- move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFelder,ableitung],(buf+bufPos)^,sizeof(extended));
- inc(bufPos,sizeof(extended));
- sX:=sX+sDX;
+ case wasSpeichern of
+ 0: // em-Feld speichern
+ for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
+ while cX>=sX do begin
+ dec(cnt);
+ move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFeld,ableitung],(buf+bufPos)^,sizeof(extended));
+ inc(bufPos,sizeof(extended));
+ sX:=sX+sDX;
+ end;
+ cX:=cX+gitter.dX;
end;
- cX:=cX+gitter.dX;
- end;
- end
- else if teilchen=-1 then begin
- for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
- while cX>=sX do begin
- dec(cnt);
- move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emQuellen[emQuellen],(buf+bufPos)^,sizeof(extended));
- inc(bufPos,sizeof(extended));
- sX:=sX+sDX;
+ 1: // em-Quelle
+ for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
+ while cX>=sX do begin
+ dec(cnt);
+ val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuelle);
+ move(val,(buf+bufPos)^,sizeof(extended));
+ inc(bufPos,sizeof(extended));
+ sX:=sX+sDX;
+ end;
+ cX:=cX+gitter.dX;
end;
- cX:=cX+gitter.dX;
- end;
- end
- else begin
- for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
- while cX>=sX do begin
- dec(cnt);
- val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuellen);
- move(val,(buf+bufPos)^,sizeof(extended));
- inc(bufPos,sizeof(extended));
- sX:=sX+sDX;
+ 2: // Materiefeld
+ for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
+ while cX>=sX do begin
+ dec(cnt);
+ val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,matFeld);
+ move(val,(buf+bufPos)^,sizeof(extended));
+ inc(bufPos,sizeof(extended));
+ sX:=sX+sDX;
+ end;
+ cX:=cX+gitter.dX;
end;
- cX:=cX+gitter.dX;
- end;
- end;
+ end{of case};
if cnt<>0 then begin
pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true);
pro.destroyall;
@@ -590,17 +621,18 @@ end;
function tImpulsPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen
var
- i,richtung,entfernung: longint;
defizit: extended;
+ i: longint;
pro: tProtokollant;
- jemandErreicht,vorZurueck: boolean;
+(* i,richtung,entfernung: longint;
+ jemandErreicht,vorZurueck: boolean; *)
begin
result:=0;
- for i:=0 to length(werte)-1 do
+ for i:=0 to length(werte)-1 do begin
if werte[i,false]<0 then begin
defizit:=-werte[i,false];
result:=result+defizit;
- werte[i,false]:=0;
+ werte[i,false]:=0; (*
entfernung:=1;
jemandErreicht:=true;
while (defizit>0) and jemandErreicht do begin
@@ -616,8 +648,15 @@ begin
pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'WertePunkt.nichtnegativieren');
pro.schreibe('Ich konnte ein Teilchendefizit der Sorte '+inttostr(i+1)+' nicht auffüllen!',true);
raumpunkt.felder.gitter.abbrechen;
- end;
+ end; *)
end;
+ if werte[i,false]>1E10 then begin
+ pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'impulsPunkt.nichtnegativieren');
+ pro.schreibe('An einer Stelle im Phasenraum sind bereits mehr als 10^100 Teilchen, ich breche ab!',true);
+ pro.free;
+ raumpunkt.felder.gitter.abbrechen;
+ end;
+ end;
end;
procedure tImpulsPunkt.akkumuliereEMQuellen(var emQuellen: tEMQuellen; pX,dVP: extended);
@@ -688,16 +727,16 @@ begin
+ raumpunkt.matFelder[i,mfPY,false] * iMGamma[i] * raumpunkt.emFelder[efBZ,false]);
werte[i,true]:=
tX * (
- ( nachbarn[0,true].werte[i,false] - nachbarn[0,false].werte[i,false])
- - sign(tX)*raumpunkt.felder.teilchen[i].diffusion[0] * (
+ (nachbarn[0,true].werte[i,false] - nachbarn[0,false].werte[i,false])
+ + sign(tX)*raumpunkt.felder.teilchen[i].diffusion[0] * (
nachbarn[0,true].werte[i,false] +
nachbarn[0,false].werte[i,false] -
2*werte[i,false]
)
) * iDX/2
+ tP * (
- ( nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false])
- - sign(tP)*raumpunkt.felder.teilchen[i].diffusion[1] * (
+ (nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false])
+ + sign(tP)*raumpunkt.felder.teilchen[i].diffusion[1] * (
nachbarn[1,true].werte[i,false] +
nachbarn[1,false].werte[i,false] -
2*werte[i,false]
@@ -846,8 +885,6 @@ begin
end;
function tRaumPunkt.impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended;
-var
- i: longint;
begin
if teilchen<0 then begin
result:=emQuellen[emQ]; // das ist leicht :-)
@@ -857,12 +894,50 @@ begin
result:=0;
case emQ of
- eqRho:
+ eqRho: result:=impulsIntegral(teilchen,msN);
+ eqJX: result:=impulsIntegral(teilchen,msVX);
+ eqJY: result:=impulsIntegral(teilchen,msVY);
+ end{of case};
+ result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung];
+end;
+
+function tRaumPunkt.impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended;
+var
+ i: longint;
+begin
+ if teilchen<0 then begin
+ result:=0;
+ for i:=0 to length(felder.teilchen)-1 do
+ result:=result+impulsIntegral(i,maF);
+ exit;
+ end;
+
+ result:=0;
+
+ case maF of
+ msN:
for i:=0 to aP-1 do
result:=
result +
phasenraum[i].werte[teilchen,false]*dP;
- eqJX: begin
+ msPX:
+ for i:=0 to aP-1 do
+ result:=
+ result +
+ phasenraum[i].werte[teilchen,false]*dP * (i-aP/2)*dP;
+ msPXSqr:
+ for i:=0 to aP-1 do
+ result:=
+ result +
+ phasenraum[i].werte[teilchen,false]*dP * sqr((i-aP/2)*dP);
+ msPY: begin
+ for i:=0 to aP-1 do
+ result:=
+ result +
+ phasenraum[i].werte[teilchen,false]*dP;
+ result:=result * matFelder[teilchen,mfPY,false];
+ end;
+ msVX: begin
for i:=0 to aP-1 do
result:=
result +
@@ -870,7 +945,7 @@ begin
/ sqrt(1 + (sqr((i-aP/2)*dP) + sqr(matFelder[teilchen,mfPY,false]))*felder.teilchen[teilchen].eigenschaften[tsgIQdrMasse]);
result:=result / felder.teilchen[teilchen].eigenschaften[tsgMasse];
end;
- eqJY: begin
+ msVY: begin
for i:=0 to aP-1 do
result:=
result +
@@ -878,24 +953,11 @@ begin
/ sqrt(1 + (sqr((i-aP/2)*dP) + sqr(matFelder[teilchen,mfPY,false]))*felder.teilchen[teilchen].eigenschaften[tsgIQdrMasse]);
result:=result * matFelder[teilchen,mfPY,false] / felder.teilchen[teilchen].eigenschaften[tsgMasse];
end;
+ msPXRipple:
+ for i:=1 to aP-1 do
+ result:=
+ result + abs(phasenraum[i].werte[teilchen,false]-phasenraum[i-1].werte[teilchen,false])*dP;
end{of case};
- result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung];
-end;
-
-function tRaumPunkt.impulsIntegral(teilchen: longint): extended;
-var
- i,j: longint;
-begin
- result:=0;
-
- if teilchen<0 then begin
- for j:=0 to length(matFelder)-1 do
- for i:=0 to aP-1 do
- result:=result + phasenraum[i].werte[j,false]*dP;
- end
- else
- for i:=0 to aP-1 do
- result:=result + phasenraum[i].werte[teilchen,false]*dP;
end;
procedure tRaumPunkt.initialisiereDichte(teilchen: longint; breite,n: extended);
@@ -1086,7 +1148,7 @@ begin
for i:=0 to length(massen)-1 do begin
dens:=0;
for j:=0 to length(inhalt)-1 do
- dens:=dens+inhalt[j].impulsIntegral(i);
+ dens:=dens+inhalt[j].impulsIntegral(i,msN);
dens:=dens*teilchen[i].eigenschaften[tsgMasse]*gitter.dX;
prot.schreibe('n['+inttostr(i+1)+'] = '+floattostr(dens)+' (relative Abweichung: '+floattostr(dens/massen[i]-1)+')',true);
end;
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-11-21 22:52:45 +0000