diff options
Diffstat (limited to 'Physikunit.pas')
| -rw-r--r-- | Physikunit.pas | 1044 |
1 files changed, 395 insertions, 649 deletions
diff --git a/Physikunit.pas b/Physikunit.pas index 556334d..b49675b 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -13,7 +13,7 @@ unit Physikunit; interface uses - Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix; + Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix, fftw_l; type tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier,zfRungeKuttaZehn,zfRungeKuttaZwoelf,zfRungeKuttaVierzehn); @@ -21,7 +21,7 @@ type tEMFeldGroesze = ( efAX,efAY, efDAXDT,efDAYDT, - efEX,efBY,efBZ + efEX,efBZ ); tEMQuellGroesze = ( eqRho,eqJX,eqJY @@ -37,7 +37,6 @@ const type - tRaumPunkt = class; tFelder = class; tGitter = class; tSimulation = class; @@ -87,87 +86,51 @@ type procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant); end; - { tImpulsPunkt } - - tImpulsPunkt = class(tObject) // Besetzungsdichten verschiedener Teilchenspezies für ein Paar (r;p) - private - raumpunkt: tRaumPunkt; - nachbarn: array[0..1,boolean] of tImpulsPunkt; // x- und px-Nachbarn (jeweil in Richtung - und +) - istNullPkt: boolean; - public - werte: array of array[boolean] of extended; // Besetzungsdichten und deren Ableitungen für verschiedene Teilchenspezies - iMGamma: array of extended; // 1/m/gamma (= v/p) - grad: array of array[0..1] of extended; // Ableitung der Besetzungsdichte nach x und px - constructor create(anzTeilchen: longint); overload; // NullPunkt erstellen - constructor create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); overload; - destructor destroy; override; - {$DEFINE LiKotImpulsPunktHeader} - {$INCLUDE linearkombinationen.inc} - {$UNDEF LiKotImpulsPunktHeader} - function nichtnegativieren: extended; inline; - procedure akkumuliereEMQuellen(var emQuellen: tEMQuellen; pX,dVP: extended); inline; - function gibMirWas(var defizit: extended; entfernung, teilchen,richtung: longint; positiv: boolean): boolean; inline; // ich habe jemanden erreicht - procedure setzeNull; inline; - procedure berechneAbleitungen(pX: extended); inline; - property istNullPunkt: boolean read istNullPkt; - end; - - { tRaumPunkt } - - tRaumPunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Raumgitter und deren Zeitableitungen - private - felder: tFelder; - public - phasenraum: array of tImpulsPunkt; // Phasenraumbesetzungsdichten - emFelder: array[tEMFeldGroesze,boolean] of extended; // EM-Felder und deren Zeitableitungen - // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier! - emQuellen: array[tEMQuellGroesze] of extended; // Quellen für die EM-Felder - matFelder: array of array[tMaterieFeldGroesze,boolean] of extended; // impulsunabhängige Materiefelder getrennt nach Teilchenspezies - lN,rN: tRaumPunkt; // Nachbarpunkte - dP: extended; // Impulsdiskretisierung - aP: longint; // Anzahl der Impulsbins - - constructor create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aP: longint; _dP: extended); - destructor destroy; override; - procedure berechneEMFelder(dX: extended; rechts: boolean); overload; inline; - procedure berechneEMFelder(dX,iDX: extended); overload; inline; - procedure berechneEMQuellen; - procedure berechneEMAbleitungen(iDX: extended); inline; - procedure berechnePhasenraumAbleitungen; inline; - {$DEFINE LiKotRaumPunktHeader} - {$INCLUDE linearkombinationen.inc} - {$UNDEF LiKotRaumPunktHeader} - function nichtnegativieren: extended; inline; - function impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; overload; inline; - function impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended; overload; inline; - procedure initialisiereDichte(teilchen: longint; breite,n: extended); inline; - procedure reflektiereTeilchen(rechts: boolean); inline; - procedure setzeNull; inline; - procedure berechnePY; inline; - end; - { TFelder } +(* + procedure berechneAbleitungen(pX: extended); inline; + function nichtnegativieren: extended; inline; +*) tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen private - teilchen: array of tTeilchenSpezies; - lichters: array[boolean] of tMyStringlist; - massen: array of extended; // Gesamtmassen je Teilchenspezies zur Überwachung der Erhaltungsgrößen - gesamtDefizit: extended; - procedure setzeRaender(iDX: extended); + + // 2d-Arrays (x & px) sind wie folgt indiziert: [iP + iX*aP] + + teilchen: array of tTeilchenSpezies; + lichters: array[boolean] of tMyStringlist; + massen: array of extended; // Gesamtmassen je Teilchenspezies zur Überwachung der Erhaltungsgrößen + gesamtDefizit: extended; + impulsRaumGradient: array of array[0..1] of pExtended; // Ableitung des Impulsraumes nach x und px + iMGamma: array of pExtended; // 1/m/gamma (= v/p) + _aX,_aP: longint; + _dX,_dP,iDX: extended; + ffts: array of array[0..1,boolean] of fftw_plan_extended; // fft-Algorithmen + fftTmp: Pcomplex_extended; // Zwischenspeicher für ffts + procedure initialisiereDichte(ort,tlc: longint; breite,n: extended); inline; public - inhalt: array of tRaumPunkt; + emFelder: array[tEMFeldGroesze,boolean] of pExtended; // EM-Felder und deren Zeitableitungen + // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier! + emQuellen: array[tEMQuellGroesze] of pExtended; // Quellen für die EM-Felder + matFelder: array of array[tMaterieFeldGroesze] of pExtended; // px-unabhängige Felder, nach Teilchenspezies sortiert + impulsraum: array of array[boolean] of pExtended; // px-abhängige Felder, nach Teilchenspezies sortiert, sowie deren Ableitung + gitter: tGitter; - constructor create(groesse: longint; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aP: longint; dP: extended); + constructor create(anzX,anzP: longint; deltaX,deltaP: extended; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); destructor destroy; override; - procedure berechneAbleitungen(dX,iDX: extended); inline; - procedure berechneGradienten(iDX,iDP: extended); - {$DEFINE LiKotFelderHeader} - {$INCLUDE linearkombinationen.inc} - {$UNDEF LiKotFelderHeader} + procedure berechneAbleitungen; procedure setzeNull; inline; procedure dumpErhaltungsgroeszen(prot: tProtokollant); + property aX: longint read _aX; + property aP: longint read _aP; + property dX: extended read _dX; + property dP: extended read _dP; + function impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): extended; overload; inline; + function impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): extended; overload; inline; + {$DEFINE LiKoInterface} + {$INCLUDE linearkombinationen.inc} + {$UNDEF LiKoInterface} end; { tGitter } @@ -186,7 +149,7 @@ type aktuelleFelder: longint; felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern - constructor create(derBesitzer: tSimulation; size: longint; deltaX: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string; aP: longint; dP: extended); + constructor create(derBesitzer: tSimulation; aX,aP: longint; deltaX,deltaP: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); destructor destroy; override; procedure iteriereSchritt(dT: extended); procedure dumpErhaltungsgroeszen; @@ -214,7 +177,7 @@ implementation const emFeldNamen: array[tEMFeldGroesze] of string = ( - 'AX','AY','DAXDT','DAYDT','EX','BY','BZ' + 'AX','AY','DAXDT','DAYDT','EX','BZ' ); emQuellNamen: array[tEMQuellGroesze] of string = ( 'RHO','JX','JY' @@ -399,7 +362,7 @@ begin if sT>=nNum+sDT/2 then schreibeKopf; - cnt:=floor(((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)*gitter.dX+Min(gitter.dX,sDX)/2)/sDX+1); + 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; @@ -423,20 +386,20 @@ begin sX:=cX-Min(gitter.dX,sDX)/2; case wasSpeichern of 0: // em-Feld speichern - for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin + for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin while cX>=sX do begin dec(cnt); - move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFeld,ableitung],(buf+bufPos)^,sizeof(extended)); + move((gitter.felders[gitter.aktuelleFelder].emFelder[emFeld,ableitung]+i)^,(buf+bufPos)^,sizeof(extended)); inc(bufPos,sizeof(extended)); sX:=sX+sDX; end; cX:=cX+gitter.dX; end; 1: // em-Quelle - for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin + for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin while cX>=sX do begin dec(cnt); - val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuelle); + val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,emQuelle); move(val,(buf+bufPos)^,sizeof(extended)); inc(bufPos,sizeof(extended)); sX:=sX+sDX; @@ -444,10 +407,10 @@ begin cX:=cX+gitter.dX; end; 2: // Materiefeld - for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin + for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin while cX>=sX do begin dec(cnt); - val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,matFeld); + val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,matFeld); move(val,(buf+bufPos)^,sizeof(extended)); inc(bufPos,sizeof(extended)); sX:=sX+sDX; @@ -583,365 +546,402 @@ begin end; end; -// tImpulsPunkt **************************************************************** - -constructor tImpulsPunkt.create(anzTeilchen: longint); +(* +function tRaumPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen var - i,j: longint; - b: boolean; + i: longint; begin - inherited create; - istNullPkt:=true; // als NullPunkt erstellen - raumPunkt:=nil; - fillchar(werte,sizeof(werte),#0); - setlength(werte,anzTeilchen); - for i:=0 to length(werte)-1 do - for b:=false to true do - werte[i,b]:=0; - fillchar(grad,sizeof(grad),#0); - setlength(grad,anzTeilchen); - for i:=0 to length(grad)-1 do - for j:=0 to 1 do - grad[i,j]:=0; - fillchar(iMGamma,sizeof(iMGamma),#0); - setlength(iMGamma,length(werte)); - for i:=0 to length(iMGamma)-1 do - iMGamma[i]:=1; - nachbarn[0,false]:=self; - nachbarn[0,true]:=self; - nachbarn[1,false]:=self; - nachbarn[1,true]:=self; + result:=0; + for i:=0 to length(phasenraum)-1 do + result:=result+phasenraum[i].nichtnegativieren; end; -constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); -var - i,j: longint; - b: boolean; - np: tImpulsPunkt; -begin - inherited create; - istNullPkt:=false; - raumPunkt:=rp; - fillchar(werte,sizeof(werte),#0); - setlength(werte,anzTeilchen); - for i:=0 to length(werte)-1 do - for b:=false to true do - werte[i,b]:=0; - fillchar(grad,sizeof(grad),#0); - setlength(grad,anzTeilchen); - for i:=0 to length(grad)-1 do - for j:=0 to 1 do - grad[i,j]:=0; - fillchar(iMGamma,sizeof(iMGamma),#0); - setlength(iMGamma,length(werte)); - for i:=0 to length(iMGamma)-1 do - iMGamma[i]:=0; - nachbarn[0,false]:=linksRaum; - nachbarn[0,true]:=nil; - nachbarn[1,false]:=linksImpuls; - nachbarn[1,true]:=nil; - np:=nil; - for i:=0 to 1 do - if assigned(nachbarn[i,false]) then begin - if nachbarn[i,false].nachbarn[i,true].istNullPunkt then - np:=nachbarn[i,false].nachbarn[i,true]; - nachbarn[i,false].nachbarn[i,true]:=self; - end; - if not assigned(np) then - np:=tImpulsPunkt.create(anzTeilchen); - for i:=0 to 1 do - for b:=false to true do - if not assigned(nachbarn[i,b]) then - nachbarn[i,b]:=np; -end; +*) +// tFelder ********************************************************************* -destructor tImpulsPunkt.destroy; +constructor tFelder.create(anzX,anzP: longint; deltaX,deltaP: extended; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter); var - i: longint; - b: boolean; + i,j: longint; + rechts,abl: boolean; + dens: extended; + emF: tEMFeldGroesze; + emQ: tEMQuellGroesze; + maF: tMaterieFeldGroesze; begin - for i:=0 to 1 do - for b:=false to true do - if assigned(nachbarn[i,b]) and - not nachbarn[i,b].istNullPunkt then - nachbarn[i,b].nachbarn[i,not b]:=nil; - setlength(werte,0); - setlength(iMGamma,0); - inherited destroy; -end; + 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]); -function tImpulsPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen -var - defizit: extended; - i: longint; - pro: tProtokollant; -(* i,richtung,entfernung: longint; - jemandErreicht,vorZurueck: boolean; *) -begin - result:=0; - for i:=0 to length(werte)-1 do begin - if werte[i,false]<0 then begin - defizit:=-werte[i,false]; - result:=result+defizit; - werte[i,false]:=0; (* - entfernung:=1; - jemandErreicht:=true; - while (defizit>0) and jemandErreicht do begin - jemandErreicht:=false; - for richtung:=0 to 1 do - for vorZurueck:=false to true do - if assigned(nachbarn[richtung,vorZurueck]) then - jemandErreicht:=nachbarn[richtung,vorZurueck].gibMirWas(defizit,entfernung,i,richtung,vorZurueck) or jemandErreicht; - - inc(entfernung); - end; - if defizit>0 then begin - pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'WertePunkt.nichtnegativieren'); - pro.schreibe('Ich konnte ein Teilchendefizit der Sorte '+inttostr(i+1)+' nicht auffüllen!',true); - raumpunkt.felder.gitter.abbrechen; - end; *) + 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*sizeof(extended)); + fillchar(impulsRaumGradient[i,j]^,aP*aX*sizeof(extended),#0); end; - if werte[i,false]>1E100 then begin - pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'impulsPunkt.nichtnegativieren'); - pro.schreibe('An einer Stelle im Phasenraum sind bereits mehr als 10^100 Teilchen, ich breche ab (t = '+floattostr(raumpunkt.felder.gitter.t)+' T)!',true); - pro.free; - raumpunkt.felder.gitter.abbrechen; + 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(extended)); + fillchar(iMGamma[i]^,aP*aX*sizeof(extended),#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(extended)); + fillchar(emFelder[emF,abl]^,aX*sizeof(extended),#0); end; end; -end; - -procedure tImpulsPunkt.akkumuliereEMQuellen(var emQuellen: tEMQuellen; pX,dVP: extended); -var - i: longint; - tmp: extended; -begin - for i:=0 to length(werte)-1 do begin - iMGamma[i]:=1/(sqrt(1 + (sqr(pX) + sqr(raumpunkt.matFelder[i,mfPY,false]))*raumpunkt.felder.teilchen[i].eigenschaften[tsgIQdrMasse])*raumpunkt.felder.teilchen[i].eigenschaften[tsgMasse]); - tmp:=werte[i,false] * raumpunkt.felder.teilchen[i].eigenschaften[tsgLadung] * dVP; - emQuellen[eqRho]:=emQuellen[eqRho] + tmp; // * iGamma[i] ??? TODO ... falls, dann aber ohne M! - emQuellen[eqJX]:= emQuellen[eqJX] + tmp * pX * iMGamma[i]; - emQuellen[eqJY]:= emQuellen[eqJY] + tmp * raumpunkt.matFelder[i,mfPY,false] * iMGamma[i]; + for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do begin + fftw_getmem(emQuellen[emQ],aX*sizeof(extended)); + fillchar(emQuellen[emQ]^,aX*sizeof(extended),#0); end; -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(extended)); + fillchar(matFelder[i,maF]^,aX*sizeof(extended),#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*sizeof(extended)); + fillchar(impulsraum[i,abl]^,aX*aP*sizeof(extended),#0); + end; -function tImpulsPunkt.gibMirWas(var defizit: extended; entfernung,teilchen,richtung: longint; positiv: boolean): boolean; // ich habe jemanden erreicht -var - i: longint; - b: boolean; -begin - result:=entfernung<=0; - if result then begin - if (defizit>0) and (werte[teilchen,false]>0) then begin - if defizit>werte[teilchen,false] then begin - defizit:=defizit-werte[teilchen,false]; - werte[teilchen,false]:=0; - end - else begin - werte[teilchen,false]:=werte[teilchen,false]-defizit; - defizit:=0; - end; + 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; - end - else - if defizit>0 then - for i:=richtung to 1 do - for b:=false to i>richtung do - if assigned(nachbarn[i,positiv xor b]) then - result:=nachbarn[i,positiv xor b].gibMirWas(defizit,entfernung-1,teilchen,i,positiv xor b) or result; + 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 *',''); + + fftw_getmem(fftTmp,(aP+1)*(aX+1)*sizeof(complex_extended)); + + fillchar(ffts,sizeof(ffts),#0); + setlength(ffts,length(teilchen)); write('.'); + for i:=0 to length(ffts)-1 do begin +// ffts[i,0,false]:= // Planung der Hintransformationen über x +// fftw_plan_many_dft_r2c(1,@_aX,aP,impulsraum[i,false],nil,aP,1,fftTmp,nil,aP,1,[fftw_preserve_input,fftw_measure]); write('1'); + ffts[i,0,true]:= // Planung der Rücktransformationen über x + fftw_plan_many_dft_c2r(1,@_aX,aP,fftTmp,nil,aP,1,impulsRaumGradient[i,0],nil,aP,1,[fftw_measure]); write('2'); + +// ffts[i,1,false]:= // Planung der Hintransformationen über p +// fftw_plan_many_dft_r2c(1,@_aP,aX,impulsraum[i,false],nil,1,aP,fftTmp,nil,1,aP,[fftw_preserve_input,fftw_measure]); write('3'); +// ffts[i,1,true]:= // Planung der Rücktransformationen über p +// fftw_plan_many_dft_c2r(1,@_aP,aX,fftTmp,nil,1,aP,impulsRaumGradient[i,1],nil,1,aP,[fftw_measure]); write('4'); + end; write('>'); + (* + 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); + *) end; -procedure tImpulsPunkt.setzeNull; +destructor tFelder.destroy; var - i: longint; + i,j: longint; abl: boolean; + emF: tEMFeldGroesze; + emQ: tEMQuellGroesze; + maF: tMaterieFeldGroesze; begin - for i:=0 to length(werte)-1 do 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 - werte[i,abl]:=0; - iMGamma[i]:=0; - end; + 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); + inherited destroy; end; -procedure tImpulsPunkt.berechneAbleitungen(pX: extended); +procedure tFelder.initialisiereDichte(ort,tlc: longint; breite,n: extended); var - i: longint; - tx,tp: extended; + i: longint; + ges: extended; begin - // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f - - for i:=0 to length(werte)-1 do begin - tX:= // der Teil vor "Nabla f" - - pX * iMGamma[i]; - tP:= // der Teil vor "Nabla_p f" - - raumpunkt.felder.teilchen[i].eigenschaften[tsgLadung] * - ( raumpunkt.emFelder[efEX,false] - + raumpunkt.matFelder[i,mfPY,false] * iMGamma[i] * raumpunkt.emFelder[efBZ,false]); - werte[i,true]:= - grad[i,0]*tX + - grad[i,1]*tP; + 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; -// tRaumPunkt ****************************************************************** - -constructor tRaumPunkt.create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aP: longint; _dP: extended); +procedure tFelder.berechneAbleitungen; var - emF: tEMFeldGroesze; - emQ: tEMQuellGroesze; - maF: tMaterieFeldGroesze; - abl: boolean; - i: longint; - xN,pN: tImpulsPunkt; + i,j,k: longint; + emQ: tEMQuellGroesze; + emF: tEMFeldGroesze; + rechts: boolean; + tmp,err:extended; begin - inherited create; - felder:=derBesitzer; + + // 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; + + // 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 - emQuellen[emQ]:=0; - for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do - for abl:=false to true do - emFelder[emF,abl]:=0; - fillchar(matFelder,sizeof(matFelder),#0); - setlength(matFelder,teilchen); + 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 abl:=false to true do - matFelder[i,maF,abl]:=0; - lN:=linkerNachbar; - rN:=nil; - if assigned(lN) then - lN.rN:=self; - aP:=_aP; - dP:=_dP; - fillchar(phasenraum,sizeof(phasenraum),#0); - setlength(phasenraum,aP); - for i:=0 to aP-1 do begin - if assigned(linkerNachbar) then - xN:=linkerNachbar.phasenraum[i] - else - xN:=nil; - if i>0 then - pN:=phasenraum[i-1] - else - pN:=nil; - phasenraum[i]:=tImpulsPunkt.create(self,xN,pN,teilchen); - end; -end; + 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; -destructor tRaumPunkt.destroy; -var - i: longint; -begin - for i:=0 to length(phasenraum)-1 do - phasenraum[i].free; - setlength(phasenraum,0); - if assigned(lN) then - lN.rN:=nil; - if assigned(rN) then - rN.lN:=nil; - inherited destroy; -end; + // E und B aus Rho und A berechnen -procedure tRaumPunkt.berechneEMFelder(dX: extended; rechts: boolean); -begin // dEX/dx = rho ... hier muss integriert werden: // EX = EX(x- deltax) + <rho> * deltax - if rechts then - emFelder[efEX,false]:= - lN.emFelder[efEX,false] + - (lN.emQuellen[eqRho] + emQuellen[eqRho]) * dX/2 - else - emFelder[efEX,false]:= - emQuellen[eqRho] * dX; + 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[efBY,false]:=0; // kann man vllt auch hübscher machen, ich bin aber gerade zu faul dafür - emFelder[efBZ,false]:=0; -end; + 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 -procedure tRaumPunkt.berechneEMFelder(dX,iDX: extended); -begin - // dEX/dx = rho ... hier muss integriert werden: - // EX = EX(x- deltax) + <rho> * deltax - - emFelder[efEX,false]:= - lN.emFelder[efEX,false] + - (lN.emQuellen[eqRho] + emQuellen[eqRho]) * dX/2; + // 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) - // B = rot A - emFelder[efBZ,false]:= - (rN.emFelder[efAY,false]-lN.emFelder[efAY,false])*iDX/2; -end; + 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); -procedure tRaumPunkt.berechneEMQuellen; -var - i: longint; - emQ: tEMQuellGroesze; -begin - for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do - emQuellen[emQ]:=0; - for i:=0 to aP-1 do - phasenraum[i].akkumuliereEMQuellen(emQuellen,(i-aP/2)*dP,dP); -end; + // sonstige Ableitungen des A-Feldes berechnen -procedure tRaumPunkt.berechneEMAbleitungen(iDX: extended); // Zeitableitungen der EM-Potentiale berechnen -begin // d2A/dt2 = Laplace(A) - j ( - dE/dt wird auf dA/dt direkt aufgeschlagen !!! ) - emFelder[efDAXDT,true]:= - ( rN.emFelder[efAX,false] - 2*emFelder[efAX,false] + lN.emFelder[efAX,false] )*sqr(iDX) - - emQuellen[eqJX]; - emFelder[efDAYDT,true]:= - ( rN.emFelder[efAY,false] - 2*emFelder[efAY,false] + lN.emFelder[efAY,false] )*sqr(iDX) - - emQuellen[eqJY]; + 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)^; // dA/dt = dA/dt - E - emFelder[efAX,true]:=emFelder[efDAXDT,false] - emFelder[efEX,false]; - emFelder[efAY,true]:=emFelder[efDAYDT,false]; + 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(extended)); + + // Gradienten der Phasenraumbesetzungsdichte berechnen + + for i:=0 to length(ffts)-1 do begin +// fftw_execute(ffts[i,0,false]); + // *i*k*2*pi; +// fftw_execute(ffts[i,0,true]); + write('A'); + fftw_execute(ffts[i,1,false]); + write('B'); + // *i*k*2*pi; + fftw_execute(ffts[i,1,true]); + write('C'); + end; + + tmp:=0; + err:=0; + for i:=0 to length(impulsraum)-1 do + for j:=0 to aP*aX-1 do begin + tmp:=tmp+abs((impulsraum[i,false]+j)^)+abs((impulsRaumGradient[i,1]+j)^); + err:=err+abs((impulsraum[i,false]+j)^-(impulsRaumGradient[i,1]+j)^); +// tmp:=tmp+abs((impulsraum[i,false]+j)^)+abs((impulsRaumGradient[i,0]+j)^)+abs((impulsRaumGradient[i,1]+j)^); +// err:=err+abs((impulsraum[i,false]+j)^-(impulsRaumGradient[i,0]+j)^) +// +abs((impulsraum[i,false]+j)^-(impulsRaumGradient[i,1]+j)^); + end; + writeln(tmp,' ',err); + halt; + + // 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 tRaumPunkt.berechnePhasenraumAbleitungen; +procedure tFelder.setzeNull; var - i: longint; + i,j,k: longint; + abl: boolean; + emF: tEMFeldGroesze; + emQ: tEMQuellGroesze; + maF: tMaterieFeldGroesze; begin - // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f - for i:=0 to aP-1 do - phasenraum[i].berechneAbleitungen((i-(aP div 2))*dP); + 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; -function tRaumPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen +procedure tFelder.dumpErhaltungsgroeszen(prot: tProtokollant); var - i: longint; + i,j: longint; + dens: extended; begin - result:=0; - for i:=0 to length(phasenraum)-1 do - result:=result+phasenraum[i].nichtnegativieren; + dens:=0; + for i:=0 to length(massen)-1 do + dens:=dens+massen[i]/teilchen[i].eigenschaften[tsgMasse]; + prot.schreibe('Gesamtdefizit: '+floattostr(gesamtDefizit)+' (Anteil '+floattostr(gesamtDefizit/dens)+')',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; + prot.schreibe('n['+inttostr(i+1)+'] = '+floattostr(dens)+' (relative Abweichung: '+floattostr(dens/massen[i]-1)+')',true); + end; end; -function tRaumPunkt.impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; +function tFelder.impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): extended; begin - if teilchen<0 then begin - result:=emQuellen[emQ]; // das ist leicht :-) + if tlc<0 then begin + result:=(emQuellen[emQ]+ort)^; // das ist leicht :-) exit; end; result:=0; case emQ of - eqRho: result:=impulsIntegral(teilchen,msN); - eqJX: result:=impulsIntegral(teilchen,msVX); - eqJY: result:=impulsIntegral(teilchen,msVY); + eqRho: result:=impulsIntegral(ort,tlc,msN); + eqJX: result:=impulsIntegral(ort,tlc,msVX); + eqJY: result:=impulsIntegral(ort,tlc,msVY); end{of case}; - result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung]; + result:=result * teilchen[tlc].eigenschaften[tsgLadung]; end; -function tRaumPunkt.impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended; +function tFelder.impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): extended; var i: longint; begin - if teilchen<0 then begin + if tlc<0 then begin result:=0; - for i:=0 to length(felder.teilchen)-1 do - result:=result+impulsIntegral(i,maF); + for i:=0 to length(teilchen)-1 do + result:=result+impulsIntegral(ort,i,maF); exit; end; @@ -952,308 +952,54 @@ begin for i:=0 to aP-1 do result:= result + - phasenraum[i].werte[teilchen,false]*dP; + (impulsraum[tlc,false]+i+ort*aP)^ * dP; msPX: for i:=0 to aP-1 do result:= result + - phasenraum[i].werte[teilchen,false]*dP * (i-aP/2)*dP; + (impulsraum[tlc,false]+i+ort*aP)^ * dP * (i-aP/2)*dP; msPXSqr: for i:=0 to aP-1 do result:= result + - phasenraum[i].werte[teilchen,false]*dP * sqr((i-aP/2)*dP); + (impulsraum[tlc,false]+i+ort*aP)^ * dP * sqr((i-aP/2)*dP); msPY: begin for i:=0 to aP-1 do result:= result + - phasenraum[i].werte[teilchen,false]*dP; - result:=result * matFelder[teilchen,mfPY,false]; + (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 + - phasenraum[i].werte[teilchen,false]*dP * (i-aP/2)*dP - / sqrt(1 + (sqr((i-aP/2)*dP) + sqr(matFelder[teilchen,mfPY,false]))*felder.teilchen[teilchen].eigenschaften[tsgIQdrMasse]); - result:=result / felder.teilchen[teilchen].eigenschaften[tsgMasse]; + (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 + - phasenraum[i].werte[teilchen,false]*dP - / sqrt(1 + (sqr((i-aP/2)*dP) + sqr(matFelder[teilchen,mfPY,false]))*felder.teilchen[teilchen].eigenschaften[tsgIQdrMasse]); - result:=result * matFelder[teilchen,mfPY,false] / felder.teilchen[teilchen].eigenschaften[tsgMasse]; + (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(phasenraum[i].werte[teilchen,false]-phasenraum[i-1].werte[teilchen,false])*dP; + result + abs((impulsraum[tlc,false]+i+ort*aP)^-(impulsraum[tlc,false]+i-1+ort*aP)^)*dP; end{of case}; end; -procedure tRaumPunkt.initialisiereDichte(teilchen: longint; breite,n: extended); -var - i: longint; - ges: extended; -begin - ges:=0; - for i:=0 to length(phasenraum)-1 do begin - phasenraum[i].werte[teilchen,false]:=power(1/2,sqr((i-aP/2)*dP*2/breite)); - ges:=ges + phasenraum[i].werte[teilchen,false]; - end; - ges:=n/ges/dP; - for i:=0 to length(phasenraum)-1 do - phasenraum[i].werte[teilchen,false]:=phasenraum[i].werte[teilchen,false] * ges; -end; - -procedure tRaumPunkt.reflektiereTeilchen(rechts: boolean); -var - i,j: longint; -begin - for i:=(aP div 2 + 1)*byte(rechts) to (aP div 2)-1 + (aP-(aP div 2))*byte(rechts) do - for j:=0 to length(felder.teilchen)-1 do - if phasenraum[i].werte[j,false]>0 then begin - phasenraum[aP-1-i].werte[j,false]:=phasenraum[aP-1-i].werte[j,false]+phasenraum[i].werte[j,false]; - phasenraum[i].werte[j,false]:=0; - end; -end; - -procedure tRaumPunkt.setzeNull; -var - i: longint; - emF: tEMFeldGroesze; - emQ: tEMQuellGroesze; - maF: tMaterieFeldGroesze; - abl: boolean; -begin - for i:=0 to length(phasenraum)-1 do - phasenraum[i].setzeNull; - for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do - emQuellen[emQ]:=0; - for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do - for abl:=false to true do - emFelder[emF,abl]:=0; - for i:=0 to length(matFelder)-1 do - for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do - for abl:=false to true do - matFelder[i,maF,abl]:=0; -end; - -procedure tRaumPunkt.berechnePY; -var - i: longint; -begin - for i:=0 to length(matFelder)-1 do - matFelder[i,mfPY,false]:=-felder.teilchen[i].eigenschaften[tsgLadung] * emFelder[efAY,false]; -end; - -// linearkombinationsspezifische Methoden von tImpulsPunkt, tRaumPunkt und tFelder - -{$DEFINE LiKoimplementation} +{$DEFINE LiKoImplementation} {$INCLUDE linearkombinationen.inc} -{$UNDEF LiKoimplementation} - -// tFelder ********************************************************************* - -constructor tFelder.create(groesse: longint; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aP: longint; dP: extended); -var - i,j: longint; - rechts: boolean; - dens: extended; -begin - inherited create; - gitter:=parent; - gesamtDefizit:=0; - fillchar(teilchen,sizeof(teilchen),#0); - setlength(teilchen,length(_teilchen)); - for i:=0 to length(teilchen)-1 do - teilchen[i]:=tTeilchenSpezies.create(_teilchen[i]); - fillchar(inhalt,sizeof(inhalt),#0); - setlength(inhalt,groesse); - inhalt[0]:=tRaumPunkt.create(nil,self,length(teilchen),aP,dP); - for i:=1 to length(inhalt)-1 do - inhalt[i]:=tRaumPunkt.create(inhalt[i-1],self,length(teilchen),aP,dP); - fillchar(massen,sizeof(massen),#0); - setlength(massen,length(teilchen)); - for i:=0 to length(massen)-1 do - massen[i]:=0; - for i:=0 to length(inhalt)-1 do - for j:=0 to length(teilchen)-1 do begin - dens:=_teilchen[j].gibDichte(gitter.xl+(i-2)*gitter.dX,parent.kvs); - inhalt[i].initialisiereDichte(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: longint; -begin - setlength(massen,0); - for i:=0 to length(teilchen)-1 do - teilchen[i].free; - setlength(teilchen,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(iDX: extended); -var - emF: tEMFeldGroesze; - rechts: boolean; - i: longint; -begin - for emF:=efAX to efAY do begin // Vakuumrandbedingungen für das A-Feld - inhalt[0].emFelder[emF,true]:= - (inhalt[1].emFelder[emF,false] - - inhalt[0].emFelder[emF,false])*iDX; - inhalt[length(inhalt)-1].emFelder[emF,true]:= - (inhalt[length(inhalt)-2].emFelder[emF,false] - - inhalt[length(inhalt)-1].emFelder[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 - inhalt[(length(inhalt)-1)*byte(rechts)].emFelder[efAY,true]:= - inhalt[(length(inhalt)-1)*byte(rechts)].emFelder[efAY,true] + - exprToFloat(false,lichters[rechts][i],gitter.kvs,nil); -end; - -procedure tFelder.berechneGradienten(iDX,iDP: extended); -var - res,ims: array of extended; - i,j,tlc,len: longint; - tmp: extended; -begin - - // Ableitungen nach PX - - len:=length(inhalt[0].phasenraum); - setlength(res,len); - setlength(ims,len); - for i:=0 to length(inhalt)-1 do - for tlc:=0 to length(teilchen)-1 do begin - if len<>length(inhalt[i].phasenraum) then - raise exception.create('Unterschiedliche Diskretisierung im Phasenraum an unterschiedlichen Orten!'); - for j:=0 to length(inhalt[i].phasenraum)-1 do begin - res[j]:=inhalt[i].phasenraum[j].werte[tlc,false]; - ims[j]:=0; - end; - fft(res,ims,false); - for j:=0 to (len div 2)-1 do begin - tmp:=res[j]; - res[j]:=ims[j]*j; - ims[j]:=-tmp*j; - - tmp:=res[len-1-j]; - res[len-1-j]:=ims[len-1-j]*(-j-1); - ims[len-1-j]:=-tmp*(-j-1); - end; - fft(res,ims,true); - for j:=0 to length(inhalt[i].phasenraum)-1 do - inhalt[i].phasenraum[j].grad[tlc,0]:=res[j]*iDP; - end; - - // Ableitungen nach X - - len:=length(inhalt); - setlength(res,len); - setlength(ims,len); - for i:=0 to length(inhalt[0].phasenraum)-1 do - for tlc:=0 to length(teilchen)-1 do begin - for j:=0 to length(inhalt)-1 do begin - res[j]:=inhalt[j].phasenraum[i].werte[tlc,false]; - ims[j]:=0; - end; - fft(res,ims,false); - for j:=0 to (len div 2)-1 do begin - tmp:=res[j]; - res[j]:=ims[j]*j; - ims[j]:=-tmp*j; - - tmp:=res[len-1-j]; - res[len-1-j]:=ims[len-1-j]*(-j-1); - ims[len-1-j]:=-tmp*(-j-1); - end; - fft(res,ims,true); - for j:=0 to length(inhalt)-1 do - inhalt[j].phasenraum[i].grad[tlc,1]:=res[j]*iDX; - end; -end; - -procedure tFelder.berechneAbleitungen(dX,iDX: extended); -var - i: longint; - b: boolean; -begin - for b:=false to true do // PX ggf. korrigieren - inhalt[(length(inhalt)-1)*byte(b)].reflektiereTeilchen(b); - - for i:=0 to length(inhalt)-1 do // PY berechnen - inhalt[i].berechnePY; - - for i:=0 to length(inhalt)-1 do // rho, j und nebenbei iGamma berechnen - inhalt[i].berechneEMQuellen; - - inhalt[0].berechneEMFelder(dX,false); // E und B aus Rho und A berechnen - for i:=1 to length(inhalt)-2 do - inhalt[i].berechneEMFelder(dX,iDX); - inhalt[length(inhalt)-1].berechneEMFelder(dX,false); - - setzeRaender(iDX); // Randbedingungen für die Ableitungen des A-Felds setzen - - for i:=1 to length(inhalt)-2 do // sonstige Ableitungen des A-Feldes berechnen - inhalt[i].berechneEMAbleitungen(iDX); - - berechneGradienten(iDX,1/inhalt[0].dP); // Gradienten der Phasenraumbesetzungsdichte berechnen - - for i:=1 to length(inhalt)-2 do // Ableitungen der Phasenraumbesetzungsdichten berechnen (Rand wird ausgelassen!) - inhalt[i].berechnePhasenraumAbleitungen; -end; - -procedure tFelder.setzeNull; -var - i: longint; -begin - for i:=0 to length(inhalt)-1 do - inhalt[i].setzeNull; -end; - -procedure tFelder.dumpErhaltungsgroeszen(prot: tProtokollant); -var - i,j: longint; - dens: extended; -begin - dens:=0; - for i:=0 to length(massen)-1 do - dens:=dens+massen[i]/teilchen[i].eigenschaften[tsgMasse]; - prot.schreibe('Gesamtdefizit: '+floattostr(gesamtDefizit)+' (Anteil '+floattostr(gesamtDefizit/dens)+')',true); - for i:=0 to length(massen)-1 do begin - dens:=0; - for j:=0 to length(inhalt)-1 do - dens:=dens+inhalt[j].impulsIntegral(i,msN); - dens:=dens*teilchen[i].eigenschaften[tsgMasse]*gitter.dX; - prot.schreibe('n['+inttostr(i+1)+'] = '+floattostr(dens)+' (relative Abweichung: '+floattostr(dens/massen[i]-1)+')',true); - end; -end; +{$UNDEF LiKoImplementation} // tGitter ********************************************************************* -constructor tGitter.create(derBesitzer: tSimulation; size: longint; deltaX: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string; aP: longint; dP: extended); +constructor tGitter.create(derBesitzer: tSimulation; aX,aP: longint; deltaX,deltaP: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); var i: longint; begin @@ -1261,11 +1007,11 @@ begin i:=aP; aP:=round(power(2,ceil(ln(aP-0.5)/ln(2)))); - dP:=dP*(i-1)/(aP-1); + deltaP:=deltaP*(i-1)/(aP-1); - i:=size; - size:=round(power(2,ceil(ln(size+4-0.5)/ln(2)))); // zwei Felder links und rechts extra für Randbedingungen - deltaX:=deltaX*(i-1)/(size-4-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; besitzer:=derBesitzer; @@ -1289,7 +1035,7 @@ begin xl:=dX/2; for i:=0 to length(felders)-1 do - felders[i]:=tFelder.create(size,teilchen,lichter,self,aP,dP); + felders[i]:=tFelder.create(aX,aP,deltaX,deltaP,teilchen,lichter,self); aktuelleFelder:=0; t:=0; @@ -1328,7 +1074,7 @@ begin kvs.add('t',t); repeat - felders[aktuelleFelder].berechneAbleitungen(dX,iDX); // y' = y'(t,y(t)) + felders[aktuelleFelder].berechneAbleitungen; // y' = y'(t,y(t)) case zeitverfahren of zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt @@ -1749,7 +1495,7 @@ begin setlength(sighupSimulationen,length(sighupSimulationen)+1); sighupSimulationen[length(sighupSimulationen)-1]:=self; - gitter:=tGitter.create(self,round(breite/dX)+1,dX,kvs,teilchen,lichter,zeitverfahren,prot,'gitter',2*round(pMax/dP)+1,dP); + gitter:=tGitter.create(self,round(breite/dX)+1,2*round(pMax/dP)+1,dX,dP,kvs,teilchen,lichter,zeitverfahren,prot,'gitter'); for i:=0 to length(teilchen)-1 do teilchen[i].free; setlength(teilchen,0); |