unit Physikunit; {$mode objfpc}{$H+} {$IFNDEF UNIX} {$ERROR This program can be compiled only on/for Unix/Linux based systems.} {$ENDIF} {$MACRO on} // nach einer LiKo negative Dichten entfernen: { $DEFINE DichteNichtnegativieren:=perSkalieren} // insgesamt runter skalieren um negative Dichten aufzufüllen {$DEFINE DichteNichtnegativieren:=perEinzelklau} // negative Dichten aus zufälligen Phasenraumpunkten auffüllen { $DEFINE negativeDichteueberwachung} // prüfen, ob die Dichten negativ werden { $DEFINE exzessiveArrayBereichsTests} // zeitaufwändige(!) Array-Bereichstest für fftw interface uses Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix, fftw; type tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier,zfRungeKuttaZehn,zfRungeKuttaZwoelf,zfRungeKuttaVierzehn); tVerteilungsfunktion = function(x: double): double; tEMFeldGroesze = ( efAX,efAY, efDAXDT,efDAYDT, efEX,efBZ ); tEMQuellGroesze = ( eqRho,eqJX,eqJY ); tMaterieFeldGroesze = (mfPY); tMaterieSpeicherGroesze = (msPX,msPXSqr,msPY,msVX,msVY,msN,msPXRipple); tTeilchenSpeziesGroeszen = (tsgMasse,tsgIQdrMasse,tsgLadung); tEMQuellen = array[tEMQuellGroesze] of double; const erstesEMFMitAbleitung = efAX; letztesEMFMitAbleitung = efDAYDT; type tFelder = class; tGitter = class; tSimulation = class; { tFeldAusgabeDatei } tFeldAusgabeDatei = class private ableitung: boolean; emFeld: tEMFeldGroesze; emQuelle: tEMQuellGroesze; matFeld: tMaterieSpeicherGroesze; wasSpeichern, // 0: EM-Feld; 1: EM-Quelle; 2: Mat-Feld teilchen, // -1: Gesamt-Feld; sonst: Feld des entsprechenden Teilchen nNum,tCnt: longint; pre,suf: string; datei: file; pro: tProtokollant; buf: pointer; bufLen,bufPos: longint; procedure schreibeBuffer(minBufLen: longint); public zeitAnz: longint; sDT: double; constructor create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); destructor destroy; override; function dump: string; procedure schreibeKopf; procedure speichereWerte(gitter: tGitter; sT, sDX: double); procedure verzeichnisAufraeumen; end; { tImpulsRaumAusgabeDatei } tImpulsRaumAusgabeDatei = class private teilchen,nNum,tCnt: longint; pre,suf: string; pro: tProtokollant; public zeitAnz: longint; sDT: double; constructor create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); destructor destroy; override; function dump: string; procedure naechsterZykel; procedure speichereWerte(gitter: tGitter; sT: double); procedure verzeichnisAufraeumen; end; { tTeilchenSpezies } tTeilchenSpezies = class private ortsGrenzen: array of double; dichteStuecke: array of string; public nMin,nMax,breite: double; // minimale und maximale Massendichte in nc; Breite der Impulsverteilung eigenschaften: array[tTeilchenSpeziesGroeszen] of double; // q in e; m in me diffusion: array[0..1] of double; // Diffusionsterm in x- und p-Richtung (= maximales d(ln n)/d[xp] * iD[XP]) constructor create; overload; constructor create(original: tTeilchenSpezies); overload; destructor destroy; override; procedure dump(prot: tProtokollant; prefix: string); function gibDichte(x: double; kvs: tKnownValues): double; // Massendichte procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant); end; { TFelder } (* procedure berechneAbleitungen(pX: double); inline; function nichtnegativieren: double; inline; *) tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen private // 2d-Arrays (x & px) sind wie folgt indiziert: [iP + iX*aP] teilchen: array of tTeilchenSpezies; lichters: array[boolean] of tMyStringlist; massen: array of double; // Gesamtmassen je Teilchenspezies zur Überwachung der Erhaltungsgrößen gesamtDefizit: double; impulsRaumGradient: array of array[0..1] of pDouble; // Ableitung des Impulsraumes nach x und px iMGamma: array of pDouble; // 1/m/gamma (= v/p) _aX,_aP: longint; _dX,_dP,iDX: double; ffts: array of array[0..1,boolean] of fftw_plan_double; // fft-Algorithmen fftTmp,lFftTmp: Pcomplex_double; // Zwischenspeicher für ffts transps: array of fftw_plan_double; // Algorithmen zum Transponieren des Phasenraumes transpTmp: pDouble; // Zwischenspeicher dafür procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline; {$IFDEF exzessiveArrayBereichsTests} procedure pruefeArrayEnden(fehler: string); {$ENDIF} {$IFDEF dichteNichtnegativieren} procedure nichtnegativieren; inline; {$ENDIF} procedure berechnePhasenraumAbleitungen; inline; public emFelder: array[tEMFeldGroesze,boolean] of pDouble; // EM-Felder und deren Zeitableitungen // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier! emQuellen: array[tEMQuellGroesze] of pDouble; // Quellen für die EM-Felder matFelder: array of array[tMaterieFeldGroesze] of pDouble; // px-unabhängige Felder, nach Teilchenspezies sortiert impulsraum: array of array[boolean] of pDouble; // px-abhängige Felder, nach Teilchenspezies sortiert, sowie deren Ableitung gitter: tGitter; constructor create(anzX,anzP: longint; deltaX,deltaP: double; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); destructor destroy; override; procedure berechneAbleitungen; procedure setzeNull; inline; procedure dumpErhaltungsgroeszen; property aX: longint read _aX; property aP: longint read _aP; property dX: double read _dX; property dP: double read _dP; function impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): double; overload; inline; function impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): double; overload; inline; {$DEFINE LiKoInterface} {$INCLUDE linearkombinationen.inc} {$UNDEF LiKoInterface} procedure schreibePhasenraum(tlc: longint; datNamPre, datNamSuf: string); function transponiereFuer(tlc: longint): pDouble; end; { tGitter } tGitter = class(tObject) private simulation: tSimulation; prot: tProtokollant; zeitverfahren: tZeitverfahren; abbruch: boolean; dX,iDX,xl,t: double; kvs: tKnownValues; {$IF dichteNichtnegativieren = perEinzelklau} impulsRaumPermutationen: array of tLongintArray; {$ENDIF} procedure abbrechen; public aktuelleFelder: longint; felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern constructor create(besitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; name: string); destructor destroy; override; procedure iteriereSchritt(dT: double); procedure dumpErhaltungsgroeszen; {$IF dichteNichtnegativieren = perEinzelklau} function gibImpulsRaumPermutation: pTLongintArray; {$ENDIF} end; { tSimulation } tSimulation = class(tObject) private prot: tProtokollant; gitter: tGitter; dT,tEnde,sT,sDT,sDX: double; fortschrittsAnzeige: boolean; feldAusgabeDateien: array of tFeldAusgabeDatei; impulsRaumAusgabeDateien: array of tImpulsRaumAusgabeDatei; transpTmp: array of double; public gotSigusr1,gotSigterm, gotSigint: boolean; constructor create(inName: string; protokollant: tProtokollant; name: string); destructor destroy; override; function iteriereSchritt(start: double; var zeitPhysik,zeitDatei: double): boolean; // noch nicht zu Ende? end; procedure SignalCapture(signal : longint); cdecl; function fourierWindow(x: extended): extended; inline; implementation const emFeldNamen: array[tEMFeldGroesze] of string = ( 'AX','AY','DAXDT','DAYDT','EX','BZ' ); emQuellNamen: array[tEMQuellGroesze] of string = ( 'RHO','JX','JY' ); matSpeicherNamen: array[tMaterieSpeicherGroesze] of string = ( 'PX','PXSQR','PY','VX','VY','N','PXRIPPLE' ); minAusgabeBuffer = 1024*1024; var simulationen: array of tSimulation; // tFeldAusgabeDatei *********************************************************** constructor tFeldAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); var emF: tEMFeldGroesze; emQ: tEMQuellGroesze; maF: tMaterieSpeicherGroesze; 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:=-1; emFeld:=efAX; emQuelle:=eqRho; matFeld:=msN; wasSpeichern:=0; 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(tEMFeldGroesze) to high(tEMFeldGroesze) do if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin emFeld:=emF; wasSpeichern:=0; ableitung:=abl; teilchen:=-2; gef:=true; break; end; end; if not gef then for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do if emQuellNamen[emQ] = feldName then begin ableitung:=false; emQuelle:=emQ; wasSpeichern:=1; gef:=true; break; end; if not gef then for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do if matSpeicherNamen[maF] = feldName then begin ableitung:=false; matFeld:=maF; wasSpeichern:=2; gef:=true; break; end; if not gef then begin pro.schreibe('tFeldAusgabeDatei.create kennt Größe '''+feldName+''' nicht!',true); nam:=''; for abl:=false to true do for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do nam:=nam+', '+copy('D',1,byte(abl))+emFeldNamen[emF]; for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do nam:=nam+', '+emQuellNamen[emQ]; for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do nam:=nam+', '+matSpeicherNamen[maF]; delete(nam,1,2); pro.schreibe('Ich kenne nur: '+nam,true); raise exception.create('Fehler in tFeldAusgabeDatei.create!'); end; case wasSpeichern of 0: pre:=emFeldNamen[emFeld]; 1: pre:=emQuellNamen[emQuelle]; 2: pre:=matSpeicherNamen[matFeld]; end; if teilchen>=0 then pre:=pre+inttostr(teilchen+1); if ableitung then pre:='d'+pre; bufLen:=0; buf:=nil; bufPos:=0; pre:=prefix+pre; suf:=suffix; end; destructor tFeldAusgabeDatei.destroy; begin if bufPos<>0 then schreibeBuffer(0); freeMem(buf,bufLen); 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); raise exception.create('Fehler in tFeldAusgabeDatei.destroy!'); end; inherited destroy; end; procedure tFeldAusgabeDatei.schreibeBuffer(minBufLen: longint); begin if buf=nil then begin bufLen:=max(minAusgabeBuffer,minBufLen); getmem(buf,bufLen); end else begin reset(datei,1); seek(datei,fileSize(datei)); blockWrite(datei,buf^,bufPos); end; bufPos:=0; end; function tFeldAusgabeDatei.dump: string; begin case wasSpeichern of 0: result:=emFeldNamen[emFeld]; 1: result:=emQuellNamen[emQuelle]; 2: result:=matSpeicherNamen[matFeld]; end; if teilchen>=0 then result:=result+'['+inttostr(teilchen+1)+']'; if ableitung then result:=result+''''; result:=result+' -> '+pre+'-$i'+suf; end; procedure tFeldAusgabeDatei.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); raise exception.create('Fehler in tFeldAusgabeDatei.schreibeKopf!'); end; if bufPos<>0 then schreibeBuffer(0); tCnt:=0; assignfile(datei,pre+'-'+inttostr(nNum)+suf); inc(nNum); rewrite(datei,1); blockwrite(datei,nNum,sizeof(longint)); blockwrite(datei,zeitAnz,sizeof(longint)); closefile(datei); end; procedure tFeldAusgabeDatei.speichereWerte(gitter: tGitter; sT, sDX: double); var i,cnt: longint; sX,cX,val: double; begin if (teilchen>=length(gitter.felders[gitter.aktuelleFelder].teilchen)) then begin pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true); raise exception.create('Fehler in tFeldAusgabeDatei.speichereWerte!'); end; if sT>=nNum then schreibeKopf; cnt:=floor(((gitter.felders[gitter.aktuelleFelder].aX-1)*gitter.dX+Min(gitter.dX,sDX)/2)/sDX+1); cX:=gitter.xl; sX:=cX+(cnt-1)*sDX; inc(tCnt); if bufLen-bufPos < sizeof(longint) + (2+cnt)*sizeof(double) then schreibeBuffer(sizeof(longint) + (2+cnt)*sizeof(double)); if bufLen-bufPos < sizeof(longint) + (2+cnt)*sizeof(double) then begin pro.schreibe('Schreibbuffer ist zu klein ('+inttostr(bufLen)+' statt mindestens '+inttostr(sizeof(longint) + (2+cnt)*sizeof(double))+' Bytes)!'); raise exception.create('Fehler in tFeldAusgabeDatei.speichereWerte!'); end; move(cX,(buf+bufPos)^,sizeof(double)); inc(bufPos,sizeof(double)); move(sX,(buf+bufPos)^,sizeof(double)); inc(bufPos,sizeof(double)); move(cnt,(buf+bufPos)^,sizeof(longint)); inc(bufPos,sizeof(longint)); sX:=cX-Min(gitter.dX,sDX)/2; case wasSpeichern of 0: // em-Feld speichern for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin while cX>=sX do begin dec(cnt); move((gitter.felders[gitter.aktuelleFelder].emFelder[emFeld,ableitung]+i)^,(buf+bufPos)^,sizeof(double)); inc(bufPos,sizeof(double)); sX:=sX+sDX; end; cX:=cX+gitter.dX; end; 1: // em-Quelle for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin while cX>=sX do begin dec(cnt); val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,emQuelle); move(val,(buf+bufPos)^,sizeof(double)); inc(bufPos,sizeof(double)); sX:=sX+sDX; end; cX:=cX+gitter.dX; end; 2: // Materiefeld for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin while cX>=sX do begin dec(cnt); val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,matFeld); move(val,(buf+bufPos)^,sizeof(double)); inc(bufPos,sizeof(double)); sX:=sX+sDX; end; cX:=cX+gitter.dX; end; end{of case}; if cnt<>0 then begin pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true); raise exception.create('Fehler in tFeldAusgabeDatei.speichereWerte!'); end; end; procedure tFeldAusgabeDatei.verzeichnisAufraeumen; var sr: tSearchRec; err: longint; dats: tStringList; pfad: string; i: longint; begin dats:=tStringList.create; pfad:=extractfilepath(pre+'00'+suf); err:=findFirst(pfad+'*',$00,sr); while err=0 do begin if (sr.name<>'.') and (sr.name<>'..') then dats.add(pfad+sr.name); err:=findNext(sr); end; findClose(sr); for i:=0 to dats.count-1 do if (rightStr(dats[i],4)<>'.dat') and (rightStr(dats[i],4)<>'.png') and (rightStr(dats[i],4)<>'.bmp') then begin pro.schreibe('Zu löschende Datei ist keine .dat, .bmp oder .png!',true); raise exception.create('Zu löschende Datei ist keine .dat, .bmp oder .png!'); end; for i:=0 to dats.count-1 do deleteFile(dats[i]); dats.free; end; // tImpulsRaumAusgabeDatei ***************************************************** constructor tImpulsRaumAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string); var num: longint; 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 begin pro.schreibe('tImpulsRaumAusgabeDatei.create kann keine über Teilchen integrierten Größen abspeichern!',true); raise exception.create('Fehler in tImpulsRaumAusgabeDatei.create!'); end; if feldName<>'N' then begin pro.schreibe('tImpulsRaumAusgabeDatei.create kennt Größe '''+feldName+''' nicht!',true); pro.schreibe('Ich kenne nur N!',true); raise exception.create('Fehler in tImpulsRaumAusgabeDatei.create!'); end; nNum:=0; // Header 0 wurde also noch nicht geschrieben zeitAnz:=-1; tCnt:=-1; pre:='N'; if teilchen>=0 then pre:=pre+inttostr(teilchen+1); pre:=prefix+pre; suf:=suffix; end; destructor tImpulsRaumAusgabeDatei.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); raise exception.create('Fehler in tImpulsRaumAusgabeDatei.destroy!'); end; inherited destroy; end; function tImpulsRaumAusgabeDatei.dump: string; begin result:='N'; if teilchen>=0 then result:=result+'['+inttostr(teilchen+1)+']'; result:=result+' -> '+pre+'-$i_$j'+suf; end; procedure tImpulsRaumAusgabeDatei.naechsterZykel; 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); raise exception.create('Fehler in tImpulsRaumAusgabeDatei.schreibeKopf!'); end; tCnt:=0; inc(nNum); end; procedure tImpulsRaumAusgabeDatei.speichereWerte(gitter: tGitter; sT: double); var tmpI: longint; tmpD: double; sF: tFelder; datei: file; begin if (teilchen>=length(gitter.felders[gitter.aktuelleFelder].teilchen)) then begin pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true); raise exception.create('Fehler in tImpulsRaumAusgabeDatei.speichereWerte!'); end; if sT>=nNum then naechsterZykel; assignFile(datei,pre+'-'+inttostr(nNum-1)+'-'+myinttostr(tCnt,ceil(ln(zeitAnz)/ln(10)),'0')+suf); rewrite(datei,1); inc(tCnt); sF:=gitter.felders[gitter.aktuelleFelder]; tmpD:=-(sF.aP-1)*sF.dP/2; blockWrite(datei,tmpD,sizeof(double)); tmpD:=-tmpD; blockWrite(datei,tmpD,sizeof(double)); tmpI:=sF.aP; blockWrite(datei,tmpI,sizeof(longint)); blockWrite(datei,gitter.xl,sizeof(double)); tmpD:=gitter.xl + (sF.aX-1)*sF.dX; blockWrite(datei,tmpD,sizeof(double)); tmpI:=sF.aX; blockWrite(datei,tmpI,sizeof(longint)); blockWrite(datei,sF.transponiereFuer(teilchen)^,sizeof(double)*sF.aX*sF.aP); closeFile(datei); end; procedure tImpulsRaumAusgabeDatei.verzeichnisAufraeumen; var sr: tSearchRec; err: longint; dats: tStringList; pfad: string; i: longint; begin dats:=tStringList.create; pfad:=extractfilepath(pre+'00'+suf); err:=findFirst(pfad+'*',$00,sr); while err=0 do begin if (sr.name<>'.') and (sr.name<>'..') then dats.add(pfad+sr.name); err:=findNext(sr); end; findClose(sr); for i:=0 to dats.count-1 do if (rightStr(dats[i],4)<>'.dat') and (rightStr(dats[i],4)<>'.png') and (rightStr(dats[i],4)<>'.bmp') then begin pro.schreibe('Zu löschende Datei ist keine .dat, .bmp oder .png!',true); raise exception.create('Zu löschende Datei ist keine .dat, .bmp oder .png!'); end; for i:=0 to dats.count-1 do deleteFile(dats[i]); dats.free; end; // tTeilchenSpezies ************************************************************ constructor tTeilchenSpezies.create; var egr: tTeilchenSpeziesGroeszen; begin inherited create; fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0); setlength(ortsGrenzen,0); fillchar(dichteStuecke,sizeof(dichteStuecke),#0); setlength(dichteStuecke,0); nMin:=1E-3; nMax:=nMin; for egr:=low(tTeilchenSpeziesGroeszen) to high(tTeilchenSpeziesGroeszen) do eigenschaften[egr]:=byte(egr in [tsgMasse,tsgIQdrMasse]); diffusion[0]:=0; diffusion[1]:=0; breite:=1; end; constructor tTeilchenSpezies.create(original: tTeilchenSpezies); var egr: tTeilchenSpeziesGroeszen; i: longint; begin create; nMin:=original.nMin; nMax:=original.nMax; for egr:=low(tTeilchenSpeziesGroeszen) to high(tTeilchenSpeziesGroeszen) do eigenschaften[egr]:=original.eigenschaften[egr]; for i:=0 to 1 do diffusion[i]:=original.diffusion[i]; breite:=original.breite; fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0); setlength(ortsGrenzen,length(original.ortsGrenzen)); for i:=0 to length(ortsGrenzen)-1 do ortsGrenzen[i]:=original.ortsGrenzen[i]; fillchar(dichteStuecke,sizeof(dichteStuecke),#0); setlength(dichteStuecke,length(original.dichteStuecke)); for i:=0 to length(dichteStuecke)-1 do dichteStuecke[i]:=original.dichteStuecke[i]; 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+'m = '+floattostr(eigenschaften[tsgMasse])+' me'); prot.schreibe(prefix+'q = '+floattostr(eigenschaften[tsgLadung])+' e'); prot.schreibe(prefix+'dLnN/dX < '+floattostr(diffusion[0])+' / dX'); prot.schreibe(prefix+'dLnN/dP < '+floattostr(diffusion[1])+' / dP'); 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); raise exception.create('Fehler in tTeilchenSpezies.liesDichteFunktion!'); 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); raise exception.create('Fehler in tTeilchenSpezies.liesDichteFunktion!'); 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; // tFelder ********************************************************************* constructor tFelder.create(anzX,anzP: longint; deltaX,deltaP: double; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); var i,j: longint; rechts,abl: boolean; dens: double; emF: tEMFeldGroesze; emQ: tEMQuellGroesze; maF: tMaterieFeldGroesze; dims: Pfftw_iodim; begin inherited create; gitter:=parent; gesamtDefizit:=0; _aX:=anzX; _aP:=anzP; _dX:=deltaX; _dP:=deltaP; iDX:=1/dX; fillchar(teilchen,sizeof(teilchen),#0); setlength(teilchen,length(_teilchen)); for i:=0 to length(teilchen)-1 do teilchen[i]:=tTeilchenSpezies.create(_teilchen[i]); fillchar(impulsRaumGradient,sizeof(impulsRaumGradient),#0); setlength(impulsRaumGradient,length(teilchen)); for i:=0 to length(impulsRaumGradient)-1 do for j:=0 to length(impulsRaumGradient[i])-1 do begin fftw_getmem(impulsRaumGradient[i,j],(aP*aX+1)*sizeof(double)); fillchar(impulsRaumGradient[i,j]^,aP*aX*sizeof(double),#0); (impulsRaumGradient[i,j]+aP*aX)^:=pi; end; fillchar(iMGamma,sizeof(iMGamma),#0); setlength(iMGamma,length(teilchen)); for i:=0 to length(iMGamma)-1 do begin fftw_getmem(iMGamma[i],aP*aX*sizeof(double)); fillchar(iMGamma[i]^,aP*aX*sizeof(double),#0); end; for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do begin emFelder[emF,true]:=nil; for abl:=false to not(emF in [efEX,efBZ]) do begin fftw_getmem(emFelder[emF,abl],aX*sizeof(double)); fillchar(emFelder[emF,abl]^,aX*sizeof(double),#0); end; end; for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do begin fftw_getmem(emQuellen[emQ],aX*sizeof(double)); fillchar(emQuellen[emQ]^,aX*sizeof(double),#0); end; fillchar(matFelder,sizeof(matFelder),#0); setlength(matFelder,length(teilchen)); for i:=0 to length(matFelder)-1 do for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do begin fftw_getmem(matFelder[i,maF],aX*sizeof(double)); fillchar(matFelder[i,maF]^,aX*sizeof(double),#0); end; fillchar(impulsraum,sizeof(impulsraum),#0); setlength(impulsraum,length(teilchen)); for i:=0 to length(impulsraum)-1 do for abl:=false to true do begin fftw_getmem(impulsraum[i,abl],(aX*aP+1)*sizeof(double)); fillchar(impulsraum[i,abl]^,aX*aP*sizeof(double),#0); (impulsraum[i,abl]+aX*aP)^:=pi; end; i:=max((aP div 2 + 1)*aX,(aX div 2 + 1)*aP); fftw_getmem(fftTmp,(i+1)*sizeof(complex_double)); fillchar(fftTmp^,(i+1)*sizeof(complex_double),#0); lFftTmp:=fftTmp+i; lFftTmp^.re:=pi; lFftTmp^.im:=sqr(pi); fillchar(ffts,sizeof(ffts),#0); // Vorsicht, die Planung der FFTs überschreibt den Input, auch wenn die eigentlichen FFTs das nicht tun! setlength(ffts,length(teilchen)); for i:=0 to length(ffts)-1 do begin ffts[i,0,false]:= // Planung der Hintransformationen über x: fftw_plan_many_dft_r2c( 1, // jeweils 1-d @_aX, // jeweils _aX Werte aP, // aP Mal impulsraum[i,false], // von dort nil, // kein Padding aP, // Abstand benachbarter Werte 1, // Verschiebung der Werte zwischen den Transformationen fftTmp, // nach dort nil, // kein Padding aP, // Abstand benachbarter Werte 1, // Verschiebung der Werte zwischen den Transformationen [fftw_preserve_input,fftw_exhaustive] // input nicht überschreiben, erschöpfende Suche nach FFT ); ffts[i,0,true]:= // Planung der Rücktransformationen über x: fftw_plan_many_dft_c2r( 1, // jeweils 1-d @_aX, // jeweils _aX Werte aP, // aP Mal fftTmp, // von hier nil, // kein Padding aP, // Abstand benachbarter Werte 1, // Verschiebung der Werte zwischen den Transformationen impulsRaumGradient[i,0], // nach dort nil, // kein Padding aP, // Abstand benachbarter Werte 1, // Verschiebung der Werte zwischen den Transformationen [fftw_destroy_input,fftw_exhaustive] // input überschreiben, erschöpfende Suche nach FFT ); ffts[i,1,false]:= // Planung der Hintransformationen über p: fftw_plan_many_dft_r2c( 1, // jeweils 1-d @_aP, // jeweils _aP Werte aX, // aX Mal impulsraum[i,false], // von dort nil, // kein Padding 1, // Abstand benachbarter Werte aP, // Verschiebung der Werte zwischen den Transformationen fftTmp, // nach dort nil, // kein Padding 1, // Abstand benachbarter Werte aP div 2 + 1, // Verschiebung der Werte zwischen den Transformationen [fftw_preserve_input,fftw_exhaustive] // input nicht überschreiben, erschöpfende Suche nach FFT ); ffts[i,1,true]:= // Planung der Rücktransformationen über p: fftw_plan_many_dft_c2r( 1, // jeweils 1-d @_aP, // jeweils _aP Werte aX, // aX Mal fftTmp, // von dort nil, // kein Padding 1, // Abstand benachbarter Werte aP div 2 + 1, // Verschiebung der Werte zwischen den Transformationen impulsRaumGradient[i,1], // nach dort nil, // kein Padding 1, // Abstand benachbarter Werte aP, // Verschiebung der Werte zwischen den Transformationen [fftw_destroy_input,fftw_exhaustive] // input überschreiben, erschöpfende Suche nach FFT ); end; (* fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, double *in, const int *inembed, int istride, int idist, fftw_complex *out, const int *onembed, int ostride, int odist, unsigned flags); fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany, fftw_complex *in, const int *inembed, int istride, int idist, double *out, const int *onembed, int ostride, int odist, unsigned flags); *) fillchar(transps,sizeof(transps),#0); setlength(transps,length(teilchen)); fftw_getmem(transpTmp,(aX*aP+1)*sizeof(double)); (transpTmp+aX*aP)^:=pi; getmem(dims,2*sizeof(fftw_iodim)); dims^.n:=aX; // aX Spalten dims^.ins:=aP; // Werte im Input um jeweils aP versetzt dims^.ous:=1; // und im Output um jeweils 1 (dims+1)^.n:=aP; // aP Zeilen (dims+1)^.ins:=1; // Werte im Input um jeweils 1 versetzt (dims+1)^.ous:=aX; // und im Output um jeweils aX for i:=0 to length(transps)-1 do begin transps[i]:= fftw_plan_guru_r2r( 0, // null Dimensionen für die FFT, also nur kopieren! nil, // Bockwurst 2, // 2 Dimensionen dims, // @(aX;aP) impulsraum[i,false], // von hier transpTmp, // nach dort nil, // Bockwurst [fftw_estimate] // einfach irgendwas nehmen ); end; freemem(dims); (* fftw_plan fftw_plan_guru_r2r(int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags); *) {$IFDEF exzessiveArrayBereichsTests} pruefeArrayEnden('FFT/Transponieren schreibt beim Planen über das Ende hinaus!'); {$ENDIF} fillchar(massen,sizeof(massen),#0); setlength(massen,length(teilchen)); for i:=0 to length(massen)-1 do massen[i]:=0; for i:=0 to aX-1 do for j:=0 to length(teilchen)-1 do begin dens:=teilchen[j].gibDichte(gitter.xl+(i-2)*gitter.dX,parent.kvs); initialisiereDichte(i,j,teilchen[j].breite,dens); massen[j]:=massen[j]+dens*gitter.dX*teilchen[j].eigenschaften[tsgMasse]; 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,j: longint; abl: boolean; emF: tEMFeldGroesze; emQ: tEMQuellGroesze; maF: tMaterieFeldGroesze; begin setlength(massen,0); for i:=0 to length(teilchen)-1 do teilchen[i].free; setlength(teilchen,0); for i:=0 to length(impulsRaumGradient)-1 do for j:=0 to length(impulsRaumGradient[i])-1 do fftw_freemem(impulsRaumGradient[i,j]); setlength(iMGamma,0); for i:=0 to length(iMGamma)-1 do fftw_freemem(iMGamma[i]); setlength(iMGamma,0); for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do for abl:=false to not(emF in [efEX,efBZ]) do fftw_freemem(emFelder[emF,abl]); for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do fftw_freemem(emQuellen[emQ]); for i:=0 to length(matFelder)-1 do for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do fftw_freemem(matFelder[i,maF]); setlength(matFelder,0); for i:=0 to length(impulsraum)-1 do for abl:=false to true do fftw_freemem(impulsraum[i,abl]); setlength(impulsraum,0); lichters[false].free; lichters[true].free; for i:=0 to length(ffts)-1 do for j:=0 to 1 do for abl:=false to true do fftw_destroy_plan(ffts[i,j,abl]); setlength(ffts,0); fftw_freemem(fftTmp); for i:=0 to length(transps)-1 do fftw_destroy_plan(transps[i]); setlength(transps,0); fftw_freemem(transpTmp); inherited destroy; end; procedure tFelder.initialisiereDichte(ort,tlc: longint; breite,n: double); var i: longint; ges: double; begin ges:=0; for i:=0 to aP-1 do begin (impulsraum[tlc,false]+i+ort*aP)^:=power(1/2,sqr((i-aP/2)*dP*2/breite)); ges:=ges + (impulsraum[tlc,false]+i+ort*aP)^; end; ges:=n/ges/dP; for i:=0 to aP-1 do (impulsraum[tlc,false]+i+ort*aP)^:=(impulsraum[tlc,false]+i+ort*aP)^ * ges; end; {$IFDEF exzessiveArrayBereichsTests} procedure tFelder.pruefeArrayEnden(fehler: string); var i,j: longint; abl: boolean; begin if abs((transpTmp+aP*aX)^-pi)>1e-10 then begin gitter.prot.schreibe(fehler); gitter.prot.schreibe(floattostr((transpTmp+aP*aX)^)+' vs. '+floattostr(pi)); raise exception.create('Fehler in tFelder.pruefeArrayEnden!'); end; if (abs(lFftTmp^.re-pi)>1e-10) or (abs(lFftTmp^.im-sqr(pi))>1e-10) then begin gitter.prot.schreibe(fehler); gitter.prot.schreibe(floattostr(lFftTmp^.re)+' vs. '+floattostr(pi)+' und '+floattostr(lFftTmp^.im)+' vs. '+floattostr(sqr(pi))); raise exception.create('Fehler in tFelder.pruefeArrayEnden!'); end; for i:=0 to length(impulsRaumGradient)-1 do for j:=0 to 1 do if abs((impulsRaumGradient[i,j]+aP*aX)^-pi)>1e-10 then begin gitter.prot.schreibe(fehler); gitter.prot.schreibe(floattostr((impulsRaumGradient[i,j]+aP*aX)^)+' vs. '+floattostr(pi)); raise exception.create('Fehler in tFelder.pruefeArrayEnden!'); end; for i:=0 to length(impulsRaum)-1 do for abl:=false to true do if abs((impulsRaum[i,abl]+aP*aX)^-pi)>1e-10 then begin gitter.prot.schreibe(fehler); gitter.prot.schreibe(floattostr((impulsRaum[i,abl]+aP*aX)^)+' vs. '+floattostr(pi)); raise exception.create('Fehler in tFelder.pruefeArrayEnden!'); end; end; {$ENDIF} procedure tFelder.berechneAbleitungen; var i,j,k: longint; emQ: tEMQuellGroesze; emF: tEMFeldGroesze; rechts: boolean; tmp: double; begin // PX ggf. korrigieren for i:=0 to length(teilchen)-1 do // linker Rand for j:=0 to (aP div 2)-1 do // negativer Impuls if (impulsraum[i,false]+j)^>0 then begin (impulsraum[i,false]+aP-j-byte(j=0))^:=(impulsraum[i,false]+aP-j-byte(j=0))^ + (impulsraum[i,false]+j)^; (impulsraum[i,false]+j)^:=0; end; for i:=0 to length(teilchen)-1 do // rechter Rand for j:=(aP div 2)+1 to aP-1 do // positiver Impuls if (impulsraum[i,false]+j+(aX-1)*aP)^>0 then begin (impulsraum[i,false]+aP-j+(aX-1)*aP)^:=(impulsraum[i,false]+aP-j+(aX-1)*aP)^ + (impulsraum[i,false]+j+(aX-1)*aP)^; (impulsraum[i,false]+j+(aX-1)*aP)^:=0; end; for i:=0 to length(teilchen)-1 do for j:=0 to aX-1 do begin if (impulsraum[i,false]+j*aP)^>0 then begin // unterer Rand (impulsraum[i,false]+1+j*aP)^:=(impulsraum[i,false]+1+j*aP)^+(impulsraum[i,false]+j*aP)^; (impulsraum[i,false]+j*aP)^:=0; end; if (impulsraum[i,false]+aP-1+j*aP)^>0 then begin // oberer Rand (impulsraum[i,false]+aP-2+j*aP)^:=(impulsraum[i,false]+aP-2+j*aP)^+(impulsraum[i,false]+aP-1+j*aP)^; (impulsraum[i,false]+aP-1+j*aP)^:=0; end; end; // PY und iMGamma berechnen for i:=0 to length(matFelder)-1 do for j:=0 to aX-1 do begin (matFelder[i,mfPY]+j)^:=-teilchen[i].eigenschaften[tsgLadung] * (emFelder[efAY,false]+j)^; for k:=0 to aP-1 do (iMGamma[i]+k+j*aP)^:= 1/(sqrt(1 + (sqr((k-(aP div 2))*dP) + sqr((matFelder[i,mfPY]+j)^))*teilchen[i].eigenschaften[tsgIQdrMasse])*teilchen[i].eigenschaften[tsgMasse]); // 1/(m Gamma) = 1/sqrt(1 + (pX/m)^2 + (pY/m)^2)*m end; // rho und j berechnen for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do for i:=0 to aX-1 do (emQuellen[emQ]+i)^:=0; for i:=0 to length(matFelder)-1 do for j:=0 to aX-1 do for k:=0 to aP-1 do begin tmp:=(impulsraum[i,false]+k+j*aP)^ * teilchen[i].eigenschaften[tsgLadung] * dP; (emQuellen[eqRho]+j)^:=(emQuellen[eqRho]+j)^ + tmp; // * iGamma[i] ??? TODO ... falls, dann aber ohne M! (emQuellen[eqJX]+j)^:= (emQuellen[eqJX]+j)^ + tmp * (k-(aP div 2))*dP * (iMGamma[i]+k+j*aP)^; (emQuellen[eqJY]+j)^:= (emQuellen[eqJY]+j)^ + tmp * (matFelder[i,mfPY]+j)^ * (iMGamma[i]+k+j*aP)^; end; // E und B aus Rho und A berechnen // dEX/dx = rho ... hier muss integriert werden: // EX = EX(x- deltax) + * deltax emFelder[efEX,false]^:=emQuellen[eqRho]^ * dX; // linker Rand for i:=1 to aX-1 do // Rest (emFelder[efEX,false]+i)^:= (emFelder[efEX,false]+i-1)^ + ((emQuellen[eqRho]+i-1)^ + (emQuellen[eqRho]+i)^) * dX/2; // B = rot A emFelder[efBZ,false]^:=0; // linker Rand for i:=1 to aX-2 do // Mitte (emFelder[efBZ,false]+i)^:= ((emFelder[efAY,false]+i+1)^ - (emFelder[efAY,false]+i-1)^)*iDX/2; (emFelder[efBZ,false]+aX-1)^:=0; // rechter Rand // erste Ableitung des A-Feldes berechnen // dA/dt = dA/dt - E for i:=1 to aX-2 do (emFelder[efAX,true]+i)^:=(emFelder[efDAXDT,false]+i)^ - (emFelder[efEX,false]+i)^; move(emFelder[efDAYDT,false]^,emFelder[efAY,true]^,aX*sizeof(double)); // Randbedingungen für die Ableitungen des A-Felds setzen for emF:=efAX to efAY do begin // Vakuumrandbedingungen für das A-Feld emFelder[emF,true]^:= ((emFelder[emF,false]+1)^ - emFelder[emF,false]^)*iDX; (emFelder[emF,true]+aX-1)^:= ((emFelder[emF,false]+aX-2)^ - (emFelder[emF,false]+aX-1)^)*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 (emFelder[efAY,true]+(aX-1)*byte(rechts))^:= (emFelder[efAY,true]+(aX-1)*byte(rechts))^ + exprToFloat(false,lichters[rechts][i],gitter.kvs,nil); // zweite Ableitung des A-Feldes berechnen // d2A/dt2 = Laplace(A) - j ( - dE/dt wird auf dA/dt direkt aufgeschlagen !!! ) for i:=1 to aX-2 do (emFelder[efDAXDT,true]+i)^:= ( (emFelder[efAX,false]+i+1)^ - 2*(emFelder[efAX,false]+i)^ + (emFelder[efAX,false]+i-1)^ )*sqr(iDX) - (emQuellen[eqJX]+i)^; for i:=1 to aX-2 do (emFelder[efDAYDT,true]+i)^:= ( (emFelder[efAY,false]+i+1)^ - 2*(emFelder[efAY,false]+i)^ + (emFelder[efAY,false]+i-1)^ )*sqr(iDX) - (emQuellen[eqJY]+i)^; // Gradienten der Phasenraumbesetzungsdichte berechnen berechnePhasenraumAbleitungen; // Zeitableitung der Phasenraumbesetzungsdichte berechnen for i:=0 to length(teilchen)-1 do for j:=1 to aX-2 do for k:=1 to aP-2 do (impulsraum[i,true]+k+j*aP)^:= // df/dt = - (k-(aP div 2))*dP * (iMGamma[i]+k+j*aP)^ // - v * (impulsRaumGradient[i,0]+k+j*aP)^ // * Nabla f - teilchen[i].eigenschaften[tsgLadung] * // - q( ( (emFelder[efEX,false]+j)^ // E + (matFelder[i,mfPY]+j)^ * (iMGamma[i]+k+j*aP)^ * (emFelder[efBZ,false]+j)^) // + v/c x B) * (impulsRaumGradient[i,1]+k+j*aP)^; // * Nabla_p f end; procedure tFelder.setzeNull; var i,j,k: longint; abl: boolean; emF: tEMFeldGroesze; emQ: tEMQuellGroesze; maF: tMaterieFeldGroesze; begin for i:=0 to length(impulsRaumGradient)-1 do for j:=0 to length(impulsRaumGradient[i])-1 do for k:=0 to aX*aP-1 do (impulsRaumGradient[i,j]+k)^:=0; for i:=0 to length(iMGamma)-1 do for j:=0 to aX*aP-1 do (iMGamma[i]+j)^:=1; for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do for abl:=false to not(emF in [efEX,efBZ]) do for i:=0 to aX-1 do (emFelder[emF,abl]+i)^:=0; for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do for i:=0 to aX-1 do (emQuellen[emQ]+i)^:=0; for i:=0 to length(matFelder)-1 do for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do for j:=0 to aX-1 do (matFelder[i,maF]+j)^:=0; for i:=0 to length(impulsraum)-1 do for abl:=false to true do for j:=0 to aX*aP-1 do (impulsraum[i,abl]+j)^:=0; end; procedure tFelder.dumpErhaltungsgroeszen; var i,j: longint; dens: double; begin if length(massen)>0 then begin dens:=0; for i:=0 to length(massen)-1 do dens:=dens+massen[i]/teilchen[i].eigenschaften[tsgMasse]; gitter.prot.schreibe('Gesamtdefizit: '+floattostr(gesamtDefizit)+' (Anteil '+floattostr(gesamtDefizit/dens*dX*dP)+')',true); for i:=0 to length(massen)-1 do begin dens:=0; for j:=0 to aX-1 do dens:=dens+impulsIntegral(j,i,msN); dens:=dens*teilchen[i].eigenschaften[tsgMasse]*gitter.dX; gitter.prot.schreibe('n['+inttostr(i+1)+'] = '+floattostr(dens)+' (relative Abweichung: '+floattostr(dens/massen[i]-1)+')',true); end; end; end; function tFelder.impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): double; begin if tlc<0 then begin result:=(emQuellen[emQ]+ort)^; // das ist leicht :-) exit; end; result:=0; case emQ of eqRho: result:=impulsIntegral(ort,tlc,msN); eqJX: result:=impulsIntegral(ort,tlc,msVX); eqJY: result:=impulsIntegral(ort,tlc,msVY); else begin gitter.prot.schreibe('unbekannte EMQuellgröße in Impulsintegral!'); raise exception.create('Fehler in tFelder.impulsIntegral!'); end; end{of case}; result:=result * teilchen[tlc].eigenschaften[tsgLadung]; end; function tFelder.impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): double; var i: longint; begin if tlc<0 then begin result:=0; for i:=0 to length(teilchen)-1 do result:=result+impulsIntegral(ort,i,maF); exit; end; result:=0; case maF of msN: for i:=0 to aP-1 do result:= result + (impulsraum[tlc,false]+i+ort*aP)^ * dP; msPX: for i:=0 to aP-1 do result:= result + (impulsraum[tlc,false]+i+ort*aP)^ * dP * (i-aP/2)*dP; msPXSqr: for i:=0 to aP-1 do result:= result + (impulsraum[tlc,false]+i+ort*aP)^ * dP * sqr((i-aP/2)*dP); msPY: begin for i:=0 to aP-1 do result:= result + (impulsraum[tlc,false]+i+ort*aP)^ * dP; result:=result * (matFelder[tlc,mfPY]+ort)^; end; msVX: begin for i:=0 to aP-1 do result:= result + (impulsraum[tlc,false]+i+ort*aP)^ * dP * (i-aP/2)*dP / sqrt(1 + (sqr((i-aP/2)*dP) + sqr((matFelder[tlc,mfPY]+ort)^))*teilchen[tlc].eigenschaften[tsgIQdrMasse]); result:=result / teilchen[tlc].eigenschaften[tsgMasse]; end; msVY: begin for i:=0 to aP-1 do result:= result + (impulsraum[tlc,false]+i+ort*aP)^ * dP / sqrt(1 + (sqr((i-aP/2)*dP) + sqr((matFelder[tlc,mfPY]+ort)^))*teilchen[tlc].eigenschaften[tsgIQdrMasse]); result:=result * (matFelder[tlc,mfPY]+ort)^ / teilchen[tlc].eigenschaften[tsgMasse]; end; msPXRipple: for i:=1 to aP-1 do result:= result + abs((impulsraum[tlc,false]+i+ort*aP)^-(impulsraum[tlc,false]+i-1+ort*aP)^)*dP; else begin gitter.prot.schreibe('unbekannte Materiefeldgröße in Impulsintegral!'); raise exception.create('Fehler in tFelder.impulsIntegral!'); end; end{of case}; end; {$DEFINE LiKoImplementation} {$INCLUDE linearkombinationen.inc} {$UNDEF LiKoImplementation} procedure tFelder.schreibePhasenraum(tlc: longint; datNamPre, datNamSuf: string); var f: textfile; i,j: longint; typ: byte; begin assignfile(f,datNamPre+datNamSuf); rewrite(f); for i:=0 to aP-1 do begin for j:=0 to aX-1 do begin if j>0 then write(f,#9); write(f,(impulsraum[tlc,false]+i+j*aP)^); end; writeln(f); end; closefile(f); berechnePhasenraumAbleitungen; for typ:=0 to 1 do begin assignfile(f,datNamPre+'_d'+char(ord('x')+(ord('p')-ord('x'))*typ)+datNamSuf); rewrite(f); for i:=0 to aP-1 do begin for j:=0 to aX-1 do begin if j>0 then write(f,#9); write(f,(impulsRaumGradient[tlc,typ]+i+j*aP)^); end; writeln(f); end; closefile(f); end; end; {$IFDEF dichteNichtnegativieren} procedure tFelder.nichtnegativieren; // Dichten nicht negativ machen var i,j: longint; defizit {$IF dichteNichtnegativieren = perSkalieren} ,enfizit {$ENDIF}: double; {$IF dichteNichtnegativieren = perEinzelklau} perm: pTLongintArray; {$ENDIF} begin for i:=0 to length(teilchen)-1 do begin defizit:=0; {$IF dichteNichtnegativieren = perSkalieren} enfizit:=0; {$ENDIF} for j:=0 to aP*aX-1 do if (impulsraum[i,false]+j)^ < 0 then begin defizit:=defizit-(impulsraum[i,false]+j)^; (impulsraum[i,false]+j)^:=0; end {$IF dichteNichtnegativieren = perSkalieren} else enfizit:=enfizit+(impulsraum[i,false]+j)^ {$ENDIF}; gesamtDefizit:=gesamtDefizit+defizit; {$IF dichteNichtnegativieren = perEinzelklau} if defizit>0 then begin perm:=gitter.gibImpulsRaumPermutation; for j:=0 to aP*aX-1 do begin if (impulsraum[i,false]+perm^[j])^>=defizit then begin (impulsraum[i,false]+perm^[j])^:= (impulsraum[i,false]+perm^[j])^-defizit; defizit:=0; break; end else if (impulsraum[i,false]+perm^[j])^>0 then begin defizit:=defizit-(impulsraum[i,false]+perm^[j])^; (impulsraum[i,false]+perm^[j])^:=0; end; end; end; if defizit>0 then begin {$ENDIF} {$IF dichteNichtnegativieren = perSkalieren} if enfizit=0 then begin {$ENDIF} gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, '+floattostr(defizit)+' bleibt übrig!',true); gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, es ist nichts positives mehr übrig!',true); gitter.abbrechen; continue; end; {$IF dichteNichtnegativieren = perSkalieren} enfizit:=1-defizit/enfizit; if defizit>0 then for j:=0 to aP*aX-1 do (impulsraum[i,false]+j)^:=(impulsraum[i,false]+j)^ * enfizit; {$ENDIF} end; end; {$ENDIF} procedure tFelder.berechnePhasenraumAbleitungen; var i,j,k: longint; tmp: double; begin {$IFDEF exzessiveArrayBereichsTests} pruefeArrayEnden('Jemand hat das FFT-Array überschrieben!'); {$ENDIF} for i:=0 to length(ffts)-1 do begin {$IFDEF exzessiveArrayBereichsTests} for j:=0 to aP*(aX div 2 + 1) - 1 do begin (fftTmp+j)^.re:=pi; (fftTmp+j)^.im:=-pi; end; {$ENDIF} fftw_execute(ffts[i,0,false]); // FFT in x {$IFDEF exzessiveArrayBereichsTests} for j:=0 to aP*(aX div 2 + 1) - 1 do if ((fftTmp+j)^.re=pi) or ((fftTmp+j)^.im=-pi) then raise Exception.create('Vorwärts-FFT ('+inttostr(i)+',0) schreibt nicht auf Position '+inttostr(j)+'!'); pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!'); {$ENDIF} for k:=0 to aX div 2 do // Iteration über die Werte einer Transformation for j:=0 to aP-1 do begin // Iteration über die Transformationen tmp:=(fftTmp+j+k*aP)^.re; // *i*k*2*pi; (fftTmp+j+k*aP)^.re:=- 2*pi/dX/aX * k * (fftTmp+j+k*aP)^.im/aX * fourierWindow(k/aX); (fftTmp+j+k*aP)^.im:= 2*pi/dX/aX * k * tmp/aX * fourierWindow(k/aX); end; {$IFDEF exzessiveArrayBereichsTests} for j:=0 to aP*aX - 1 do (impulsRaumGradient[i,0]+j)^:=pi+1; {$ENDIF} fftw_execute(ffts[i,0,true]); {$IFDEF exzessiveArrayBereichsTests} for j:=0 to aP*(aX div 2 + 1) - 1 do if (impulsRaumGradient[i,0]+j)^=pi+1 then raise Exception.create('Rückwärts-FFT ('+inttostr(i)+',0) schreibt nicht auf Position '+inttostr(j)+'!'); pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!'); for j:=0 to aX*(aP div 2 + 1) - 1 do begin (fftTmp+j)^.re:=pi; (fftTmp+j)^.im:=-pi; end; {$ENDIF} fftw_execute(ffts[i,1,false]); // FFT in p {$IFDEF exzessiveArrayBereichsTests} for j:=0 to aX*(aP div 2 + 1) - 1 do if ((fftTmp+j)^.re=pi) or ((fftTmp+j)^.im=-pi) then raise Exception.create('Vorwärts-FFT ('+inttostr(i)+',0) schreibt nicht auf Position '+inttostr(j)+'!'); pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!'); {$ENDIF} for k:=0 to aX-1 do // Iteration über die Transformationen for j:=0 to aP div 2 do begin // Iteration über die Werte einer Transformation tmp:=(fftTmp+j+k*(aP div 2 + 1))^.re; // *i*k*2*pi; (fftTmp+j+k*(aP div 2 + 1))^.re:=- 2*pi/dP/aP * j * (fftTmp+j+k*(aP div 2 + 1))^.im/aP * fourierWindow(j/aP); (fftTmp+j+k*(aP div 2 + 1))^.im:= 2*pi/dP/aP * j * tmp/aP * fourierWindow(j/aP); end; {$IFDEF exzessiveArrayBereichsTests} for j:=0 to aP*aX - 1 do (impulsRaumGradient[i,1]+j)^:=pi+1; {$ENDIF} fftw_execute(ffts[i,1,true]); {$IFDEF exzessiveArrayBereichsTests} for j:=0 to aP*(aX div 2 + 1) - 1 do if (impulsRaumGradient[i,1]+j)^=pi+1 then raise Exception.create('Rückwärts-FFT ('+inttostr(i)+',1) schreibt nicht auf Position '+inttostr(j)+'!'); pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!'); {$ENDIF} end; end; function tFelder.transponiereFuer(tlc: longint): pDouble; {$IFDEF exzessiveArrayBereichsTests} var i: longint; {$ENDIF} begin {$IFDEF exzessiveArrayBereichsTests} for i:=0 to aP*aX-1 do (transpTmp+i)^:=sqrt(2); {$ENDIF} fftw_execute(transps[tlc]); {$IFDEF exzessiveArrayBereichsTests} pruefeArrayEnden('TransporniereFuer('+inttostr(tlc)+') schreibt über das Ende hinaus!'); for i:=0 to aP*aX-1 do if (transpTmp+i)^=sqrt(2) then begin gitter.prot.schreibe('TransporniereFuer('+inttostr(tlc)+') schreibt nicht alle Werte!'); gitter.prot.schreibe('Position '+inttostr(i)+' wurde ausgelassen!'); raise exception.create('Fehler in tFelder.transponiereFuer!'); end; {$ENDIF} result:=transpTmp; end; // tGitter ********************************************************************* constructor tGitter.create(besitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; name: string); var i: longint; begin inherited create; i:=aP; aP:=round(power(2,ceil(ln(aP-0.5)/ln(2)))); deltaP:=deltaP*(i-1)/(aP-1); i:=aX; aX:=round(power(2,ceil(ln(aX+4-0.5)/ln(2)))); // zwei Felder links und rechts extra für Randbedingungen deltaX:=deltaX*(i-1)/(aX-4-1); abbruch:=false; simulation:=besitzer; zeitverfahren:=zv; kvs:=bekannteWerte; prot:=tProtokollant.create(simulation.prot,name); dX:=deltaX; iDX:=1/dX; 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; {$IF dichteNichtnegativieren = perEinzelklau} setlength(impulsRaumPermutationen,100); for i:=0 to length(impulsRaumPermutationen)-1 do impulsRaumPermutationen[i]:=permutation(aX*aP); {$ENDIF} for i:=0 to length(felders)-1 do felders[i]:=tFelder.create(aX,aP,deltaX,deltaP,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.abbrechen; var i: longint; begin if not abbruch then for i:=0 to length(felders[aktuelleFelder].matFelder)-1 do felders[aktuelleFelder].schreibePhasenraum(i,extractfilepath(paramstr(0))+'phasenraum_'+inttostr(i+1),'.dump'); abbruch:=true; end; procedure tGitter.iteriereSchritt(dT: double); {$IFDEF negativeDichteueberwachung} var i,j: longint; pro: tProtokollant; {$ENDIF} begin if simulation.gotSigterm or simulation.gotSigint then begin abbrechen; simulation.gotSigterm:=false; simulation.gotSigint:=false; end; if abbruch then begin t:=t+dT; exit; end; kvs.add('t',t); repeat felders[aktuelleFelder].berechneAbleitungen; // y' = y'(t,y(t)) case zeitverfahren of zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT); zfRungeKuttaDreiAchtel: begin {$INCLUDE rk3_8.inc} end; zfRungeKuttaVier: begin {$INCLUDE rk4.inc} end; zfRungeKuttaZehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk108.txt {$INCLUDE rk108.inc} end; zfRungeKuttaZwoelf: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1210.txt {$INCLUDE rk1210.inc} end; zfRungeKuttaVierzehn: begin // Quelle: http://sce.uhcl.edu/rungekutta/rk1412.txt {$INCLUDE rk1412.inc} end; end{of case}; break; until abbruch; aktuelleFelder:=1-aktuelleFelder; {$IFDEF negativeDichteueberwachung} for i:=0 to length(felders[aktuelleFelder].impulsraum)-1 do for j:=0 to felders[aktuelleFelder].aX*felders[aktuelleFelder].aP-1 do if (felders[aktuelleFelder].impulsraum[i,false]+j)^<0 then begin pro:=tProtokollant.create(prot,'iteriereSchritt'); pro.schreibe('n<0 bei:',true); pro.schreibe(' t = ' +floattostr(t)+' (aktuelleFelder = '+inttostr(aktuelleFelder)+'),',true); pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true); pro.schreibe(' n = ' +floattostr((felders[aktuelleFelder].impulsraum[i,false]+j)^),true); abbrechen; end; {$ENDIF} t:=t+dT; end; procedure tGitter.dumpErhaltungsgroeszen; begin felders[aktuelleFelder].dumpErhaltungsgroeszen; end; {$IF dichteNichtnegativieren = perEinzelklau} function tGitter.gibImpulsRaumPermutation: pTLongintArray; begin result:=@(impulsRaumPermutationen[random(length(impulsRaumPermutationen))]); end; {$ENDIF} // tSimulation ***************************************************************** constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string); var ifile: tMyStringlist; zeitverfahren: tZeitverfahren; s,t,aSuffix,aPrefix: string; dX,breite,dP,pMax: double; i: longint; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; pro: tProtokollant; na: pSigActionRec; ausgabeverzeichnisAufraeumen: boolean; 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); raise exception.create('Fehler in tSimulation.create!'); end; // Standardeinstellungen Bereich 'allgemein' zeitverfahren:=zfRungeKuttaVier; dX:=1e-2; kvs.add('λ',1); kvs.add('T',1); kvs.add('ω',2*pi); kvs.add('dX',dX); kvs.add('qe',1); kvs.add('me',1); dT:=-1; sDT:=-1; sDX:=-1; tEnde:=100; breite:=10.0; fortschrittsAnzeige:=false; ausgabeverzeichnisAufraeumen:=false; gotSigusr1:=false; gotSigterm:=false; gotSigint:=false; dP:=-1; pMax:=3; // Standardeinstellungen Bereich 'ausgaben' aPrefix:=extractfilepath(paramstr(0)); aSuffix:='.dat'; setlength(feldAusgabeDateien,0); setlength(impulsRaumAusgabeDateien,0); fillchar(transpTmp,sizeof(transpTmp),#0); // Transpositions-Puffer initialisieren setlength(transpTmp,0); // Standardeinstellungen Bereich 'lichtVon...' lichter:=tMyStringlist.create(prot,'lichter'); // Standardeinstellungen Bereich 'teilchen...' setlength(teilchen,0); repeat if not ifile.readln(s) then begin prot.schreibe('Unerwartetes Dateiende in '''+inName+'''!',true); raise exception.create('Fehler in tSimulation.create!'); 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); raise exception.create('Fehler in tSimulation.create!'); 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 dX:=exprtofloat(false,s,kvs,nil); kvs.add('dX',dX); continue; end; if startetMit('zeitschritt ',s) then begin dT:=exprtofloat(false,s,kvs,nil); kvs.add('dT',dT); continue; end; if startetMit('impulsschritt ',s) then begin dP:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('maximalimpuls ',s) then begin pMax:=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; if s='Ausgabeverzeichnis aufräumen' then begin ausgabeverzeichnisAufraeumen:=true; continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich allgemein!',true); raise exception.create('Fehler in tSimulation.create!'); 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); raise exception.create('Fehler in tSimulation.create!'); 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(feldAusgabeDateien,length(feldAusgabeDateien)+1); feldAusgabeDateien[length(feldAusgabeDateien)-1]:=tFeldAusgabeDatei.create(t,aPrefix,aSuffix,prot,'feldAusgabeDateien['+inttostr(length(feldAusgabeDateien)-1)+']'); end; continue; end; if startetMit('phasenraum ',s) then begin s:=s+','; while s<>'' do begin t:=erstesArgument(s,','); setlength(impulsRaumAusgabeDateien,length(impulsRaumAusgabeDateien)+1); impulsRaumAusgabeDateien[length(impulsRaumAusgabeDateien)-1]:=tImpulsRaumAusgabeDatei.create(t,aPrefix,aSuffix,prot,'impulsRaumAusgabeDateien['+inttostr(length(impulsRaumAusgabeDateien)-1)+']'); end; continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true); raise exception.create('Fehler in tSimulation.create!'); 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); raise exception.create('Fehler in tSimulation.create!'); 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); raise exception.create('Fehler in tSimulation.create!'); 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); raise exception.create('Fehler in tSimulation.create!'); end; if s='teilchen'+inttostr(length(teilchen))+'Ende' then break; if startetMit('ladung ',s) then begin teilchen[length(teilchen)-1].eigenschaften[tsgLadung]:=-exprToFloat(false,s,kvs,nil); kvs.add('q'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].eigenschaften[tsgLadung]); continue; end; if startetMit('masse ',s) then begin teilchen[length(teilchen)-1].eigenschaften[tsgMasse]:=exprToFloat(false,s,kvs,nil); teilchen[length(teilchen)-1].eigenschaften[tsgIQdrMasse]:=1/sqr(teilchen[length(teilchen)-1].eigenschaften[tsgMasse]); kvs.add('m'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].eigenschaften[tsgMasse]); 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; if startetMit('maximales dLnN/dX',s) then begin teilchen[length(teilchen)-1].diffusion[0]:=exprToFloat(false,s,kvs,nil); continue; end; if startetMit('maximales dLnN/dP',s) then begin teilchen[length(teilchen)-1].diffusion[1]:=exprToFloat(false,s,kvs,nil); continue; end; if startetMit('impulsbreite ',s) then begin teilchen[length(teilchen)-1].breite:=exprToFloat(false,s,kvs,nil); continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich teilchen!',true); raise exception.create('Fehler in tSimulation.create!'); until false; continue; end; prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true); raise exception.create('Fehler in tSimulation.create!'); until false; ifile.free; if (length(feldAusgabeDateien)=0) and (length(impulsRaumAusgabeDateien)=0) then begin prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true); raise exception.create('Fehler in tSimulation.create!'); end; if ausgabeverzeichnisAufraeumen then for i:=0 to length(feldAusgabedateien)-1 do feldAusgabeDateien[i].verzeichnisAufraeumen; if dT<0 then dT:=dX/10; if sDT<0 then sDT:=dT; sDT:=1/round(1/sDT); sT:=sDT/2; if sDX<0 then sDX:=dX; if dP<0 then dP:=dX; for i:=0 to length(teilchen)-1 do begin teilchen[i].diffusion[0]:=teilchen[i].diffusion[0] * dX; teilchen[i].diffusion[1]:=teilchen[i].diffusion[1] * dP; end; 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(dX)); pro.schreibe('dT = '+floattostr(dT)); pro.schreibe('breite = '+floattostr(breite)+' (= '+inttostr(round(breite/dX)+1)+' Punkte)'); pro.schreibe('dP = '+floattostr(dP)); pro.schreibe('pMax = '+floattostr(pMax)+' (= '+inttostr(2*round(pMax/dP)+1)+' Punkte)'); 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(feldAusgabeDateien)>0 then begin pro.schreibe('feldAusgaben:'); for i:=0 to length(feldAusgabeDateien)-1 do begin feldAusgabeDateien[i].zeitAnz:=round(1/sDT); feldAusgabeDateien[i].sDT:=sDT; pro.schreibe(' '+feldAusgabeDateien[i].dump); end; end; if length(impulsRaumAusgabeDateien)>0 then begin pro.schreibe('impulsRaumAusgaben:'); for i:=0 to length(impulsRaumAusgabeDateien)-1 do begin impulsRaumAusgabeDateien[i].zeitAnz:=round(1/sDT); impulsRaumAusgabeDateien[i].sDT:=sDT; pro.schreibe(' '+impulsRaumAusgabeDateien[i].dump); end; end; pro.free; if length(simulationen)=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(SIGUSR1, na, nil) <> 0) or (fPSigaction(SIGTERM, na, nil) <> 0) or (fPSigaction(SIGINT, na, nil) <> 0) then begin pro.schreibe('Fehler: '+inttostr(fpgeterrno)+'.'); raise exception.create('Fehler in tSimulation.create!'); end; dispose(na); end; setlength(simulationen,length(simulationen)+1); simulationen[length(simulationen)-1]:=self; gitter:=tGitter.create(self,round(breite/dX)+1,2*round(pMax/dP)+1,dX,dP,kvs,teilchen,lichter,zeitverfahren,'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 for i:=0 to length(feldAusgabeDateien)-1 do feldAusgabeDateien[i].free; setlength(feldAusgabeDateien,0); gitter.free; prot.free; for i:=length(simulationen)-1 downto 0 do if simulationen[i]=self then begin for j:=i+1 to length(simulationen)-1 do simulationen[j-1]:=simulationen[j]; setlength(simulationen,length(simulationen)-1); end; setlength(transpTmp,0); inherited destroy; end; function tSimulation.iteriereSchritt(start: double; var zeitPhysik,zeitDatei: double): boolean; // noch nicht zu Ende? var i: longint; begin result:=false; if gitter.abbruch or (dT>sDT) then dT:=sDT; zeitPhysik:=zeitPhysik-now; if errorCode<2 then try gitter.iteriereSchritt(dT); except zeitPhysik:=zeitPhysik+now; exit; end; zeitPhysik:=zeitPhysik+now; zeitDatei:=zeitDatei-now; while (gitter.t>=sT) and (sT