unit Physikunit; {$mode objfpc}{$H+} {$IFNDEF UNIX} {$ERROR This program can be compiled only on/for Unix/Linux based systems.} {$ENDIF} {$DEFINE Zeitschrittueberwachung} {$DEFINE Dichteueberwachung} interface uses Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix; type tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier,zfRungeKuttaZehn,zfRungeKuttaZwoelf,zfRungeKuttaVierzehn); tVerteilungsfunktion = function(x: extended): extended; tEMFeldInhalt = ( efA,efAX,efAY,efAZ, efDAXDT,efDAYDT,efDAZDT, efDPhiDX ); tMaterieFeldInhalt = ( mfN,mfDPsiDX, mfP,mfPX,mfPY,mfPZ, mfGamma,mfIGamma ); tGitter = class; { tAusgabeDatei } tAusgabeDatei = class private ableitung: boolean; emInhalt: tEMFeldInhalt; matInhalt: tMaterieFeldInhalt; teilchen,nNum,tCnt: longint; // wenn teilchen < 0, dann EM-Feld pre,suf: string; datei: file; pro: tProtokollant; public zeitAnz: longint; sDT: extended; constructor create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); destructor destroy; override; function dump: string; procedure schreibeKopf; procedure speichereWerte(gitter: tGitter; sDX: extended); end; { tTeilchenSpezies } tTeilchenSpezies = class private ortsGrenzen: array of extended; dichteStuecke: array of string; public nMin,nMax, // minimale und maximale Massendichte in nc spezifischeLadung: extended; // q/m in e/me constructor create; destructor destroy; override; procedure dump(prot: tProtokollant; prefix: string); function gibDichte(x: extended; kvs: tKnownValues): extended; // Massendichte procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant); end; { tWertePunkt } 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; constructor create(linkerNachbar: tWertePunkt; materien: longint); destructor destroy; override; procedure berechneGammaUndP; procedure berechneNAbleitung(iDX,pDNMax: extended); overload; procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload; procedure berechneAbleitungen(dX,iDX: extended); procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt procedure liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); overload; procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload; procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6: tWertePunkt; fak2,fak3,fak4,fak5,fak6: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; function maxDT: extended; end; { TFelder } tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen private matAnz: longint; spezLadungen: tExtendedArray; lichters: array[boolean] of tMyStringlist; procedure setzeRaender(dT,iDX,iDT: extended); public inhalt: array of tWertePunkt; par: tGitter; constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); destructor destroy; override; procedure berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended); procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt procedure liKo(in1,in2,in3: tFelder; fak2,fak3: extended); overload; procedure liKo(in1,in2,in3,in4: tFelder; fak2,fak3,fak4: extended); overload; procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6: tFelder; fak2,fak3,fak4,fak5,fak6: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7: tFelder; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; function maxDT: extended; end; { tGitter } tGitter = class(tObject) private prot: tProtokollant; zeitverfahren: tZeitverfahren; pDNMaxDynamisch: boolean; dX,iDX,pDNMax,dTMaximum,xl,t: extended; kvs: tKnownValues; procedure diffusionsTermAnpassen(pro: tProtokollant); function pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert? public aktuelleFelder: longint; felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); destructor destroy; override; procedure iteriereSchritt(var dT: extended); function dumpErhaltungsgroessen: boolean; procedure berechne(was: char; teilchen: longint); end; { tSimulation } tSimulation = class(tObject) private prot: tProtokollant; gitter: tGitter; dT,tEnde,sT,sDT,sDX: extended; fortschrittsAnzeige: boolean; ausgabeDateien: array of tAusgabeDatei; public gotSighup: boolean; constructor create(inName: string; protokollant: tProtokollant; name: string); destructor destroy; override; function iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? end; procedure SignalCapture(signal : longint); cdecl; implementation const 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' ); var sighupSimulationen: array of tSimulation; { tAusgabeDatei } constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); var emF: tEMFeldInhalt; maF: tMaterieFeldInhalt; num: longint; abl,gef: boolean; begin inherited create; pro:=tProtokollant.create(prot,nam); 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)))-1; delete(feldName,num,length(feldName)); end else teilchen:=0; emInhalt:=efA; matInhalt:=mfN; feldName:=ansiUpperCase(feldName); nNum:=0; // Header 0 wurde also noch nicht geschrieben zeitAnz:=-1; tCnt:=-1; 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 pro.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true); pro.destroyAll; halt(1); end; if teilchen>=0 then pre:=matFeldNamen[maF]+inttostr(teilchen+1) else pre:=emFeldNamen[emF]; if ableitung then pre:='d'+pre; pre:=prefix+pre; suf:=suffix; end; destructor tAusgabeDatei.destroy; begin if (tCnt<>zeitAnz) and ((nNum<>0) or (tCnt<>-1)) then begin pro.schreibe('Falsche Anzahl an Zeitschritten in '+pre+'-'+inttostr(nNum-1)+suf+' geschrieben ('+inttostr(tCnt)+' statt '+inttostr(zeitAnz)+')!',true); pro.destroyall; halt(1); end; inherited destroy; end; function tAusgabeDatei.dump: string; begin if teilchen<0 then result:=emFeldNamen[emInhalt] else result:=matFeldNamen[matInhalt]+inttostr(teilchen); if ableitung then result:=result+''''; result:=result+' -> '+pre+'-$i'+suf; end; procedure tAusgabeDatei.schreibeKopf; begin if (tCnt<>zeitAnz) and ((nNum<>0) or (tCnt<>-1)) then begin pro.schreibe('Falsche Anzahl an Zeitschritten in '+pre+'-'+inttostr(nNum-1)+suf+' geschrieben ('+inttostr(tCnt)+' statt '+inttostr(zeitAnz)+')!',true); pro.destroyall; halt(1); end; tCnt:=0; assignfile(datei,pre+'-'+inttostr(nNum)+suf); inc(nNum); rewrite(datei,1); blockwrite(datei,nNum,sizeof(integer)); blockwrite(datei,zeitAnz,sizeof(integer)); closefile(datei); end; procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sDX: extended); var i,cnt: longint; sX,cX: extended; begin if (teilchen>=gitter.felders[gitter.aktuelleFelder].matAnz) then begin pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true); pro.destroyall; halt(1); end; if (teilchen<0) and (emInhalt=efA) then gitter.berechne('A',teilchen); if (teilchen>=0) and (matInhalt=mfP) then gitter.berechne('P',teilchen); if gitter.t>=nNum+sDT/2 then schreibeKopf; cnt:=floor((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)/sDX*gitter.dX+1); cX:=gitter.xl; sX:=cX+(cnt-1)*sDX; inc(tCnt); reset(datei,1); seek(datei,fileSize(datei)); blockWrite(datei,cX,sizeof(extended)); blockWrite(datei,sX,sizeof(extended)); blockWrite(datei,cnt,sizeof(integer)); sX:=cX-Min(gitter.dX,sDX)/2; if teilchen<0 then begin for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin while cX>=sX do begin dec(cnt); blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],sizeof(extended)); sX:=sX+sDX; 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); blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],sizeof(extended)); sX:=sX+sDX; end; cX:=cX+gitter.dX; end; end; closeFile(datei); if cnt<>0 then begin pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true); pro.destroyall; halt(1); end; end; { tTeilchenSpezies } constructor tTeilchenSpezies.create; begin inherited create; fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0); setlength(ortsGrenzen,0); fillchar(dichteStuecke,sizeof(dichteStuecke),#0); setlength(dichteStuecke,0); nMin:=1E-3; nMax:=nMin; spezifischeLadung:=0; end; destructor tTeilchenSpezies.destroy; begin setlength(ortsGrenzen,0); setlength(dichteStuecke,0); inherited destroy; end; procedure tTeilchenSpezies.dump(prot: tProtokollant; prefix: string); var i: longint; begin prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc'); prot.schreibe(prefix+'q/m = '+floattostr(spezifischeLadung)+' e/me'); prot.schreibe(prefix+'n(x):'); for i:=0 to length(dichteStuecke)-1 do begin prot.schreibe(prefix+' '+dichteStuecke[i]); if i verteilung stückweise!',true); prot.destroyall; halt(1); end; if length(dichteStuecke)=length(ortsGrenzen) then begin // neues Funktionsstück wird definiert setlength(dichteStuecke,length(dichteStuecke)+1); dichteStuecke[length(dichteStuecke)-1]:=s; continue; end; // kein neues Funktionsstück => if s='stückweiseEnde' then break; // Ende oder // neue Ortsgrenze setlength(ortsGrenzen,length(ortsGrenzen)+1); ortsGrenzen[length(ortsGrenzen)-1]:=exprToFloat(false,s,kvs,nil); until false; end else if startetMit('wie teilchen',rest) then begin ori:=strtoint(rest)-1; if (ori<0) or (ori>=length(teilchen)-1) then begin prot.schreibe('Kann mich nicht auf die Dichte von Teilchen '+inttostr(ori+1)+' beziehen, weil es noch nicht definiert wurde!',true); prot.destroyall; halt(1); end; setlength(ortsGrenzen,length(teilchen[ori].ortsGrenzen)); for i:=0 to length(ortsGrenzen)-1 do ortsGrenzen[i]:=teilchen[ori].ortsGrenzen[i]; setlength(dichteStuecke,length(teilchen[ori].dichteStuecke)); for i:=0 to length(dichteStuecke)-1 do dichteStuecke[i]:=teilchen[ori].dichteStuecke[i]; end else begin setlength(ortsGrenzen,0); setlength(dichteStuecke,1); dichteStuecke[0]:=rest; end; end; { tWertePunkt } constructor tWertePunkt.create(linkerNachbar: tWertePunkt; materien: longint); var emF: tEMFeldInhalt; maF: tMaterieFeldInhalt; abl: boolean; i: longint; begin inherited create; 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; 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 tWertePunkt.berechneGammaUndP; // aus aktuellem dPsi/dX und A var i: longint; begin 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; end; 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]:= -( rN^.rN^.Groessen[false,fiN]*rN^.rN^.Groessen[false,fiiGamma]*rN^.rN^.Groessen[false,fiPX] + rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*rN^.Groessen[false,fiPX] - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*lN^.Groessen[false,fiPX] - lN^.lN^.Groessen[false,fiN]*lN^.lN^.Groessen[false,fiiGamma]*lN^.lN^.Groessen[false,fiPX] )*iDX/6; // es wird über die beiden nächsten linken bzw. rechten Nachbarn gemittelt *) 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 + lN.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 tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); var i: longint; begin if rechts then begin for i:=0 to length(matWerte)-1 do // rechter Randterm von dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n matWerte[i,mfN,true]:= ( lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false]) - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax - matWerte[i,mfPX,false])) * iDX/2; end else begin for i:=0 to length(matWerte)-1 do // linker Randterm von 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]) - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax + matWerte[i,mfPX,false])) * iDX/2; end end; procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen var i: longint; begin // d2Psi/dxdt = dPhi/dx - dGamma/dx 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 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 emWerte[efAX,true]:=emWerte[efDAXDT,false];// + emWerte[efDAXDT,true]*dT; emWerte[efAY,true]:=emWerte[efDAYDT,false];// + emWerte[efDAYDT,true]*dT; emWerte[efAZ,true]:=emWerte[efDAZDT,false];// + emWerte[efDAZDT,true]*dT; end; function tWertePunkt.maxDT: extended; var i: longint; begin result:=infinity; for i:=0 to length(matWerte)-1 do if (matWerte[i,mfN,true]<0) and (matWerte[i,mfN,false]>0) then result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]); end; { linearkombinationsspezifische Methoden von tWertePunkt und tFelder } {$INCLUDE linearkombination.inc} {$DEFINE lkA3} {$INCLUDE linearkombination.inc} {$DEFINE lkA4} {$INCLUDE linearkombination.inc} {$DEFINE lkA5} {$INCLUDE linearkombination.inc} {$DEFINE lkA6} {$INCLUDE linearkombination.inc} {$DEFINE lkA7} {$INCLUDE linearkombination.inc} {$DEFINE lkA8} {$INCLUDE linearkombination.inc} {$DEFINE lkA9} {$INCLUDE linearkombination.inc} {$DEFINE lkA10} {$INCLUDE linearkombination.inc} {$DEFINE lkA11} {$INCLUDE linearkombination.inc} {$DEFINE lkA12} {$INCLUDE linearkombination.inc} {$DEFINE lkA14} {$INCLUDE linearkombination.inc} {$DEFINE lkA15} {$INCLUDE linearkombination.inc} {$DEFINE lkA16} {$INCLUDE linearkombination.inc} {$DEFINE lkA22} {$INCLUDE linearkombination.inc} {$DEFINE lkA23} {$INCLUDE linearkombination.inc} { tFelder } constructor tFelder.create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); var i,j: longint; rechts: boolean; begin inherited create; par:=parent; matAnz:=length(teilchen); fillchar(spezLadungen,sizeof(spezLadungen),#0); setlength(spezLadungen,matAnz); for i:=0 to length(spezLadungen)-1 do spezLadungen[i]:=teilchen[i].spezifischeLadung; 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 for j:=0 to length(teilchen)-1 do begin inhalt[i].matWerte[j,mfN,false]:=teilchen[j].gibDichte(par.xl+(i-2)*par.dX,parent.kvs); inhalt[i].matWerte[j,mfDPsiDX,false]:=0; end; for rechts:=false to true do begin lichters[rechts]:=tMyStringlist.create(nil,''); lichters[rechts].text:=lichter.text; end; lichters[false].grep('^links .*'); lichters[false].subst('^links *',''); lichters[true].grep('^rechts .*'); lichters[true].subst('^rechts *',''); end; destructor tFelder.destroy; var i: longint; begin setlength(spezLadungen,0); for i:=0 to length(inhalt)-1 do inhalt[i].free; setlength(inhalt,0); lichters[false].free; lichters[true].free; inherited destroy; end; procedure tFelder.setzeRaender(dT,iDX,iDT: extended); var emF: tEMFeldInhalt; rechts: boolean; i: longint; nVal: extended; 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) for rechts:=false to true do for i:=0 to lichters[rechts].count-1 do begin par.kvs.add('t',par.t+dT); nVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil); par.kvs.add('t',par.t); inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true]:= inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true] + (nVal - exprToFloat(false,lichters[rechts][i],par.kvs,nil))*iDT; end; end; procedure tFelder.berechneAbleitungen(dT,dX,iDT,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; for i:=2 to length(inhalt)-3 do inhalt[i].berechneNAbleitung(iDX,pDNMax); inhalt[1].berechneNAbleitung(iDX,pDNMax,false); inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true); setzeRaender(dT,iDX,iDT); for i:=1 to length(inhalt)-2 do inhalt[i].berechneAbleitungen(dX,iDX); end; function tFelder.maxDT: extended; var i: longint; begin result:=inhalt[0].maxDT; for i:=1 to length(inhalt)-1 do result:=min(result,inhalt[i].maxDT); end; { tGitter } constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); var i: longint; begin inherited create; zeitverfahren:=ZV; kvs:=bekannteWerte; prot:=tProtokollant.create(protokollant,name); dTMaximum:=deltaT; dX:=deltaX; iDX:=1/dX; if pDNMa < 0 then begin pDNMax:=0; pDNMaxDynamisch:=true; end else pDNMax:=pDNMa; case Zeitverfahren of zfEulerVorwaerts: Setlength(Felders,2); zfRungeKuttaDreiAchtel,zfRungeKuttaVier: Setlength(Felders,5); zfRungeKuttaZehn: Setlength(Felders,18); zfRungeKuttaZwoelf: Setlength(Felders,26); zfRungeKuttaVierzehn: Setlength(Felders,36); end{of Case}; xl:=dX/2; for i:=0 to length(felders)-1 do felders[i]:=tFelder.create(size,teilchen,lichter,self); aktuelleFelder:=0; t:=0; kvs.add('t',t); end; destructor tGitter.destroy; var i: longint; begin for i:=0 to length(Felders)-1 do felders[i].free; setlength(felders,0); inherited destroy; 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; function tGitter.pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert? begin mDT:=min(mDT,felders[aF].maxDT); // das maximale dT, welches bei den Momentanen Ableitungen nicht zu // unphysikalischen Effekten (z.B. negativen Dichten) führt result:=mDT '+floattostr(dT)); if dT<1E-30 then begin prot.schreibe('pruefeMaxDT: Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < 10^-30) - irgendwas scheint grundsätzlich kaputt zu sein!',true); prot.destroyall; halt(1); end; {$ENDIF} exit; end; end; procedure tGitter.iteriereSchritt(var dT: extended); var i: longint; {$IFDEF Zeitschrittueberwachung} j: longint; {$ENDIF} {$IFDEF Dichteueberwachung} pro: tProtokollant; {$ENDIF} iDT,mDT: extended; begin kvs.add('t',t); diffusionsTermAnpassen(prot); mDT:=infinity; repeat iDT:=1/dT; felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t)) if pruefeMaxDT(aktuelleFelder,mDT,dT) then continue; case zeitverfahren of zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT); zfRungeKuttaDreiAchtel: begin felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/3); // ya = y(t) + y' dt/3 felders[2].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/3,ya) if pruefeMaxDT(2,mDT,dT) then continue; felders[3].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],-dT/3,dT); // yb = y(t) - y' dt/3 + ya' dt felders[3].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // yb' = y'(t+2dt/3,yb) if pruefeMaxDT(3,mDT,dT) then continue; felders[4].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],felders[3],dT,-dT,dT); // yc = y(t) + a' dt - ya' dt + yb' dt felders[4].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc) if pruefeMaxDT(4,mDT,dT) then continue; felders[1-aktuelleFelder].liKo( felders[aktuelleFelder], felders[aktuelleFelder], felders[2], felders[3], felders[4], dT/8, dT*3/8, dT*3/8, dT/8 ); // y(t+dt) = y(t) + (y' + 3(ya' + yb') + yc') dt/8 end; zfRungeKuttaVier: begin felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2 felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/2,ya) if pruefeMaxDT(2,mDT,dT) then continue; felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2 felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yb' = y'(t+dt/2,yb) if pruefeMaxDT(3,mDT,dT) then continue; felders[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc) if pruefeMaxDT(4,mDT,dT) then continue; 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; zfRungeKuttaZehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk108.txt felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT*0.1); felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(2,mDT,dT) then continue; felders[3].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2], -0.915176561375291440520015019275342154318951387664369720564660 * dT, 1.45453440217827322805250021715664459117622483736537873607016 * dT); felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(3,mDT,dT) then continue; felders[4].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3], 0.202259190301118170324681949205488413821477543637878380814562 * dT, 0.606777570903354510974045847616465241464432630913635142443687 * dT); felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(4,mDT,dT) then continue; felders[5].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3],felders[4], 0.184024714708643575149100693471120664216774047979591417844635 * dT, 0.197966831227192369068141770510388793370637287463360401555746 * dT, -0.0729547847313632629185146671595558023015011608914382961421311* dT); felders[5].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(5,mDT,dT) then continue; felders[6].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[4],felders[5], 0.0879007340206681337319777094132125475918886824944548534041378* dT, 0.410459702520260645318174895920453426088035325902848695210406 * dT, 0.482713753678866489204726942976896106809132737721421333413261 * dT); felders[6].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(6,mDT,dT) then continue; felders[7].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[4],felders[5],felders[6], 0.0859700504902460302188480225945808401411132615636600222593880* dT, 0.330885963040722183948884057658753173648240154838402033448632 * dT, 0.489662957309450192844507011135898201178015478433790097210790 * dT, -0.0731856375070850736789057580558988816340355615025188195854775* dT); felders[7].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(7,mDT,dT) then continue; felders[8].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[5],felders[6],felders[7], 0.120930449125333720660378854927668953958938996999703678812621 * dT, 0.260124675758295622809007617838335174368108756484693361887839 * dT, 0.0325402621549091330158899334391231259332716675992700000776101* dT, -0.0595780211817361001560122202563305121444953672762930724538856* dT); felders[8].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(8,mDT,dT) then continue; felders[9].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8], 0.110854379580391483508936171010218441909425780168656559807038 * dT, -0.0605761488255005587620924953655516875526344415354339234619466* dT, 0.321763705601778390100898799049878904081404368603077129251110 * dT, 0.510485725608063031577759012285123416744672137031752354067590 * dT); felders[9].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(9,mDT,dT) then continue; felders[10].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9], 0.112054414752879004829715002761802363003717611158172229329393 * dT, -0.144942775902865915672349828340980777181668499748506838876185 * dT, -0.333269719096256706589705211415746871709467423992115497968724 * dT, 0.499269229556880061353316843969978567860276816592673201240332 * dT, 0.509504608929686104236098690045386253986643232352989602185060 * dT); felders[10].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(10,mDT,dT) then continue; felders[11].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9],felders[10], 0.113976783964185986138004186736901163890724752541486831640341 * dT, -0.0768813364203356938586214289120895270821349023390922987406384* dT, 0.239527360324390649107711455271882373019741311201004119339563 * dT, 0.397774662368094639047830462488952104564716416343454639902613 * dT, 0.0107558956873607455550609147441477450257136782823280838547024* dT, -0.327769124164018874147061087350233395378262992392394071906457 * dT); felders[11].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(11,mDT,dT) then continue; felders[12].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9],felders[10],felders[11], 0.0798314528280196046351426864486400322758737630423413945356284* dT, -0.0520329686800603076514949887612959068721311443881683526937298* dT, -0.0576954146168548881732784355283433509066159287152968723021864* dT, 0.194781915712104164976306262147382871156142921354409364738090 * dT, 0.145384923188325069727524825977071194859203467568236523866582 * dT, -0.0782942710351670777553986729725692447252077047239160551335016* dT, -0.114503299361098912184303164290554670970133218405658122674674 * dT); felders[12].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(12,mDT,dT) then continue; felders[13].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[4],felders[5],felders[6],felders[7],felders[8],felders[9],felders[10],felders[11],felders[12], 0.985115610164857280120041500306517278413646677314195559520529 * dT, 0.330885963040722183948884057658753173648240154838402033448632 * dT, 0.489662957309450192844507011135898201178015478433790097210790 * dT, -1.37896486574843567582112720930751902353904327148559471526397 * dT, -0.861164195027635666673916999665534573351026060987427093314412 * dT, 5.78428813637537220022999785486578436006872789689499172601856 * dT, 3.28807761985103566890460615937314805477268252903342356581925 * dT, -2.38633905093136384013422325215527866148401465975954104585807 * dT, -3.25479342483643918654589367587788726747711504674780680269911 * dT, -2.16343541686422982353954211300054820889678036420109999154887 * dT); felders[13].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(13,mDT,dT) then continue; felders[14].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3],felders[4],felders[6],felders[7],felders[8],felders[9],felders[10],felders[11],felders[12],felders[13], 0.895080295771632891049613132336585138148156279241561345991710 * dT, 0.197966831227192369068141770510388793370637287463360401555746 * dT, -0.0729547847313632629185146671595558023015011608914382961421311* dT, -0.851236239662007619739049371445966793289359722875702227166105 * dT, 0.398320112318533301719718614174373643336480918103773904231856 * dT, 3.63937263181035606029412920047090044132027387893977804176229 * dT, 1.54822877039830322365301663075174564919981736348973496313065 * dT, -2.12221714704053716026062427460427261025318461146260124401561 * dT, -1.58350398545326172713384349625753212757269188934434237975291 * dT, -1.71561608285936264922031819751349098912615880827551992973034 * dT, -0.0244036405750127452135415444412216875465593598370910566069132* dT); felders[14].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(14,mDT,dT) then continue; felders[15].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],felders[5],felders[7],felders[13],felders[14], -0.915176561375291440520015019275342154318951387664369720564660 * dT, 1.45453440217827322805250021715664459117622483736537873607016 * dT, -0.777333643644968233538931228575302137803351053629547286334469 * dT, -0.0910895662155176069593203555807484200111889091770101799647985* dT, 0.0910895662155176069593203555807484200111889091770101799647985* dT, 0.777333643644968233538931228575302137803351053629547286334469 * dT); felders[15].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(15,mDT,dT) then continue; felders[16].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3],felders[15], 0.1 * dT, -0.157178665799771163367058998273128921867183754126709419409654 * dT, 0.157178665799771163367058998273128921867183754126709419409654 * dT); felders[16].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(16,mDT,dT) then continue; felders[17].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],felders[3],felders[5],felders[6],felders[7],felders[8],felders[9], felders[10],felders[11],felders[12],felders[13],felders[14],felders[15],felders[16], 0.181781300700095283888472062582262379650443831463199521664945 * dT, 0.675 * dT, 0.342758159847189839942220553413850871742338734703958919937260 * dT, 0.259111214548322744512977076191767379267783684543182428778156 * dT, -0.358278966717952089048961276721979397739750634673268802484271 * dT, -1.04594895940883306095050068756409905131588123172378489286080 * dT, 0.930327845415626983292300564432428777137601651182965794680397 * dT, 1.77950959431708102446142106794824453926275743243327790536000 * dT, 0.1 * dT, -0.282547569539044081612477785222287276408489375976211189952877 * dT, -0.159327350119972549169261984373485859278031542127551931461821 * dT, -0.145515894647001510860991961081084111308650130578626404945571 * dT, -0.259111214548322744512977076191767379267783684543182428778156 * dT, -0.342758159847189839942220553413850871742338734703958919937260 * dT, -0.675 * dT); felders[17].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(17,mDT,dT) then continue; felders[1-aktuelleFelder].liKo( felders[aktuelleFelder], felders[aktuelleFelder], felders[2], felders[3], felders[5], felders[7], felders[9], felders[10], felders[11], felders[12], felders[13], felders[14], felders[15], felders[16], felders[17], dT/30, dT/40, dT/30, dT/20, dT/25, 0.189237478148923490158306404106012326238162346948625830327194 * dT, 0.277429188517743176508360262560654340428504319718040836339472 * dT, 0.277429188517743176508360262560654340428504319718040836339472 * dT, 0.189237478148923490158306404106012326238162346948625830327194 * dT, -dT/25, -dT/20, -dT/30, -dT/40, -dT/30 ); end; zfRungeKuttaZwoelf: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1210.txt felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT*0.2); felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(2,mDT,dT) then continue; felders[3].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2], -0.216049382716049382716049382716049382716049382716049382716049 * dT, 0.771604938271604938271604938271604938271604938271604938271605 * dT); felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(3,mDT,dT) then continue; felders[4].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3], dT/4.8, 0.625 * dT); felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(4,mDT,dT) then continue; felders[5].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3],felders[4], 0.193333333333333333333333333333333333333333333333333333333333 * dT, 0.22 * dT, -0.08 * dT); felders[5].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(5,mDT,dT) then continue; felders[6].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[4],felders[5], 0.1 * dT, 0.4 * dT, 0.5 * dT); felders[6].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(6,mDT,dT) then continue; felders[7].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[4],felders[5],felders[6], 0.103364471650010477570395435690481791543342708330349879244197 * dT, 0.124053094528946761061581889237115328211074784955180298044074 * dT, 0.483171167561032899288836480451962508724109257517289177302380 * dT, -0.0387530245694763252085681443767620580395733302341368038804290* dT); felders[7].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(7,mDT,dT) then continue; felders[8].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[5],felders[6],felders[7], 0.124038261431833324081904585980175168140024670698633612292480 * dT, 0.217050632197958486317846256953159942875916353757734167684657 * dT, 0.0137455792075966759812907801835048190594443990939408530842918* dT, -0.0661095317267682844455831341498149531672668252085016565917546* dT); felders[8].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(8,mDT,dT) then continue; felders[9].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8], 0.0914774894856882983144991846980432197088832099976660100090486 * dT, -0.00544348523717469689965754944144838611346156873847009178068318* dT, 0.0680716801688453518578515120895103863112751730758794372203952 * dT, 0.408394315582641046727306852653894780093303185664924644551239 * dT); felders[9].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(9,mDT,dT) then continue; felders[10].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9], 0.0890013652502551018954509355423841780143232697403434118692699 * dT, 0.00499528226645532360197793408420692800405891149406814091955810* dT, 0.397918238819828997341739603001347156083435060931424970826304 * dT, 0.427930210752576611068192608300897981558240730580396406312359 * dT, -0.0865117637557827005740277475955029103267246394128995965941585 * dT); felders[10].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(10,mDT,dT) then continue; felders[11].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9],felders[10], 0.0695087624134907543112693906409809822706021061685544615255758 * dT, 0.129146941900176461970759579482746551122871751501482634045487 * dT, 1.53073638102311295076342566143214939031177504112433874313011 * dT, 0.577874761129140052546751349454576715334892100418571882718036 * dT, -0.951294772321088980532340837388859453930924498799228648050949 * dT, -0.408276642965631951497484981519757463459627174520978426909934 * dT); felders[11].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(11,mDT,dT) then continue; felders[12].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9],felders[10],felders[11], 0.0444861403295135866269453507092463581620165501018684152933313 * dT, -0.00380476867056961731984232686574547203016331563626856065717964 * dT, 0.0106955064029624200721262602809059154469206077644957399593972 * dT, 0.0209616244499904333296674205928919920806734650660039898074652 * dT, -0.0233146023259321786648561431551978077665337818756053603898847 * dT, 0.00263265981064536974369934736325334761174975280887405725010964 * dT, 0.00315472768977025060103545855572111407955208306374459723959783 * dT); felders[12].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(12,mDT,dT) then continue; felders[13].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[9],felders[10],felders[11],felders[12], 0.0194588815119755475588801096525317761242073762016273186231215 * dT, 0.0000678512949171812509306121653452367476194364781259165332321534 * dT, -0.0000429795859049273623271005330230162343568863387724883603675550 * dT, 0.0000176358982260285155407485928953302139937553442829975734148981 * dT, 0.0653866627415027051009595231385181033549511358787382098351924 * dT); felders[13].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(13,mDT,dT) then continue; felders[14].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[9],felders[10],felders[11],felders[12],felders[13], 0.206836835664277105916828174798272361078909196043446411598231 * dT, 0.0166796067104156472828045866664696450306326505094792505215514 * dT, -0.00879501563200710214457024178249986591130234990219959208704979 * dT, 0.00346675455362463910824462315246379209427513654098596403637231 * dT, -0.861264460105717678161432562258351242030270498966891201799225 * dT, 0.908651882074050281096239478469262145034957129939256789178785 * dT); felders[14].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(14,mDT,dT) then continue; felders[15].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[9],felders[10],felders[11],felders[12],felders[13],felders[14], 0.0203926084654484010091511314676925686038504449562413004562382 * dT, 0.0869469392016685948675400555583947505833954460930940959577347 * dT, -0.0191649630410149842286436611791405053287170076602337673587681 * dT, 0.00655629159493663287364871573244244516034828755253746024098838 * dT, 0.0987476128127434780903798528674033899738924968006632201445462 * dT, 0.00535364695524996055083260173615567408717110247274021056118319 * dT, 0.301167864010967916837091303817051676920059229784957479998077 * dT); felders[15].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(15,mDT,dT) then continue; felders[16].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[9],felders[10],felders[11],felders[12],felders[13],felders[14],felders[15], 0.228410433917778099547115412893004398779136994596948545722283 * dT, -0.498707400793025250635016567442511512138603770959682292383042 * dT, 0.134841168335724478552596703792570104791700727205981058201689 * dT, -0.0387458244055834158439904226924029230935161059142806805674360 * dT, -1.27473257473474844240388430824908952380979292713250350199641 * dT, 1.43916364462877165201184452437038081875299303577911839630524 * dT, -0.214007467967990254219503540827349569639028092344812795499026 * dT, 0.958202417754430239892724139109781371059908874605153648768037 * dT); felders[16].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(16,mDT,dT) then continue; felders[17].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[9],felders[10],felders[11],felders[12],felders[13],felders[14],felders[15],felders[16], 2.00222477655974203614249646012506747121440306225711721209798 * dT, 2.06701809961524912091954656438138595825411859673341600679555 * dT, 0.623978136086139541957471279831494466155292316167021080663140 * dT, -0.0462283685500311430283203554129062069391947101880112723185773 * dT, -8.84973288362649614860075246727118949286604835457092701094630 * dT, 7.74257707850855976227437225791835589560188590785037197433615 * dT, -0.588358519250869210993353314127711745644125882130941202896436 * dT, -1.10683733362380649395704708016953056176195769617014899442903 * dT, -0.929529037579203999778397238291233214220788057511899747507074 * dT); felders[17].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(17,mDT,dT) then continue; felders[18].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9], felders[10],felders[11],felders[12],felders[13],felders[14],felders[15],felders[16],felders[17], 3.13789533412073442934451608989888796808161259330322100268310 * dT, 0.129146941900176461970759579482746551122871751501482634045487 * dT, 1.53073638102311295076342566143214939031177504112433874313011 * dT, 0.577874761129140052546751349454576715334892100418571882718036 * dT, 5.42088263055126683050056840891857421941300558851862156403363 * dT, 0.231546926034829304872663800877643660904880180835945693836936 * dT, 0.0759292995578913560162301311785251873561801342333194895292058 * dT, -12.3729973380186513287414553402595806591349822617535905976253 * dT, 9.85455883464769543935957209317369202080367765721777101906955 * dT, 0.0859111431370436529579357709052367772889980495122329601159540 * dT, -5.65242752862643921117182090081762761180392602644189218673969 * dT, -1.94300935242819610883833776782364287728724899124166920477873 * dT, -0.128352601849404542018428714319344620742146491335612353559923 * dT); felders[18].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(18,mDT,dT) then continue; felders[19].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[6],felders[7],felders[8],felders[9], felders[10],felders[11],felders[12],felders[13],felders[14],felders[15],felders[16],felders[17],felders[18], 1.38360054432196014878538118298167716825163268489922519995564 * dT, 0.00499528226645532360197793408420692800405891149406814091955810 * dT, 0.397918238819828997341739603001347156083435060931424970826304 * dT, 0.427930210752576611068192608300897981558240730580396406312359 * dT, -1.30299107424475770916551439123047573342071475998399645982146 * dT, 0.661292278669377029097112528107513072734573412294008071500699 * dT, -0.144559774306954349765969393688703463900585822441545655530145 * dT, -6.96576034731798203467853867461083919356792248105919255460819 * dT, 6.65808543235991748353408295542210450632193197576935120716437 * dT, -1.66997375108841486404695805725510845049807969199236227575796 * dT, 2.06413702318035263832289040301832647130604651223986452170089 * dT, -0.674743962644306471862958129570837723192079875998405058648892 * dT, -0.00115618834794939500490703608435907610059605754935305582045729 * dT, -0.00544057908677007389319819914241631024660726585015012485938593 * dT); felders[19].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(19,mDT,dT) then continue; felders[20].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[5],felders[6],felders[7],felders[9], felders[10],felders[11],felders[12],felders[13],felders[14],felders[15],felders[16],felders[17],felders[18],felders[19], 0.951236297048287669474637975894973552166903378983475425758226 * dT, 0.217050632197958486317846256953159942875916353757734167684657 * dT, 0.0137455792075966759812907801835048190594443990939408530842918 * dT, -0.0661095317267682844455831341498149531672668252085016565917546 * dT, 0.152281696736414447136604697040747131921486432699422112099617 * dT, -0.337741018357599840802300793133998004354643424457539667670080 * dT, -0.0192825981633995781534949199286824400469353110630787982121133 * dT, -3.68259269696866809932409015535499603576312120746888880201882 * dT, 3.16197870406982063541533528419683854018352080342887002331312 * dT, -0.370462522106885290716991856022051125477943482284080569177386 * dT, -0.0514974200365440434996434456698127984941168616474316871020314 * dT, -0.000829625532120152946787043541792848416659382675202720677536554 * dT, 0.00000279801041419278598986586589070027583961355402640879503213503 * dT, 0.0418603916412360287969841020776788461794119440689356178942252 * dT, 0.279084255090877355915660874555379649966282167560126269290222 * dT); felders[20].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(20,mDT,dT) then continue; felders[21].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[4],felders[5],felders[6],felders[8], felders[10],felders[11],felders[18],felders[19],felders[20], 0.103364471650010477570395435690481791543342708330349879244197 * dT, 0.124053094528946761061581889237115328211074784955180298044074 * dT, 0.483171167561032899288836480451962508724109257517289177302380 * dT, -0.0387530245694763252085681443767620580395733302341368038804290 * dT, -0.438313820361122420391059788940960176420682836652600698580091 * dT, -0.218636633721676647685111485017151199362509373698288330593486 * dT, -0.0312334764394719229981634995206440349766174759626578122323015 * dT, 0.0312334764394719229981634995206440349766174759626578122323015 * dT, 0.218636633721676647685111485017151199362509373698288330593486 * dT, 0.438313820361122420391059788940960176420682836652600698580091 * dT); felders[21].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(21,mDT,dT) then continue; felders[22].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3],felders[4],felders[7],felders[8], felders[10],felders[11],felders[18],felders[19],felders[20],felders[21], 0.193333333333333333333333333333333333333333333333333333333333 * dT, 0.22 * dT, -0.08 * dT, 0.0984256130499315928152900286856048243348202521491288575952143 * dT, -0.196410889223054653446526504390100417677539095340135532418849 * dT, 0.436457930493068729391826122587949137609670676712525034763317 * dT, 0.0652613721675721098560370939805555698350543810708414716730270 * dT, -0.0652613721675721098560370939805555698350543810708414716730270 * dT, -0.436457930493068729391826122587949137609670676712525034763317 * dT, 0.196410889223054653446526504390100417677539095340135532418849 * dT, -0.0984256130499315928152900286856048243348202521491288575952143 * dT); felders[22].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(22,mDT,dT) then continue; felders[23].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],felders[5],felders[7],felders[21], felders[22], -0.216049382716049382716049382716049382716049382716049382716049 * dT, 0.771604938271604938271604938271604938271604938271604938271605 * dT, -0.666666666666666666666666666666666666666666666666666666666667 * dT, -0.390696469295978451446999802258495981249099665294395945559163 * dT, 0.390696469295978451446999802258495981249099665294395945559163 * dT, 0.666666666666666666666666666666666666666666666666666666666667 * dT); felders[23].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(23,mDT,dT) then continue; felders[24].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[3],felders[23], 0.2 * dT, -0.164609053497942386831275720164609053497942386831275720164609 * dT, 0.164609053497942386831275720164609053497942386831275720164609 * dT); felders[24].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(24,mDT,dT) then continue; felders[25].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],felders[3],felders[5],felders[7], felders[8],felders[9],felders[10],felders[11],felders[12],felders[13],felders[14],felders[15],felders[16], felders[17],felders[18],felders[19],felders[20],felders[21],felders[22],felders[23],felders[24], 1.47178724881110408452949550989023611293535315518571691939396 * dT, 0.7875 * dT, 0.421296296296296296296296296296296296296296296296296296296296 * dT, 0.291666666666666666666666666666666666666666666666666666666667 * dT, 0.348600717628329563206854421629657569274689947367847465753757 * dT, 0.229499544768994849582890233710555447073823569666506700662510 * dT, 5.79046485790481979159831978177003471098279506036722411333192 * dT, 0.418587511856506868874073759426596207226461447604248151080016 * dT, 0.307039880222474002649653817490106690389251482313213999386651 * dT, -4.68700905350603332214256344683853248065574415794742040470287 * dT, 3.13571665593802262152038152399873856554395436199962915429076 * dT, 1.40134829710965720817510506275620441055845017313930508348898 * dT, -5.52931101439499023629010306005764336421276055777658156400910 * dT, -0.853138235508063349309546894974784906188927508039552519557498 * dT, 0.103575780373610140411804607167772795518293914458500175573749 * dT, -0.140474416950600941142546901202132534870665923700034957196546 * dT, -0.418587511856506868874073759426596207226461447604248151080016 * dT, -0.229499544768994849582890233710555447073823569666506700662510 * dT, -0.348600717628329563206854421629657569274689947367847465753757 * dT, -0.291666666666666666666666666666666666666666666666666666666667 * dT, -0.421296296296296296296296296296296296296296296296296296296296 * dT, -0.7875 * dT); felders[25].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); if pruefeMaxDT(25,mDT,dT) then continue; felders[1-aktuelleFelder].liKo( felders[aktuelleFelder], felders[aktuelleFelder], felders[2], felders[3], felders[5], felders[7], felders[8], felders[10], felders[11], felders[13], felders[14], felders[15], felders[16], felders[17], felders[18], felders[19], felders[20], felders[21], felders[22], felders[23], felders[24], felders[25], 0.0238095238095238095238095238095238095238095238095238095238095 * dT, 0.0234375 * dT, 0.03125 * dT, 0.0416666666666666666666666666666666666666666666666666666666667 * dT, 0.05 * dT, 0.05 * dT, 0.1 * dT, 0.0714285714285714285714285714285714285714285714285714285714286 * dT, 0.138413023680782974005350203145033146748813640089941234591267 * dT, 0.215872690604931311708935511140681138965472074195773051123019 * dT, 0.243809523809523809523809523809523809523809523809523809523810 * dT, 0.215872690604931311708935511140681138965472074195773051123019 * dT, 0.138413023680782974005350203145033146748813640089941234591267 * dT, -0.0714285714285714285714285714285714285714285714285714285714286 * dT, -0.1 * dT, -0.05 * dT, -0.05 * dT, -0.0416666666666666666666666666666666666666666666666666666666667 * dT, -0.03125 * dT, -0.0234375 * dT, 0.0238095238095238095238095238095238095238095238095238095238095 * dT ); end; zfRungeKuttaVierzehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1412.txt end; end{of case}; break; until false; aktuelleFelder:=1-aktuelleFelder; {$IFDEF Dichteueberwachung} for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do for j:=0 to felders[aktuelleFelder].matAnz-1 do if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then with felders[1-aktuelleFelder] do begin pro:=tProtokollant.create(prot,'iteriereSchritt'); pro.schreibe('n<0 bei:',true); pro.schreibe(' t = ' +floattostr(t)+',',true); pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true); pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* inhalt[i-1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true); pro.schreibe(' v = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]* inhalt[i].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true); pro.schreibe(' n'' = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true); pro.schreibe(' v_r = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]* inhalt[i+1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); pro.destroyall; halt(1); end; {$ENDIF} t:=t+dT; dT:=max(dT,min(1,mDT)/100); end; function tGitter.dumpErhaltungsgroessen: boolean; var i,j: integer; ns: tExtendedArray; pro: tProtokollant; s: string; begin pro:=tProtokollant.create(prot,'dumpErhaltungsgroessen'); setlength(ns,felders[aktuelleFelder].matAnz); result:=true; s:=''; 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 ns[i]:=ns[i]*dX; s:=s+ ' n['+inttostr(i)+']='+floattostr(ns[i]); {$IFDEF Dichteueberwachung} if ns[i]>1000 then begin result:=false; pro.schreibe(' n['+inttostr(i)+'] > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true); end; {$ENDIF} end; delete(s,1,1); prot.schreibe(s); pro.free; end; procedure tGitter.berechne(was: char; teilchen: longint); var i: longint; tmp: extended; emF: tEMFeldInhalt; maF: tMaterieFeldInhalt; begin if was='A' then begin for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin tmp:=0; for emF:=efAX to efAZ do tmp:=tmp+sqr(felders[aktuelleFelder].inhalt[i].emWerte[emF,false]); felders[aktuelleFelder].inhalt[i].emWerte[efA,false]:=sqrt(tmp); end; exit; end; if was='P' then begin 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[teilchen,maF,false]); felders[aktuelleFelder].inhalt[i].matWerte[teilchen,maF,false]:=sqrt(tmp); end; exit; end; prot.schreibe('Kann '''+was+''' nicht berechnen, weil ich es nicht verstehe!',true); prot.destroyall; halt; end; { tSimulation } constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string); var ifile: tMyStringlist; zeitverfahren: tZeitverfahren; s,t,aSuffix,aPrefix: string; deltaX,breite,pDNMax: extended; i: longint; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; pro: tProtokollant; na: pSigActionRec; begin inherited create; prot:=tProtokollant.create(protokollant,name); kvs:=tKnownValues.create; ifile:=tMyStringlist.create(prot,'ifile'); ifile.loadfromfile(inName); if not ifile.unfoldMacros then begin prot.schreibe('Fehlerhafte Macros in Parameterdatei '''+inName+'''!',true); prot.destroyall; halt(1); end; // Standardeinstellungen Bereich 'allgemein' zeitverfahren:=zfRungeKuttaVier; deltaX:=1e-2; kvs.add('λ',1); kvs.add('T',1); kvs.add('ω',2*pi); kvs.add('dX',deltaX); kvs.add('q',1); kvs.add('me',1); dT:=-1; sDT:=-1; sDX:=-1; tEnde:=100; breite:=10.0; pDNMax:=-1; fortschrittsAnzeige:=false; gotSighup:=false; // Standardeinstellungen Bereich 'ausgaben' aPrefix:=extractfilepath(paramstr(0)); aSuffix:='.dat'; setlength(ausgabeDateien,0); // Standardeinstellungen Breich 'lichtVon...' lichter:=tMyStringlist.create(prot,'lichter'); // Standardeinstellungen Breich 'teilchen...' setlength(teilchen,0); repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in '''+inName+'''!',true); prot.destroyall; halt(1); end; if s='Dateiende' then break; if s='allgemein' then begin repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich allgemein!',true); prot.destroyall; halt(1); end; if s='allgemeinEnde' then break; if s='runge-Kutta-3/8' then begin Zeitverfahren:=zfRungeKuttaDreiAchtel; continue; end; if s='runge-Kutta-4' then begin Zeitverfahren:=zfRungeKuttaVier; continue; end; if s='runge-Kutta-10' then begin Zeitverfahren:=zfRungeKuttaZehn; continue; end; if s='runge-Kutta-12' then begin Zeitverfahren:=zfRungeKuttaZwoelf; continue; end; if s='runge-Kutta-14' then begin Zeitverfahren:=zfRungeKuttaVierzehn; continue; end; if s='euler-Vorwärts' then begin Zeitverfahren:=zfEulerVorwaerts; continue; end; if startetMit('ortsschritt ',s) then begin deltaX:=exprtofloat(false,s,kvs,nil); kvs.add('dX',deltaX); continue; end; if startetMit('zeitschritt ',s) then begin dT:=exprtofloat(false,s,kvs,nil); kvs.add('dT',dT); continue; end; if startetMit('diffusionsterm ',s) then begin pDNMax:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('zeit ',s) then begin 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); prot.destroyall; halt(1); until false; continue; end; if s='ausgaben' then begin repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich ausgaben!',true); prot.destroyall; halt(1); end; if s='ausgabenEnde' then break; if startetMit('suffix ',s) then begin aSuffix:=s; continue; end; if startetMit('prefix ',s) then begin aPrefix:=s; continue; end; if startetMit('zeitschritt ',s) then 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,','); setlength(ausgabeDateien,length(ausgabeDateien)+1); ausgabeDateien[length(ausgabeDateien)-1]:=tAusgabeDatei.create(t,aPrefix,aSuffix,prot,'ausgabeDateien['+inttostr(length(ausgabeDateien)-1)+']'); end; continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true); prot.destroyall; halt(1); until false; continue; end; if startetMit('licht von ',s) then begin if startetMit('links ',s) then begin lichter.add('links '+s); continue; end; if startetMit('rechts ',s) then begin lichter.add('rechts '+s); continue; end; prot.schreibe('Licht kann nicht von '''+erstesArgument(s)+''' kommen!',true); prot.destroyall; halt(1); end; if startetMit('teilchen',s) then begin if (s<>'') and (strtoint(s)<>length(teilchen)+1) then begin prot.schreibe('Ich erwarte die Teilchen beginnend bei 1 der Reihe nach in Parameterdatei '''+inName+'''!',true); prot.destroyall; halt(1); end; setlength(teilchen,length(teilchen)+1); teilchen[length(teilchen)-1]:=tTeilchenSpezies.create; repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich teilchen!',true); prot.destroyall; halt(1); end; if s='teilchen'+inttostr(length(teilchen))+'Ende' then break; if startetMit('spezifische Ladung ',s) then begin teilchen[length(teilchen)-1].spezifischeLadung:=exprToFloat(false,s,kvs,nil); kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung); continue; end; if startetMit('maximaldichte ',s) then begin teilchen[length(teilchen)-1].nMax:=exprToFloat(false,s,kvs,nil); kvs.add('nmax'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMax); continue; end; if startetMit('minimaldichte ',s) then begin teilchen[length(teilchen)-1].nMin:=exprToFloat(false,s,kvs,nil); kvs.add('nmin'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMin); continue; end; if startetMit('verteilung ',s) then begin teilchen[length(teilchen)-1].liesDichteFunktion(s,ifile,kvs,teilchen,prot); continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich teilchen!',true); prot.destroyall; halt(1); until false; continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true); prot.destroyall; halt(1); until false; ifile.free; if length(ausgabedateien)=0 then begin prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true); prot.destroyall; halt(1); end; if dT<0 then dT:=deltaX/10; if sDT<0 then sDT:=dT; sDT:=1/round(1/sDT); sT:=sDT/2; if sDX<0 then sDX:=deltaX; pro:=tProtokollant.create(prot,'create'); case zeitverfahren of zfEulerVorwaerts: pro.schreibe('Iteration mittels Euler-Vorwärts'); zfRungeKuttaDreiAchtel: pro.schreibe('Iteration mittels Runge-Kutta-3/8'); zfRungeKuttaVier: pro.schreibe('Iteration mittels Runge-Kutta-4'); zfRungeKuttaZehn: pro.schreibe('Iteration mittels Runge-Kutta-10'); zfRungeKuttaZwoelf: pro.schreibe('Iteration mittels Runge-Kutta-12'); zfRungeKuttaVierzehn: pro.schreibe('Iteration mittels Runge-Kutta-14'); else pro.schreibe('Iteration mittels unbekanntem Verfahren'); end{of case}; pro.schreibe('dX = '+floattostr(deltaX)); pro.schreibe('dT = '+floattostr(dT)); pro.schreibe('breite = '+floattostr(breite)+' (= '+inttostr(round(breite/deltaX)+1)+' Punkte)'); pro.schreibe('pDNMax = '+floattostr(pDNMax)); pro.schreibe('bekannte Werte:'); kvs.dump(pro,' '); if length(teilchen)>0 then begin pro.schreibe('teilchen:'); for i:=0 to length(teilchen)-1 do teilchen[i].dump(pro,' '); end; if lichter.count>0 then begin pro.schreibe('lichter:'); lichter.dump(pro,' '); end; if length(ausgabeDateien)>0 then begin pro.schreibe('ausgaben:'); for i:=0 to length(ausgabeDateien)-1 do begin ausgabeDateien[i].zeitAnz:=round(1/sDT); ausgabeDateien[i].sDT:=sDT; pro.schreibe(' '+ausgabeDateien[i].dump); end; end; pro.free; if length(sighupSimulationen)=0 then begin new(na); na^.sa_Handler := sigActionHandler(@signalCapture); fillchar(na^.sa_Mask,sizeof(na^.sa_mask),#0); na^.sa_Flags := 0; {$ifdef Linux} // Linux specific na^.sa_Restorer := Nil; {$endif} if (fPSigaction(SIGHUP, na, nil) <> 0) then begin writeln(stderr, 'Error: ', fpgeterrno, '.'); halt(1); end; dispose(na); end; setlength(sighupSimulationen,length(sighupSimulationen)+1); sighupSimulationen[length(sighupSimulationen)-1]:=self; gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter'); for i:=0 to length(teilchen)-1 do teilchen[i].free; setlength(teilchen,0); lichter.free; end; destructor tSimulation.destroy; var i,j: longint; begin gitter.free; prot.free; for i:=length(sighupSimulationen)-1 downto 0 do if sighupSimulationen[i]=self then begin for j:=i+1 to length(sighupSimulationen)-1 do sighupSimulationen[j-1]:=sighupSimulationen[j]; setlength(sighupSimulationen,length(sighupSimulationen)-1); end; inherited destroy; end; function tSimulation.iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende? var i: longint; begin result:=false; zeitPhysik:=zeitPhysik-now; if errorCode<2 then gitter.iteriereSchritt(dT); zeitPhysik:=zeitPhysik+now; zeitDatei:=zeitDatei-now; while (gitter.t>=sT) and (sT