From a6d7fdc578d83624d351ee0c07391472d7142676 Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Mon, 27 Jul 2015 16:18:15 +0200 Subject: Tagesendstand, compöilierfaehig aber nicht lauffaehig MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Physikunit.pas | 1016 ++++++++++++++++++++++++++++++------------------- Plasmapropagation | Bin 3344888 -> 3387904 bytes Plasmapropagation.lpr | 73 +--- Plasmapropagation.lps | 186 ++++----- error | 521 ------------------------- input.plap | 16 +- 6 files changed, 750 insertions(+), 1062 deletions(-) delete mode 100644 error diff --git a/Physikunit.pas b/Physikunit.pas index 71e6a79..492866c 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -2,75 +2,99 @@ unit Physikunit; {$mode objfpc}{$H+} -{ $DEFINE Zeitschrittueberwachung} -{ $DEFINE Dichteueberwachung} +{$DEFINE Zeitschrittueberwachung} +{$DEFINE Dichteueberwachung} interface uses - Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit; + Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, crt; type - TZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier); - TVerteilungsfunktion = function(x: extended): extended; - TFeldInhalt = ( - fiA,fiAX,fiAY,fiAZ, - fidAXDT,fidAYDT,fidAZDT, - fiN,fidPhiDX,fidPsiDX, - fiP,fiPX,fiPY,fiPZ, - fiGamma,fiiGamma); - - TAusgabeDatei = record + tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier); + tVerteilungsfunktion = function(x: extended): extended; + tEMFeldInhalt = ( + efA,efAX,efAY,efAZ, + efDAXDT,efDAYDT,efDAZDT, + efDPhiDX + ); + tMaterieFeldInhalt = ( + mfN,mfDPsiDX, + mfP,mfPX,mfPY,mfPZ, + mfGamma,mfIGamma + ); + + { tAusgabeDatei } + + tAusgabeDatei = class ableitung: boolean; - inhalt: tFeldInhalt; + emInhalt: tEMFeldInhalt; + matInhalt: tMaterieFeldInhalt; + teilchen: longint; // wenn < 0, dann EM-Feld name: string; datei: file; + constructor create(feldName,prefix,suffix: string; prot: tProtokollant); + destructor destroy; override; end; - tFelder = class; + tGitter = class; + + { tWertePunkt } - { tWerteGitter } + tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen + matWerte: array of array[tMaterieFeldInhalt,boolean] of extended; + // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen + emWerte: array[tEMFeldInhalt,boolean] of extended; + // EM-Felder und deren Zeitableitungen + // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! + lN,rN: tWertePunkt; - tWerteGitter = class(tObject) // repräsentiert ein Gitter von Werten und deren Zeitableitungen - werte: array[boolean] of array of extended; - par: tFelder; + constructor create(linkerNachbar: tWertePunkt; materien: longint); + destructor destroy; override; + procedure berechneGammaUndP; + procedure berechneNAbleitung(iDX,pDNMax: extended); + procedure berechneAbleitungen(dT,dX,iDX: extended); + procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt + procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt end; { TFelder } tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen + private + matAnz: longint; + lichters: array[boolean] of tMyStringlist; + procedure setzeRaender(iDX: extended); public - materieFelder: array of tWerteGitter; // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen - elektroMagnetischeFelder: tWerteGitter; // EM-Felder und deren Zeitableitungen - // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! - dX,iDX,x,pDNMax: extended; + inhalt: array of tWertePunkt; + par: tGitter; - constructor create(deltaX,pDNMa: extended); - procedure berechneGammaUndP; - procedure berechneNAbleitung; - procedure berechneAbleitungen(dT: extended); + constructor create(groesse: longint; dichten,lichter: tMyStringList; parent: tGitter); + destructor destroy; override; + procedure berechneAbleitungen(dT,dX,iDX,pDNMax: extended); + procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt + procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt end; { tGitter } tGitter = class(tObject) private - groesse: longint; - prot: tProtokollant; - procedure berechneAbleitungen(Felder: longint; dt: extended); overload; - procedure berechneAbleitungen(Felder: longint; dt: extended; out dTMax: extended); overload; - procedure setzeRaender(Felder: longint); + prot: tProtokollant; + zeitverfahren: tZeitverfahren; + pDNMaxDynamisch: boolean; + dX,iDX,pDNMax,dTMaximum,xl,t: extended; + kvs: tKnownValues; + procedure diffusionsTermAnpassen(pro: tProtokollant); public + aktuelleFelder: longint; felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern - dTMaximum,dX,iDX,xl,xr,t,pDNMax: extended; - pDNMaxDynamisch: boolean; - Zeitverfahren: TZeitverfahren; - Al,Ar: TVerteilungsfunktion; - constructor create(size: longint; deltaT,deltaX,pDNMa: extended; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); + + constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); destructor destroy; override; - procedure iteriereSchritt(var dT: extended); - procedure macheAusgabe(AD: TAusgabedatei; sDX: extended); + procedure iteriereSchritt(dT: extended); + procedure macheAusgabe(aD: tAusgabedatei; sDX: extended); function gibErhaltungsgroessen: string; end; @@ -78,59 +102,157 @@ type tSimulation = class(tObject) private - prot: tProtokollant; - gitter: tGitter; - kvs: tKnownValues; + prot: tProtokollant; + gitter: tGitter; + dT,tEnde,sT,sDT,sDX: extended; + fortschrittsAnzeige: boolean; + ausgabeDateien: array of tAusgabeDatei; public constructor create(inName: string; Protokollant: tProtokollant); destructor destroy; override; + function iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? end; +implementation + const - OptionNamen: array[TFeldInhalt] of string = ( - 'A','AX','AY','AZ','dAXDT','dAYDT','dAZDT','N','dPhiDX','dPsiDX','P','PX','PY','PZ','Gamma','iGamma'); - AusgabeNamen: array[TFeldInhalt] of string = ( - 'fiA','fiAX','fiAY','fiAZ','fidAXDT','fidAYDT','fidAZDT','fiN','fidPhiDX','fidPsiDX','fiP','fiPX','fiPY','fiPZ','fiGamma','fiiGamma'); - AusgabeHilfe = '--d?A[XYZ]? --d?dA[XYZ]dt --d?N --d?dPhiDX --d?dPsiDX --P[XYZ]? --Gamma --iGamma'; + emFeldNamen: array[tEMFeldInhalt] of string = ( + 'A','AX','AY','AZ','DAXDT','DAYDT','DAZDT','DPHIDX' + ); + matFeldNamen: array[tMaterieFeldInhalt] of string = ( + 'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA' + ); -implementation +{ tAusgabeDatei } + +constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant); +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + num: longint; + abl,gef: boolean; +begin + inherited create; + num:=length(feldName); + while (num>0) and (feldName[num] in ['0'..'9']) do + dec(num); + inc(num); + if num<=length(feldName) then begin + teilchen:=strtoint(copy(feldName,num,length(feldName))); + delete(feldName,num,length(feldName)); + end + else + teilchen:=0; + + emInhalt:=efA; + matInhalt:=mfN; + feldName:=ansiUpperCase(feldName); + + gef:=false; + for abl:=false to true do begin + if not gef then + for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do + if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin + emInhalt:=emF; + teilchen:=-1; + gef:=true; + break; + end; + if not gef then + for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do + if copy('D',1,byte(abl))+matFeldNamen[maF] = feldName then begin + matInhalt:=maF; + gef:=true; + break; + end; + end; + + if not gef then begin + prot.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true); + halt(1); + end; + + if teilchen>=0 then + name:=matFeldNamen[maF]+inttostr(teilchen) + else + name:=emFeldNamen[emF]; -{ TFeld } + if abl then + name:='d'+name; + name:=prefix+name+suffix; -constructor TFeld.create(deltaX,pDNMa: extended); -var FI: TFeldinhalt; - Abl: Boolean; + assignFile(datei,name); + rewrite(datei); +end; + +destructor tAusgabeDatei.destroy; +begin + closeFile(datei); + inherited destroy; +end; + +{ tWertePunkt } + +constructor tWertePunkt.create(linkerNachbar: tWertePunkt; materien: longint); +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + abl: boolean; + i: longint; begin inherited create; - lN:=nil; + fillchar(matWerte,sizeof(matWerte),#0); + setlength(matWerte,materien); + for i:=0 to length(matWerte)-1 do begin + for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do + for abl:=false to true do + matWerte[i,maF,abl]:=0; + matWerte[i,mfGamma,false]:=0; + matWerte[i,mfIGamma,false]:=0; + end; + for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do + for abl:=false to true do + emWerte[emF,abl]:=0; + lN:=linkerNachbar; rN:=nil; - for Abl:=False to True do - for FI:=low(TFeldinhalt) to high(TFeldinhalt) do - Groessen[Abl,FI]:=0; - dX:=deltaX; - iDX:=1/dX; - x:=0; - pDNMax:=pDNMa; - Groessen[false,fiN]:=1; - Groessen[false,fiGamma]:=1; - Groessen[false,fiiGamma]:=1; + if assigned(lN) then + lN.rN:=self; +end; + +destructor tWertePunkt.destroy; +begin + setlength(matWerte,0); + if assigned(lN) then + lN.rN:=nil; + if assigned(rN) then + rN.lN:=nil; + inherited destroy; end; -procedure TFeld.berechneGammaUndP; // aus aktuellem dPsi/dX und A -var b: byte; - tmp: extended; +procedure tWertePunkt.berechneGammaUndP; // aus aktuellem dPsi/dX und A +var + i: longint; begin - tmp:=0; - for b:=0 to 2 do begin - Groessen[false,TFeldinhalt(Byte(fiPX)+b)]:= - Groessen[false,TFeldinhalt(Byte(fiAX)+b)] + Groessen[false,fidPsiDX]*Byte(b=0); - tmp:=tmp+sqr(Groessen[false,TFeldinhalt(Byte(fiPX)+b)]); + for i:=0 to length(matWerte)-1 do begin + matWerte[i,mfPX,false]:= emWerte[efAX,false] + matWerte[i,mfDPsiDX,false]; + matWerte[i,mfPY,false]:= emWerte[efAY,false]; + matWerte[i,mfPZ,false]:= emWerte[efAZ,false]; + + matWerte[i,mfGamma,false]:= + sqrt( + 1 + + sqr(matWerte[i,mfPX,false]) + + sqr(matWerte[i,mfPY,false]) + + sqr(matWerte[i,mfPZ,false]) + ); + matWerte[i,mfIGamma,false]:= + 1/matWerte[i,mfGamma,false]; end; - Groessen[false,fiGamma]:=sqrt(1+tmp); - Groessen[false,fiiGamma]:=1/Groessen[false,fiGamma]; end; -procedure TFeld.berechneNAbleitung; // Zeitableitung von n berechnen +procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended); // Zeitableitung von n berechnen +var + i: longint; begin (* // dn/dt = -d(n/Gamma p_x)/dx Groessen[true,fiN]:= @@ -141,104 +263,226 @@ begin )*iDX/6; // es wird über die beiden nächsten linken bzw. rechten Nachbarn gemittelt *) - // dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n - Groessen[true,fiN]:= - ( rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*(pDNMax - rN^.Groessen[false,fiPX]) - - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*(-pDNMax - lN^.Groessen[false,fiPX]) - - Groessen[false,fiN]*Groessen[false,fiiGamma]*2*pDNMax)*iDX/2 - // Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung - // und soll die numerische Stabilität bis zu p * dn/dx / n von 2*pDNMax gewährleisten. - // Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit - // alpha = pDNMax * dX + for i:=0 to length(matWerte)-1 do + // dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n + matWerte[i,mfN,true]:= + ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false]) + - ln.matWerte[i,mfN,false] * ln.matWerte[i,mfIGamma,false] * (-pDNMax - rN.matWerte[i,mfPX,false]) + - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * 2 * pDNMax) * iDX/2; + // Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung + // und soll die numerische Stabilität bis zu p * dn/dx / n von 2*pDNMax gewährleisten. + // Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit + // alpha = pDNMax * dX end; -procedure TFeld.berechneAbleitungen(dT: extended); // Zeitableitungen berechnen -var b: byte; +procedure tWertePunkt.berechneAbleitungen(dT,dX,iDX: extended); // Zeitableitungen berechnen +var + i: longint; begin // d2Psi/dxdt = dPhi/dx - dGamma/dx - Groessen[true,fidPsiDX]:= - Groessen[false,fidPhiDX] - - (rN^.Groessen[false,fiGamma] - lN^.Groessen[false,fiGamma])*iDX/2; + for i:=0 to length(matWerte)-1 do + matWerte[i,mfDPsiDX,true]:= + emWerte[efDPhiDX,false] + - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2; + // d3Phi/dx2dt = dn/dt ... hier muss integriert werden: // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + * deltax - Groessen[true,fidPhiDX]:= - lN^.Groessen[true,fidPhiDX] - + (lN^.Groessen[true,fiN] + Groessen[true,fiN])*dX/2; - // d2A/dt2 = - n/gamma p + Laplace(A) ... - for b:=0 to 2 do - Groessen[true,TFeldinhalt(Byte(fidAXDT)+b)]:= - - (Groessen[false,fiN]*Groessen[false,fiiGamma]*Groessen[false,TFeldinhalt(Byte(fiPX)+b)]) - + (rN^.Groessen[false,TFeldinhalt(Byte(fiAX)+b)] - - 2*Groessen[false,TFeldinhalt(Byte(fiAX)+b)] - + lN^.Groessen[false,TFeldinhalt(Byte(fiAX)+b)])*sqr(iDX); - Groessen[true,fidAXDT]:=Groessen[true,fidAXDT] - Groessen[false,fidPhiDX]; // ... - Nabla(phi) + emWerte[efDPhiDX,true]:= + lN.emWerte[efDPhiDX,true]; + for i:=0 to length(matWerte)-1 do + emWerte[efDPhiDX,true]:= + emWerte[efDPhiDX,true] + + (lN.matWerte[i,mfN,true] + matWerte[i,mfN,true])*dX/2; + + // d2A/dt2 = Laplace(A) - Nabla(phi) ... + emWerte[efDAXDT,true]:= + ( rN.emWerte[efAX,false] - 2*emWerte[efAX,false] + lN.emWerte[efAX,false] )*sqr(iDX) + - emWerte[efDPhiDX,false]; + emWerte[efDAYDT,true]:= + ( rN.emWerte[efAY,false] - 2*emWerte[efAY,false] + lN.emWerte[efAY,false] )*sqr(iDX); + emWerte[efDAZDT,true]:= + ( rN.emWerte[efAZ,false] - 2*emWerte[efAZ,false] + lN.emWerte[efAZ,false] )*sqr(iDX); + // ... - n/gamma p + for i:=0 to length(matWerte)-1 do begin + emWerte[efDAXDT,true]:= + emWerte[efDAXDT,true] + - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]); + emWerte[efDAYDT,true]:= + emWerte[efDAYDT,true] + - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]); + emWerte[efDAZDT,true]:= + emWerte[efDAZDT,true] + - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]); + end; + // dA/dt = dA/dt - for b:=0 to 2 do - Groessen[true,TFeldInhalt(Byte(fiAX)+b)]:= - Groessen[false,TFeldInhalt(Byte(fidAXDT)+b)] + - Groessen[true,TFeldInhalt(Byte(fidAXDT)+b)]*dT; + emWerte[efAX,true]:=emWerte[efDAXDT,false] + emWerte[efDAXDT,true]*dT; end; -{ TGitter } +procedure tWertePunkt.liKo(in1,in2: tWertePunkt; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + i: longint; +begin +(* tEMFeldInhalt = ( + efA,efAX,efAY,efAZ, + efDAXDT,efDAYDT,efDAZDT, + efDPhiDX + ); *) + for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde + emWerte[emF,false]:= in1.emWerte[emF,false] + fak2 * in2.emWerte[emF,true]; + +(* tMaterieFeldInhalt = ( + mfN,mfDPsiDX, + mfP,mfPX,mfPY,mfPZ, + mfGamma,mfIGamma + ); *) + for i:=0 to length(matWerte)-1 do // siehe oben + for maF:=mfN to mfDPsiDX do + matWerte[i,maF,false]:= in1.matWerte[i,maF,false] + fak2 * in2.matWerte[i,maF,true]; +end; -procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended; out dTMax: extended); -var i: longint; +procedure tWertePunkt.liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt +var + emF: tEMFeldInhalt; + maF: tMaterieFeldInhalt; + i: longint; begin - berechneAbleitungen(Felder,dT); - dTMax:=dTMaximum; - for i:=0 to groesse+3 do - if Felders[Felder,i].Groessen[true,fiN]<0 then begin - dTMax:=min(dTMax,Abs(Felders[Felder,i].Groessen[false,fiN]/Felders[Felder,i].Groessen[true,fiN])); - end; +(* tEMFeldInhalt = ( + efA,efAX,efAY,efAZ, + efDAXDT,efDAYDT,efDAZDT, + efDPhiDX + ); *) + for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde + emWerte[emF,false]:= + in1.emWerte[emF,false] + + fak2 * in2.emWerte[emF,true] + + fak3 * in3.emWerte[emF,true] + + fak4 * in4.emWerte[emF,true] + + fak5 * in5.emWerte[emF,true]; + +(* tMaterieFeldInhalt = ( + mfN,mfDPsiDX, + mfP,mfPX,mfPY,mfPZ, + mfGamma,mfIGamma + ); *) + for i:=0 to length(matWerte)-1 do // siehe oben + for maF:=mfN to mfDPsiDX do + matWerte[i,maF,false]:= + in1.matWerte[i,maF,false] + + fak2 * in2.matWerte[i,maF,true] + + fak3 * in3.matWerte[i,maF,true] + + fak4 * in4.matWerte[i,maF,true] + + fak5 * in5.matWerte[i,maF,true]; end; -procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended); -var i: longint; +{ tFelder } + +constructor tFelder.create(groesse: longint; dichten,lichter: tMyStringList; parent: tGitter); +var + i,j: longint; + rechts: boolean; begin - for i:=0 to groesse+3 do begin -(* Felders[Felder,i].Groessen[false,fiN]:= - max(Felders[Felder,i].Groessen[false,fiN],0); // n >= 0 *) - Felders[Felder,i].berechneGammaUndP; + inherited create; + par:=parent; + matAnz:=dichten.count; + fillchar(inhalt,sizeof(inhalt),#0); + setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen + inhalt[0]:=tWertePunkt.create(nil,matAnz); + for i:=1 to length(inhalt)-1 do + inhalt[i]:=tWertePunkt.create(inhalt[i-1],matAnz); + for i:=0 to length(inhalt)-1 do begin + parent.kvs.add('x',par.xl+(i-2)/groesse*par.dX); + for j:=0 to dichten.count-1 do + inhalt[i].matWerte[j,mfN,false]:=exprToFloat(false,dichten[j],parent.kvs,nil); + end; + parent.kvs.rem('x'); + for rechts:=false to true do begin + lichters[rechts]:=tMyStringlist.create(nil); + lichters[rechts].text:=lichter.text; end; - Felders[Felder,1].Groessen[false,fiPX]:= // Teilchen werden reflektiert - Abs(Felders[Felder,1].Groessen[false,fiPX]); - Felders[Felder,groesse+2].Groessen[false,fiPX]:= - -Abs(Felders[Felder,groesse+2].Groessen[false,fiPX]); - setzeRaender(Felder); - for i:=1 to groesse+2 do - Felders[Felder,i].berechneNAbleitung; - for i:=1 to groesse+2 do - Felders[Felder,i].berechneAbleitungen(dT); + lichters[false].grep('^links .*'); + lichters[false].subst('^links *',''); + lichters[true].grep('^rechts .*'); + lichters[true].subst('^rechts *',''); end; -procedure TGitter.setzeRaender(Felder: longint); -var FI: TFeldInhalt; +destructor tFelder.destroy; begin - for FI:=fiAX to fiAZ do begin // Vakuumrandbedingungen für das A-Feld - Felders[Felder,0].Groessen[true,FI]:= - (Felders[Felder,1].Groessen[false,FI] - - Felders[Felder,0].Groessen[false,FI])*iDX; - Felders[Felder,groesse+3].Groessen[true,FI]:= - (Felders[Felder,groesse+2].Groessen[false,FI] - - Felders[Felder,groesse+3].Groessen[false,FI])*iDX; - end; // (ein bisschen wird noch reflektiert, vmtl. durch numerische Fehler) - Felders[Felder,0].Groessen[true,fiAy]:= - Felders[Felder,0].Groessen[true,fiAy] + lichters[false].free; + lichters[true].free; + inherited destroy; +end; + +procedure tFelder.setzeRaender(iDX: extended); +var + emF: tEMFeldInhalt; +begin + for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld + inhalt[0].emWerte[emF,true]:= + (inhalt[1].emWerte[emF,false] - + inhalt[0].emWerte[emF,false])*iDX; + inhalt[length(inhalt)-1].emWerte[emF,true]:= + (inhalt[length(inhalt)-2].emWerte[emF,false] - + inhalt[length(inhalt)-1].emWerte[emF,false])*iDX; + end; // (ein bisschen wird noch reflektiert, vmtl. durch numerische Fehler bzw. nicht beachtete numerische Dispersion) + + (* inhalt[0].emWerte[efAY,true]:= + inhalt[0].emWerte[fiAY,true] + (Al(t+dTMaximum)-Al(t))/dTMaximum; Felders[Felder,groesse+3].Groessen[true,fiAy]:= Felders[Felder,groesse+3].Groessen[true,fiAy] - + (Ar(t+dTMaximum)-Ar(t))/dTMaximum; + + (Ar(t+dTMaximum)-Ar(t))/dTMaximum; *) end; -constructor TGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); -var i,j: longint; +procedure tFelder.berechneAbleitungen(dt,dx,iDX,pDNMax: extended); +var + i: longint; +begin + for i:=0 to length(inhalt)-1 do + inhalt[i].berechneGammaUndP; + for i:=0 to matAnz-1 do begin // Teilchen werden reflektiert + inhalt[1].matWerte[i,mfPX,false]:= + abs(inhalt[1].matWerte[i,mfPX,false]); + inhalt[length(inhalt)-2].matWerte[i,mfPX,false]:= + -abs(inhalt[length(inhalt)-2].matWerte[i,mfPX,false]); + end; + + setzeRaender(iDX); + + for i:=1 to length(inhalt)-2 do + inhalt[i].berechneNAbleitung(iDX,pDNMax); + for i:=1 to length(inhalt)-2 do + inhalt[i].berechneAbleitungen(dT,dX,iDX); +end; + +procedure tFelder.liKo(in1,in2: tFelder; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt +var + i: longint; +begin + for i:=0 to length(inhalt)-1 do + inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],fak2); +end; + +procedure tFelder.liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt +var + i: longint; +begin + for i:=0 to length(inhalt)-1 do + inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],in4.inhalt[i],in5.inhalt[i],fak2,fak3,fak4,fak5); +end; + +{ tGitter } + +constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); +var + i: longint; begin inherited create; - Ar:=@nullfunktion; - Al:=@nullfunktion; - Zeitverfahren:=ZV; - groesse:=size; - Prot:=TProtokollant.create(Protokollant,'TGitter'); + zeitverfahren:=ZV; + kvs:=bekannteWerte; + Prot:=TProtokollant.create(Protokollant,'tGitter'); dTMaximum:=deltaT; dX:=deltaX; iDX:=1/dX; @@ -252,252 +496,192 @@ begin zfRungeKuttaVier: Setlength(Felders,5); end{of Case}; xl:=dX/2; - for i:=0 to length(Felders)-1 do begin - setlength(Felders[i],size+4); // auf jeder Seite je zwei Felder Rand extra - for j:=0 to length(felders[i])-1 do begin - felders[i,j]:=TFeld.create(dX,pDNMax); - felders[i,j].x:=xl+j*dX; - knownValues.add('x',felders[i,j].x); - felders[i,j].Groessen[false,fiN]:=n0(felders[i,j].x); - felders[i,j].Groessen[false,fidPsiDX]:=0.0001*Cos(j/(length(felders[i])-1)*pi*15); - end; - felders[i,0].lN:=@felders[i,0]; - for j:=1 to length(felders[i])-1 do begin - felders[i,j].lN:=@felders[i,j-1]; - felders[i,j-1].rN:=@felders[i,j]; - end; - felders[i,length(felders[i])-1].rN:=@felders[i,length(felders[i])-1]; - end; - xr:=xl+dX*(length(Felders[0])-1); + + for i:=0 to length(felders)-1 do + felders[i]:=tFelder.create(size,dichten,lichter,self); + aktuelleFelder:=0; t:=0; end; -destructor TGitter.destroy; -var i,j: longint; +destructor tGitter.destroy; +var + i: longint; begin - for i:=0 to length(Felders)-1 do begin - for j:=0 to length(Felders[i])-1 do - Felders[i,j].free; - setlength(Felders[i],0); - end; - setlength(Felders,0); + for i:=0 to length(Felders)-1 do + felders[i].free; + setlength(felders,0); inherited destroy; end; -procedure TGitter.iteriereSchritt(var dT: extended); +procedure tGitter.iteriereSchritt(dT: extended); var i: longint; - FI: TFeldInhalt; - dTMax: extended; {$IFDEF Zeitschrittueberwachung} - Pro: TProtokollant; + pro: tProtokollant; + j: longint; {$ENDIF} -const zeitObergrenze = 1/4; - zeitSollwert = 1/8; - zeitUntergrenze = 1/16; begin - berechneAbleitungen(aktuelleFelder,dT,dTMax); // y' = y'(t,y(t)) - {$IFDEF Zeitschrittueberwachung} - Pro:=TProtokollant.create(Prot,'iteriereSchritt'); - {$ENDIF} - while dT>=dTMax*zeitObergrenze do begin - dT:=dTMax*zeitSollwert; - berechneAbleitungen(aktuelleFelder,dT,dTMax); // y' = y'(t,y(t)) - {$IFDEF Zeitschrittueberwachung} - Pro.schreibe('+ neues dT: '+floattostr(dT)+' (t='+floattostr(t)+')'); - {$ENDIF} - end; - {$IFDEF Zeitschrittueberwachung} - if dT<1e-30 then begin - for i:=1 to groesse+2 do - if (Felders[aktuelleFelder,i].Groessen[true,fiN]<0) and - (Abs(Felders[aktuelleFelder,i].Groessen[false,fiN]/ - Felders[aktuelleFelder,i].Groessen[true,fiN])<1e-15) then - Prot.schreibe( - floattostr(t)+' '+ - inttostr(i)+' '+ - floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fidPsiDX]* - Felders[aktuelleFelder,i-1].Groessen[false,fiiGamma])+' '+ - floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fiN])+' '+ - floattostr(Felders[aktuelleFelder,i-1].Groessen[true,fiN])+' '+ - floattostr(Felders[aktuelleFelder,i].Groessen[false,fidPsiDX]* - Felders[aktuelleFelder,i].Groessen[false,fiiGamma])+' '+ - floattostr(Felders[aktuelleFelder,i].Groessen[false,fiN])+' '+ - floattostr(Felders[aktuelleFelder,i].Groessen[true,fiN])+' '+ - floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fidPsiDX]* - Felders[aktuelleFelder,i+1].Groessen[false,fiiGamma])+' '+ - floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fiN])+' '+ - floattostr(Felders[aktuelleFelder,i+1].Groessen[true,fiN]), - true); - Pro.destroyall; - halt(1); - end; - {$ENDIF} + felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDX,pDNMax); // y' = y'(t,y(t)) - case Zeitverfahren of - zfEulerVorwaerts: - begin // y(t+dt) = y(t) + y' dt - for i:=0 to groesse+3 do - for FI:=low(TFeldinhalt) to fidPsiDX do - Felders[1-aktuelleFelder,i].Groessen[false,FI]:= - Felders[aktuelleFelder,i].Groessen[false,FI] - + Felders[aktuelleFelder,i].Groessen[true,FI]*dT; - end; + case zeitverfahren of + zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt + felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT); zfRungeKuttaVier: begin - for i:=0 to groesse+3 do // ya = y(t) + y' dt/2 - for FI:=low(TFeldinhalt) to fidPsiDX do - Felders[2,i].Groessen[false,FI]:= - Felders[aktuelleFelder,i].Groessen[false,FI] - + Felders[aktuelleFelder,i].Groessen[true,FI]*dT/2; - berechneAbleitungen(2,dT/2,dTMax); // ya' = y'(t+dt/2,ya) - if dT/2>dTMax*zeitObergrenze then begin - {$IFDEF Zeitschrittueberwachung} - Pro.schreibe('1. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true); - Pro.free; - {$ENDIF} - dT:=dTMax*zeitSollwert; - exit; - end; - for i:=0 to groesse+3 do // yb = y(t) + ya' dt/2 - for FI:=low(TFeldinhalt) to fidPsiDX do - Felders[3,i].Groessen[false,FI]:= - Felders[aktuelleFelder,i].Groessen[false,FI] - + Felders[2,i].Groessen[true,FI]*dT/2; - berechneAbleitungen(3,dT/2,dTMax); // yb' = y'(t+dt/2,yb) - if dT/2>dTMax*zeitObergrenze then begin - {$IFDEF Zeitschrittueberwachung} - Pro.schreibe('2. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true); - Pro.free; - {$ENDIF} - dT:=dTMax*zeitSollwert; - exit; - end; - for i:=0 to groesse+3 do // yc = y(t) + yb' dt - for FI:=low(TFeldinhalt) to fidPsiDX do - Felders[4,i].Groessen[false,FI]:= - Felders[aktuelleFelder,i].Groessen[false,FI] - + Felders[3,i].Groessen[true,FI]*dT; - berechneAbleitungen(4,dT/2,dTMax); // yc' = y'(t+dt,yc) - if dT/2>dTMax*zeitObergrenze then begin - {$IFDEF Zeitschrittueberwachung} - Pro.schreibe('3. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true); - Pro.free; - {$ENDIF} - dT:=dTMax*zeitSollwert; - exit; - end; - for i:=0 to groesse+3 do // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6 - for FI:=low(TFeldinhalt) to fidPsiDX do - Felders[1-aktuelleFelder,i].Groessen[false,FI]:= - Felders[aktuelleFelder,i].Groessen[false,FI] - + (Felders[aktuelleFelder,i].Groessen[true,FI] + - 2*(Felders[2,i].Groessen[true,FI] + - Felders[3,i].Groessen[true,FI]) + - Felders[4,i].Groessen[true,FI])*dT/6; + felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2 + felders[2].berechneAbleitungen(dT/2,dX,iDX,pDNMax); // ya' = y'(t+dt/2,ya) + + felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2 + felders[3].berechneAbleitungen(dT/2,dX,iDX,pDNMax); // yb' = y'(t+dt/2,yb) + + felders[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt + felders[4].berechneAbleitungen(dT/2,dX,iDX,pDNMax); // yc' = y'(t+dt,yc) + + felders[1-aktuelleFelder].liKo( + felders[aktuelleFelder], + felders[aktuelleFelder], + felders[2], + felders[3], + felders[4], + dT/6, + dT/3, + dT/3, + dT/6 + ); // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6 end; end{of case}; t:=t+dT; aktuelleFelder:=1-aktuelleFelder; - if dT=0) and (aD.matInhalt=mfP) then + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin + tmp:=0; + for maF:=mfPX to mfPZ do + tmp:=tmp+sqr(felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,maF,false]); + felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,maF,false]:=sqrt(tmp); + end; + + i:=floor((length(felders[aktuelleFelder].inhalt)-1)/sDX*dX+1); + sX:=xl+(i-1)*sDX; + blockWrite(aD.datei,xl,sizeof(extended)); + blockWrite(aD.datei,sX,sizeof(extended)); + blockWrite(aD.datei,i,sizeof(longint)); sX:=xl-Min(dX,sDX)/2; - for i:=0 to length(Felders[aktuelleFelder])-1 do - while Felders[aktuelleFelder,i].x>=sX do begin - BlockWrite(AD.Datei,Felders[aktuelleFelder,i].Groessen[AD.Ableitung,AD.Inhalt],sizeof(extended)); - sX:=sX+sDX; + cX:=xl; + if aD.teilchen<0 then begin + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin + while cX>=sX do begin + blockWrite(aD.datei,felders[aktuelleFelder].inhalt[i].emWerte[aD.emInhalt,aD.ableitung],sizeof(extended)); + sX:=sX+sDX; + end; + cX:=cX+dX; end; + end + else begin + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin + while cX>=sX do begin + blockWrite(aD.datei,felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,aD.matInhalt,aD.ableitung],sizeof(extended)); + sX:=sX+sDX; + end; + cX:=cX+dX; + end; + end; end; -function TGitter.gibErhaltungsgroessen: string; -var i,j,k: integer; - n: extended; - Pro: TProtokollant; +function tGitter.gibErhaltungsgroessen: string; +var i,j: integer; + ns: tExtendedArray; + Pro: TProtokollant; begin - Pro:=TProtokollant.create(Prot,'gibErhaltungsgroessen'); - n:=0; - for i:=0 to groesse+3 do - n:=n+Felders[aktuelleFelder,i].Groessen[false,fiN]; - for i:=1 to groesse+2 do - if abs(Felders[aktuelleFelder,i].Groessen[false,fiPx]* - (Felders[aktuelleFelder,i+1].Groessen[false,fiN]-Felders[aktuelleFelder,i-1].Groessen[false,fiN])/ - Max(Felders[aktuelleFelder,i+1].Groessen[false,fiN]+Felders[aktuelleFelder,i-1].Groessen[false,fiN],1e-10)) - > pDNMax then begin - if pDNMaxDynamisch then begin - pDNMax:= - 2* - abs(Felders[aktuelleFelder,i].Groessen[false,fiPx]* - (Felders[aktuelleFelder,i+1].Groessen[false,fiN]-Felders[aktuelleFelder,i-1].Groessen[false,fiN])/ - Max(Felders[aktuelleFelder,i+1].Groessen[false,fiN]+Felders[aktuelleFelder,i-1].Groessen[false,fiN],1e-10)); - for j:=0 to length(felders)-1 do - for k:=0 to length(felders[j])-1 do - felders[j,k].pDNMax:=pDNMax; - Pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')'); - end - else - if pDNMax>0 then begin - Pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+ - 'T), die numerische Stabilität ist nicht mehr gewährleistet!',true); - Pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true); - Pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true); - pDNMax:=-1; - end; + pro:=tProtokollant.create(prot,'gibErhaltungsgroessen'); + setlength(ns,felders[aktuelleFelder].matAnz); + result:=''; + for i:=0 to length(ns)-1 do + ns[i]:=0; + for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do + for j:=0 to length(ns)-1 do + ns[j]:=ns[j] + felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]; + for i:=0 to length(ns)-1 do begin + result:=result + ' n['+inttostr(i)+']='+floattostr(ns[i]); + {$IFDEF Dichteueberwachung} + if ns[i]>1000 then begin + errorCode:=2; + pro.schreibe(' n > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')'); end; - result:='n='+floattostr(n); - {$IFDEF Dichteueberwachung} - if n>1000 then begin - errorCode:=2; - Pro.schreibe(' n > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')'); + {$ENDIF} end; - {$ENDIF} - Pro.free; + delete(result,1,1); + pro.free; +end; + +procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant); +var + i,j: longint; +begin + for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do + for j:=0 to felders[aktuelleFelder].matAnz-1 do + if abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]* + (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/ + max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10)) + > pDNMax then begin + if pDNMaxDynamisch then begin + pDNMax:= + 2* + abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]* + (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/ + Max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10)); + pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')'); + end + else + if pDNMax>0 then begin + pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+ + 'T), die numerische Stabilität ist nicht mehr gewährleistet!',true); + pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true); + pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true); + pDNMax:=-1; + end; + end; end; { tSimulation } @@ -507,11 +691,10 @@ var ifile: tMyStringlist; zeitverfahren: tZeitverfahren; s,t,aSuffix,aPrefix: string; - deltaX,deltaT,endzeit,breite,sDT,pDNMax: extended; - ausgabeDateien: array of tAusgabeDatei; - abl: boolean; - fi: tFeldinhalt; + deltaX,deltaT,breite,pDNMax: extended; i,j: longint; + kvs: tKnownValues; + dichten,lichter: tMyStringlist; begin inherited create; prot:=tProtokollant.create(Protokollant,'tSimulation'); @@ -527,25 +710,35 @@ begin // Standardeinstellungen Bereich 'allgemein' zeitverfahren:=zfRungeKuttaVier; deltaX:=1e-2; - kvs.add('λ',1/deltaX); + kvs.add('λ',1); + kvs.add('T',1); + kvs.add('dX',deltaX); deltaT:=-1; sDT:=-1; - endzeit:=100; + sDX:=-1; + tEnde:=100; breite:=10.0; pDNMax:=-1; + fortschrittsAnzeige:=false; // Standardeinstellungen Bereich 'ausgaben' aPrefix:=extractfilepath(paramstr(0)); aSuffix:='.dat'; setlength(ausgabeDateien,0); - while ifile.readln(s) do begin + // Standardeinstellungen Breich 'lichtVon...' + lichter:=tMyStringlist.create(prot); + + // Standardeinstellungen Breich 'teilchen...' + dichten:=tMyStringlist.create(prot); + + while ifile.readln(s) do begin writeln(s); if s='allgemein' then begin repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich allgemein!',true); halt(1); - end; + end; writeln(s); if s='allgemeinEnde' then break; if s='runge-Kutta-4' then begin Zeitverfahren:=zfRungeKuttaVier; @@ -557,7 +750,7 @@ begin end; if startetMit('ortsschritt',s) then begin deltaX:=exprtofloat(false,s,kvs,nil); - kvs.add('λ',1/deltaX); + kvs.add('dX',deltaX); continue; end; if startetMit('zeitschritt',s) then begin @@ -570,13 +763,21 @@ begin continue; end; if startetMit('zeit',s) then begin - endzeit:=exprtofloat(false,s,kvs,nil); + tEnde:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('breite',s) then begin breite:=exprtofloat(false,s,kvs,nil); continue; end; + if s='mit Fortschrittsanzeige' then begin + fortschrittsAnzeige:=true; + continue; + end; + if s='ohne Fortschrittsanzeige' then begin + fortschrittsAnzeige:=false; + continue; + end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich allgemein!',true); halt(1); until false; @@ -602,20 +803,16 @@ begin sDT:=exprtofloat(false,s,kvs,nil); continue; end; + if startetMit('ortsschritt',s) then begin + sDX:=exprtofloat(false,s,kvs,nil); + continue; + end; if startetMit('felder',s) then begin s:=s+','; while s<>'' do begin t:=erstesArgument(s,','); - for abl:=false to true do - for fI:=low(tFeldInhalt) to high(tFeldInhalt) do - if not (abl and (fI in [fiGamma,fiP])) then - if t=copy('d',1,byte(abl))+optionNamen[fI] then begin - setlength(ausgabedateien,length(ausgabedateien)+1); - ausgabedateien[length(ausgabedateien)-1].inhalt:=fI; - ausgabedateien[length(ausgabedateien)-1].ableitung:=abl; - ausgabedateien[length(ausgabedateien)-1].name:=aPrefix+t+aSuffix; - assignFile(ausgabedateien[length(ausgabedateien)-1].datei,ausgabedateien[length(ausgabedateien)-1].name); - end; + setlength(ausgabeDateien,length(ausgabeDateien)+1); + ausgabeDateien[length(ausgabeDateien)-1]:=tAusgabeDatei.create(t,aPrefix,aSuffix,prot); end; continue; end; @@ -638,26 +835,79 @@ begin deltaT:=deltaX/10; if sDT<0 then sDT:=deltaT; + if sDX<0 then + sDX:=deltaX; for i:=0 to length(ausgabedateien)-1 do begin rewrite(ausgabedateien[i].datei,1); j:=0; BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); - j:=floor(Endzeit/sDT+1); + j:=floor(tEnde/sDT+1); BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); end; ifile.free; - gitter:=tGitter.create(Breite,deltaT,deltaX,pDNMax,@Anfangsdichte,Zeitverfahren,Prot); + gitter:=tGitter.create(round(breite/deltaX),deltaT,deltaX,pDNMax,kvs,dichten,lichter,zeitverfahren,prot); + dichten.free; + lichter.free; end; destructor tSimulation.destroy; begin gitter.free; - kvs.free; prot.free; inherited destroy; end; +function tSimulation.iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? +var + i: longint; + c: char; +begin + result:=false; + + zeitPhysik:=zeitPhysik-now; + if errorCode<2 then + gitter.iteriereSchritt(dT); + zeitPhysik:=zeitPhysik+now; + zeitDatei:=zeitDatei-now; + while gitter.t>=sT do begin + sT:=sT+sDT; + for i:=0 to length(Ausgabedateien)-1 do + gitter.macheAusgabe(ausgabeDateien[i],sDX); + end; + zeitDatei:=zeitDatei+now; + + gitter.iteriereSchritt(dT); + + if fortschrittsAnzeige then begin + if (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde)) or keypressed then begin + if keypressed then c:=readkey + else c:=#0; + case c of + #27,'q': begin + errorCode:=3; + exit; + end; + ' ': begin + writeln(' ... Pause (beliebige Taste drücken um fortzufahren) ...'); + readkey; + writeln(' ... weiter geht''s ...'); + end; + else begin + Prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true); + Prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true); + Prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true); + Prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true); + if errorCode<2 then + Prot.schreibe(gitter.gibErhaltungsgroessen); + end; + end{of case}; + end; + end; + + result:=gitter.t>=tEnde; +end; + end. diff --git a/Plasmapropagation b/Plasmapropagation index 013329a..98472f4 100755 Binary files a/Plasmapropagation and b/Plasmapropagation differ diff --git a/Plasmapropagation.lpr b/Plasmapropagation.lpr index b861d2e..5c914d5 100644 --- a/Plasmapropagation.lpr +++ b/Plasmapropagation.lpr @@ -28,13 +28,10 @@ type procedure TPlasmapropagation.DoRun; var - Breite,i,j: longint; simulation: tSimulation; - start,zeitDatei,zeitPhysik: extended; - FI: TFeldInhalt; - Abl: Boolean; - Prot: TProtokollant; - c: char; + start,zeitPhysik,zeitDatei: extended; + prot: tProtokollant; + s,t,u: string; begin prot:=tProtokollant.create('error'); @@ -48,71 +45,29 @@ begin zeitDatei:=0; zeitPhysik:=0; - sT:=-Min(deltaT,sDT)/2; - simulation:=tSimulation.create(paramstr(1)); -// Gitter.Al:=@InputFeld; - while Gitter.t=sT do begin - sT:=sT+sDT; - for j:=0 to length(Ausgabedateien)-1 do - Gitter.macheAusgabe(Ausgabedateien[j],sDX); - end; - zeitDatei:=zeitDatei+now; - end; + simulation:=tSimulation.create(paramstr(1),prot); + + while simulation.iteriereSchritt(start,zeitPhysik,zeitDatei) do ; case errorCode of - 2: Prot.schreibe('Simulation wurde auf halbem Wege abgebrochen wegen Überlaufs.',true); - 3: Prot.schreibe('Simulation wurde auf halbem Wege abgebrochen wegen Benutzereingriffs.',true); + 2: prot.schreibe('Simulation wurde auf halbem Wege abgebrochen wegen Überlaufs.',true); + 3: prot.schreibe('Simulation wurde auf halbem Wege abgebrochen wegen Benutzereingriffs.',true); end{of case}; - for i:=0 to length(Ausgabedateien)-1 do - CloseFile(Ausgabedateien[i].Datei); + simulation.free; + prot.schreibe('fertig!',true); - Prot.schreibe('fertig!',true); s:=timetostr(now-start); t:=timetostr(zeitDatei); u:=timetostr(zeitPhysik); while length(s) - - - + + + - - - - + + + - + @@ -31,195 +30,198 @@ - + - - - + + + + - - - - - - - + + + + + + + + + + - - - - - + + + + + - - - - + + + + + + - - - - - - - - - - - + + + - - - + + + - + + + + + + + + - - - + + - + - + - + - + - + - - + + - - + - - + + - - + + - - + + - - + + - - + + - - + + - + - + - - + + - - + + - - + + - - + + - - + + - - + + - + - + - - + + - - + + - - + + - + - - + + diff --git a/error b/error deleted file mode 100644 index 166a407..0000000 --- a/error +++ /dev/null @@ -1,521 +0,0 @@ - Simulationsbox - 1500 Felder breit (1.5 lambda) - 1 T lang (etwa 10000 Schritte) - 0.0001 Zeitschritt und 0.001 Ortsschritt - Mache Ausgabe fiAX. - Mache Ausgabe fiAY. - Mache Ausgabe fidAYDT. - Mache Ausgabe fiN. - Mache Ausgabe fidPhiDX. - Mache Ausgabe fidPsiDX. - Mache Ausgabe fiGamma. - Maximales px * dn/dx / n wird automatisch ermittelt. - Iterationsverfahren: Runge-Kutta -TGitter.gibErhaltungsgroessen Neues maximales px * dn/dx / n: 4.90867701035798E-13 (t=0.000025) -TGitter.gibErhaltungsgroessen Neues maximales px * dn/dx / n: 1.9210083367696E-6 (t=0.000025) -TGitter.gibErhaltungsgroessen Neues maximales px * dn/dx / n: 6.95677493143965E-6 (t=0.000025) -TGitter.gibErhaltungsgroessen Neues maximales px * dn/dx / n: 0.00001391665818 (t=0.000025) - 1% (t=0.0099875T) - 4s (0.93572594126268) - ETA: 6min 12s - aktueller Zeitschritt: 0.0000125T - n=1100.1504001997 - 2% (t=0.0199875T) - 7s (0.93845056986265) - ETA: 6min 4s - aktueller Zeitschritt: 0.0000125T - n=1100.15040039965 - 3% (t=0.0299875T) - 11s (0.9396443307214) - ETA: 5min 59s - aktueller Zeitschritt: 0.0000125T - n=1100.1504005996 - 4% (t=0.04T) - 15s (0.93995505725415) - ETA: 5min 55s - aktueller Zeitschritt: 0.0000125T - n=1100.1504007998 - 5% (t=0.05T) - 18s (0.94021909059327) - ETA: 5min 51s - aktueller Zeitschritt: 0.0000125T - n=1100.15040099976 - 6% (t=0.06T) - 22s (0.94028101686857) - ETA: 5min 47s - aktueller Zeitschritt: 0.0000125T - n=1100.15040119971 - 7% (t=0.07T) - 26s (0.94036230559588) - ETA: 5min 43s - aktueller Zeitschritt: 0.0000125T - n=1100.15040139966 - 8% (t=0.08T) - 30s (0.94045357808056) - ETA: 5min 39s - aktueller Zeitschritt: 0.0000125T - n=1100.15040159961 - 9% (t=0.09T) - 33s (0.94041907955508) - ETA: 5min 36s - aktueller Zeitschritt: 0.0000125T - n=1100.15040179956 - 10% (t=0.1T) - 37s (0.93998029343796) - ETA: 5min 32s - aktueller Zeitschritt: 0.0000125T - n=1100.15040199952 - 11% (t=0.11T) - 41s (0.9391274466381) - ETA: 5min 29s - aktueller Zeitschritt: 0.0000125T - n=1100.15040219947 - 12% (t=0.12T) - 44s (0.93928366492497) - ETA: 5min 25s - aktueller Zeitschritt: 0.0000125T - n=1100.15040239942 - 13% (t=0.13T) - 48s (0.93936119390127) - ETA: 5min 21s - aktueller Zeitschritt: 0.0000125T - n=1100.15040259938 - 14% (t=0.14T) - 52s (0.93947196769724) - ETA: 5min 17s - aktueller Zeitschritt: 0.0000125T - n=1100.15040279933 - 15% (t=0.1499875T) - 55s (0.93947754328833) - ETA: 5min 14s - aktueller Zeitschritt: 0.0000125T - n=1100.15040299903 - 16% (t=0.1599875T) - 59s (0.93955884945007) - ETA: 5min 10s - aktueller Zeitschritt: 0.0000125T - n=1100.15040319899 - 17% (t=0.1699875T) - 1min 3s (0.93965263560525) - ETA: 5min 6s - aktueller Zeitschritt: 0.0000125T - n=1100.15040339894 - 18% (t=0.1799875T) - 1min 6s (0.9397448828285) - ETA: 5min 2s - aktueller Zeitschritt: 0.0000125T - n=1100.15040359889 - 19% (t=0.1899875T) - 1min 10s (0.93915432836246) - ETA: 4min 59s - aktueller Zeitschritt: 0.0000125T - n=1100.15040379885 - 20% (t=0.1999875T) - 1min 14s (0.93900073347529) - ETA: 4min 55s - aktueller Zeitschritt: 0.0000125T - n=1100.15040399879 - 21% (t=0.2099875T) - 1min 17s (0.93909659887595) - ETA: 4min 51s - aktueller Zeitschritt: 0.0000125T - n=1100.15040419873 - 22% (t=0.2199875T) - 1min 21s (0.93917337012531) - ETA: 4min 48s - aktueller Zeitschritt: 0.0000125T - n=1100.15040439866 - 23% (t=0.2299875T) - 1min 25s (0.93758208132779) - ETA: 4min 44s - aktueller Zeitschritt: 0.0000125T - n=1100.15040459857 - 24% (t=0.2399875T) - 1min 29s (0.93768875699546) - ETA: 4min 41s - aktueller Zeitschritt: 0.0000125T - n=1100.15040479846 - 25% (t=0.2499875T) - 1min 32s (0.9378121327816) - ETA: 4min 37s - aktueller Zeitschritt: 0.0000125T - n=1100.15040499831 - 26% (t=0.2599875T) - 1min 36s (0.93788794049542) - ETA: 4min 33s - aktueller Zeitschritt: 0.0000125T - n=1100.15040519814 - 27% (t=0.2699875T) - 1min 39s (0.93470907399499) - ETA: 4min 29s - aktueller Zeitschritt: 0.0000125T - n=1100.15040539794 - 28% (t=0.2799875T) - 1min 43s (0.93392779228268) - ETA: 4min 24s - aktueller Zeitschritt: 0.0000125T - n=1100.15040559772 - 29% (t=0.2899875T) - 1min 46s (0.9341559070417) - ETA: 4min 21s - aktueller Zeitschritt: 0.0000125T - n=1100.15040579748 - 30% (t=0.2999875T) - 1min 50s (0.93429738634912) - ETA: 4min 17s - aktueller Zeitschritt: 0.0000125T - n=1100.15040599723 - 31% (t=0.3099875T) - 1min 54s (0.93446650346284) - ETA: 4min 13s - aktueller Zeitschritt: 0.0000125T - n=1100.15040619699 - 32% (t=0.3199875T) - 1min 57s (0.93465115233278) - ETA: 4min 10s - aktueller Zeitschritt: 0.0000125T - n=1100.15040639676 - 33% (t=0.3299875T) - 2min 1s (0.93480925657284) - ETA: 4min 6s - aktueller Zeitschritt: 0.0000125T - n=1100.15040659655 - 34% (t=0.3399875T) - 2min 5s (0.93496938818957) - ETA: 4min 2s - aktueller Zeitschritt: 0.0000125T - n=1100.15040679637 - 35% (t=0.3499875T) - 2min 8s (0.93499199359356) - ETA: 3min 59s - aktueller Zeitschritt: 0.0000125T - n=1100.15040699623 - 36% (t=0.3599875T) - 2min 12s (0.93513974146828) - ETA: 3min 55s - aktueller Zeitschritt: 0.0000125T - n=1100.15040719612 - 37% (t=0.3699875T) - 2min 16s (0.93462522712332) - ETA: 3min 51s - aktueller Zeitschritt: 0.0000125T - n=1100.15040739603 - 38% (t=0.3799875T) - 2min 20s (0.93477500648) - ETA: 3min 48s - aktueller Zeitschritt: 0.0000125T - n=1100.15040759598 - 39% (t=0.3899875T) - 2min 23s (0.93418072020113) - ETA: 3min 43s - aktueller Zeitschritt: 0.0000125T - n=1100.15040779595 - 40% (t=0.3999875T) - 2min 26s (0.93434368089496) - ETA: 3min 39s - aktueller Zeitschritt: 0.0000125T - n=1100.15040799592 - 41% (t=0.4099875T) - 2min 29s (0.93438719449992) - ETA: 3min 35s - aktueller Zeitschritt: 0.0000125T - n=1100.1504081959 - 42% (t=0.4199875T) - 2min 33s (0.93448808312355) - ETA: 3min 31s - aktueller Zeitschritt: 0.0000125T - n=1100.15040839587 - 43% (t=0.4299875T) - 2min 37s (0.93459425789043) - ETA: 3min 28s - aktueller Zeitschritt: 0.0000125T - n=1100.15040859582 - 44% (t=0.4399875T) - 2min 40s (0.93472965099255) - ETA: 3min 24s - aktueller Zeitschritt: 0.0000125T - n=1100.15040879575 - 45% (t=0.4499875T) - 2min 44s (0.93483912389575) - ETA: 3min 20s - aktueller Zeitschritt: 0.0000125T - n=1100.15040899565 - 46% (t=0.4599875T) - 2min 48s (0.93495962549421) - ETA: 3min 17s - aktueller Zeitschritt: 0.0000125T - n=1100.15040919551 - 47% (t=0.4699875T) - 2min 51s (0.93495055947613) - ETA: 3min 13s - aktueller Zeitschritt: 0.0000125T - n=1100.15040939535 - 48% (t=0.4799875T) - 2min 55s (0.93509540698992) - ETA: 3min 10s - aktueller Zeitschritt: 0.0000125T - n=1100.15040959517 - 49% (t=0.4899875T) - 2min 59s (0.93522625662536) - ETA: 3min 6s - aktueller Zeitschritt: 0.0000125T - n=1100.15040979497 - 50% (t=0.4999875T) - 3min 2s (0.93534449434932) - ETA: 3min 2s - aktueller Zeitschritt: 0.0000125T - n=1100.15040999477 - 51% (t=0.5099875T) - 3min 6s (0.93508263291742) - ETA: 2min 59s - aktueller Zeitschritt: 0.0000125T - n=1100.15041019456 - 52% (t=0.5199875T) - 3min 10s (0.93509288526884) - ETA: 2min 55s - aktueller Zeitschritt: 0.0000125T - n=1100.15041039437 - 53% (t=0.5299875T) - 3min 13s (0.93488834098026) - ETA: 2min 52s - aktueller Zeitschritt: 0.0000125T - n=1100.1504105942 - 54% (t=0.5399875T) - 3min 17s (0.93499188128487) - ETA: 2min 48s - aktueller Zeitschritt: 0.0000125T - n=1100.15041079406 - 55% (t=0.5499875T) - 3min 21s (0.93480126139187) - ETA: 2min 44s - aktueller Zeitschritt: 0.0000125T - n=1100.15041099395 - 56% (t=0.5599875T) - 3min 25s (0.93490491597985) - ETA: 2min 41s - aktueller Zeitschritt: 0.0000125T - n=1100.15041119387 - 57% (t=0.5699875T) - 3min 28s (0.93500136747388) - ETA: 2min 37s - aktueller Zeitschritt: 0.0000125T - n=1100.15041139383 - 58% (t=0.5799875T) - 3min 32s (0.93465203973783) - ETA: 2min 34s - aktueller Zeitschritt: 0.0000125T - n=1100.15041159381 - 59% (t=0.5899875T) - 3min 36s (0.93472936944875) - ETA: 2min 30s - aktueller Zeitschritt: 0.0000125T - n=1100.15041179381 - 60% (t=0.6T) - 3min 39s (0.9348087023095) - ETA: 2min 26s - aktueller Zeitschritt: 0.0000125T - n=1100.15041199407 - 61% (t=0.61T) - 3min 43s (0.93491070711469) - ETA: 2min 23s - aktueller Zeitschritt: 0.0000125T - n=1100.15041219408 - 62% (t=0.62T) - 3min 46s (0.93462517792379) - ETA: 2min 19s - aktueller Zeitschritt: 0.0000125T - n=1100.15041239408 - 63% (t=0.63T) - 3min 50s (0.93465224645404) - ETA: 2min 15s - aktueller Zeitschritt: 0.0000125T - n=1100.15041259406 - 64% (t=0.64T) - 3min 53s (0.93472351649562) - ETA: 2min 11s - aktueller Zeitschritt: 0.0000125T - n=1100.15041279402 - 65% (t=0.65T) - 3min 56s (0.9346948061847) - ETA: 2min 7s - aktueller Zeitschritt: 0.0000125T - n=1100.15041299395 - 66% (t=0.66T) - 3min 59s (0.93474317500594) - ETA: 2min 3s - aktueller Zeitschritt: 0.0000125T - n=1100.15041319385 - 67% (t=0.67T) - 4min 2s (0.93478852332337) - ETA: 1min 59s - aktueller Zeitschritt: 0.0000125T - n=1100.15041339372 - 68% (t=0.68T) - 4min 5s (0.93477508176243) - ETA: 1min 55s - aktueller Zeitschritt: 0.0000125T - n=1100.15041359357 - 69% (t=0.69T) - 4min 9s (0.93486059208118) - ETA: 1min 52s - aktueller Zeitschritt: 0.0000125T - n=1100.1504137934 - 70% (t=0.7T) - 4min 13s (0.93484459753493) - ETA: 1min 48s - aktueller Zeitschritt: 0.0000125T - n=1100.15041399322 - 71% (t=0.71T) - 4min 16s (0.93492848469232) - ETA: 1min 45s - aktueller Zeitschritt: 0.0000125T - n=1100.15041419304 - 72% (t=0.72T) - 4min 20s (0.93471191856157) - ETA: 1min 41s - aktueller Zeitschritt: 0.0000125T - n=1100.15041439287 - 73% (t=0.73T) - 4min 24s (0.9348012495454) - ETA: 1min 38s - aktueller Zeitschritt: 0.0000125T - n=1100.15041459273 - 74% (t=0.74T) - 4min 27s (0.9348810582785) - ETA: 1min 34s - aktueller Zeitschritt: 0.0000125T - n=1100.15041479261 - 75% (t=0.75T) - 4min 31s (0.93466154687072) - ETA: 1min 30s - aktueller Zeitschritt: 0.0000125T - n=1100.15041499252 - 76% (t=0.76T) - 4min 35s (0.93474800649981) - ETA: 1min 27s - aktueller Zeitschritt: 0.0000125T - n=1100.15041519246 - 77% (t=0.77T) - 4min 39s (0.93464517725866) - ETA: 1min 23s - aktueller Zeitschritt: 0.0000125T - n=1100.15041539244 - 78% (t=0.78T) - 4min 42s (0.93471589587541) - ETA: 1min 20s - aktueller Zeitschritt: 0.0000125T - n=1100.15041559244 - 79% (t=0.79T) - 4min 46s (0.93479319538545) - ETA: 1min 16s - aktueller Zeitschritt: 0.0000125T - n=1100.15041579245 - 80% (t=0.8T) - 4min 50s (0.93487046513783) - ETA: 1min 12s - aktueller Zeitschritt: 0.0000125T - n=1100.15041599248 - 81% (t=0.81T) - 4min 53s (0.93494272921598) - ETA: 1min 9s - aktueller Zeitschritt: 0.0000125T - n=1100.15041619251 - 82% (t=0.82T) - 4min 57s (0.93502453059895) - ETA: 1min 5s - aktueller Zeitschritt: 0.0000125T - n=1100.15041639253 - 83% (t=0.83T) - 5min 1s (0.93509745330369) - ETA: 1min 2s - aktueller Zeitschritt: 0.0000125T - n=1100.15041659253 - 84% (t=0.84T) - 5min 4s (0.93517156164012) - ETA: 58s - aktueller Zeitschritt: 0.0000125T - n=1100.1504167925 - 85% (t=0.85T) - 5min 8s (0.9352390216625) - ETA: 54s - aktueller Zeitschritt: 0.0000125T - n=1100.15041699245 - 86% (t=0.86T) - 5min 12s (0.9352517564134) - ETA: 51s - aktueller Zeitschritt: 0.0000125T - n=1100.15041719236 - 87% (t=0.87T) - 5min 16s (0.9351394704682) - ETA: 47s - aktueller Zeitschritt: 0.0000125T - n=1100.15041739224 - 88% (t=0.88T) - 5min 19s (0.93520735055799) - ETA: 44s - aktueller Zeitschritt: 0.0000125T - n=1100.1504175921 - 89% (t=0.89T) - 5min 23s (0.93512915953451) - ETA: 40s - aktueller Zeitschritt: 0.0000125T - n=1100.15041779194 - 90% (t=0.9T) - 5min 27s (0.93519826581303) - ETA: 36s - aktueller Zeitschritt: 0.0000125T - n=1100.15041799177 - 91% (t=0.91T) - 5min 31s (0.93483492969017) - ETA: 33s - aktueller Zeitschritt: 0.0000125T - n=1100.1504181916 - 92% (t=0.92T) - 5min 34s (0.93490573860459) - ETA: 29s - aktueller Zeitschritt: 0.0000125T - n=1100.15041839144 - 93% (t=0.93T) - 5min 38s (0.93496816016958) - ETA: 25s - aktueller Zeitschritt: 0.0000125T - n=1100.1504185913 - 94% (t=0.94T) - 5min 42s (0.9348235620777) - ETA: 22s - aktueller Zeitschritt: 0.0000125T - n=1100.15041879118 - 95% (t=0.95T) - 5min 45s (0.93488922414053) - ETA: 18s - aktueller Zeitschritt: 0.0000125T - n=1100.1504189911 - 96% (t=0.96T) - 5min 49s (0.93492794655297) - ETA: 15s - aktueller Zeitschritt: 0.0000125T - n=1100.15041919104 - 97% (t=0.97T) - 5min 53s (0.93495188615352) - ETA: 11s - aktueller Zeitschritt: 0.0000125T - n=1100.15041939102 - 98% (t=0.98T) - 5min 56s (0.93481669414105) - ETA: 7s - aktueller Zeitschritt: 0.0000125T - n=1100.15041959102 - 99% (t=0.99T) - 6min 0s (0.93488083192356) - ETA: 4s - aktueller Zeitschritt: 0.0000125T - n=1100.15041979103 - 100% (t=1T) - 6min 4s (0.93494126481698) - ETA: 0s - aktueller Zeitschritt: 0.0000125T - n=1100.15041999106 - fertig! - Das hat 6min 4s gedauert, - davon 23s für Dateizugriffe - und nur 5min 29s für die eigentliche Physik! diff --git a/input.plap b/input.plap index 2e5ecc4..d2debd1 100644 --- a/input.plap +++ b/input.plap @@ -2,19 +2,21 @@ allgemein runge-Kutta-4 - ortsschritt 10^-2 λ - zeitschritt 10^-3 T - zeit 100 T - breite 10 λ + ortsschritt 10^-2 * λ + zeitschritt 10^-3 * T + zeit 100 * T + breite 10 * λ + mit Fortschrittsanzeige allgemeinEnde ausgaben suffix _test.dat - felder AX,AY,dAYDT,N,dPhiDX,dPsiDX,Gamma + felder AX,AY,dAYDT,dPhiDX +# felder AX,AY,dAYDT,N,dPhiDX,dPsiDX,Gamma ausgabenEnde -!set $tFwhm: 5 T -!set $tMitte: 20 T +!setze $tFwhm: 5 T +!setze $tMitte: 20 T lichtVonLinks 0.5 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * Sin(ω*t) -- cgit v1.2.3-70-g09d2