From bf51ca3a16d648c5da41719b06f0f3a598a37ec3 Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Mon, 10 Aug 2015 15:08:35 +0200 Subject: wieder lauffaehig, im Phasenraum. Scheint alte numerische Probleme zu haben ... --- Physikunit.pas | 602 +++++++++++++++++++++++++++++------------------- Plasmapropagation.lpi | 1 + Plasmapropagation.lps | 193 ++++++++-------- input.epost | 10 +- input.plap | 24 +- linearkombination.inc | 2 +- linearkombinationen.inc | 84 +++---- 7 files changed, 523 insertions(+), 393 deletions(-) diff --git a/Physikunit.pas b/Physikunit.pas index 727f0dc..4861a14 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -26,6 +26,8 @@ type tEMQuellGroesze = ( eqRho,eqJX,eqJY ); + tMaterieFeldGroesze = (mfPY); + tTeilchenSpeziesGroeszen = (tsgMasse,tsgIQdrMasse,tsgLadung); tEMQuellen = array[tEMQuellGroesze] of extended; const @@ -68,12 +70,14 @@ type tTeilchenSpezies = class private - ortsGrenzen: array of extended; - dichteStuecke: array of string; + ortsGrenzen: array of extended; + dichteStuecke: array of string; public - nMin,nMax, // minimale und maximale Massendichte in nc - spezifischeLadung: extended; // q/m in e/me - constructor create; + nMin,nMax,breite: extended; // minimale und maximale Massendichte in nc; Breite der Impulsverteilung + eigenschaften: array[tTeilchenSpeziesGroeszen] of extended; // q in e; m in me + diffusion: array[0..1] of extended; // 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: extended; kvs: tKnownValues): extended; // Massendichte @@ -85,19 +89,20 @@ type tImpulsPunkt = class(tObject) // Besetzungsdichten verschiedener Teilchenspezies für ein Paar (r;p) private raumpunkt: tRaumPunkt; - nachbarn: array[0..2,boolean] of tImpulsPunkt; // x-, px-, py-Nachbarn (jeweil in Richtung - und +) + nachbarn: array[0..1,boolean] of tImpulsPunkt; // x- und px-Nachbarn (jeweil in Richtung - und +) public - werte: array of array[boolean] of extended; // Besetzungsdichten und deren Ableitungen für verschiedene Teilchenspezies - constructor create(rp: tRaumPunkt; linksRaum,linksImpuls,untenImpuls: tImpulsPunkt; anzTeilchen: longint); + werte: array of array[boolean] of extended; // Besetzungsdichten und deren Ableitungen für verschiedene Teilchenspezies + iMGamma: array of extended; // 1/m/gamma (= v/p) + constructor create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); destructor destroy; override; {$DEFINE LiKotImpulsPunktHeader} {$INCLUDE linearkombinationen.inc} {$UNDEF LiKotImpulsPunktHeader} - procedure nichtnegativieren; - procedure akkumuliereEMQuellen(var emQuellen: tEMQuellen; vx,vy,dVP: extended); - function gibMirWas(var defizit: extended; entfernung, teilchen,richtung: longint; positiv: boolean): boolean; // ich habe jemanden erreicht - procedure setzeNull; - procedure berechneAbleitungen(pX,pY,iDX,iDPX,iDPY: extended); + 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,iDX,iDPX: extended); inline; end; { tRaumPunkt } @@ -106,56 +111,55 @@ type 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 - lN,rN: tRaumPunkt; - // Nachbarpunkte - dPX,dPY: extended; - // Impulsdiskretisierung - aPX,aPY: longint; - // Anzahl der Impulsbins - - constructor create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aPX,_aPY: longint; _dPX,_dPY: extended); + 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; - procedure berechneEMFelder(dX,iDX: extended); overload; - procedure berechneEMAbleitungen(iDX: extended); - procedure berechnePhasenraumAbleitungen(iDX: extended); + procedure berechneEMFelder(dX: extended; rechts: boolean); overload; inline; + procedure berechneEMFelder(dX,iDX: extended); overload; inline; + procedure berechneEMQuellen; + procedure berechneEMAbleitungen(iDX: extended); inline; + procedure berechnePhasenraumAbleitungen(iDX: extended); inline; {$DEFINE LiKotRaumPunktHeader} {$INCLUDE linearkombinationen.inc} {$UNDEF LiKotRaumPunktHeader} - procedure nichtnegativieren; - function impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; - procedure initialisiereDichte(teilchen: longint; n: extended); - procedure reflektiereTeilchen(rechts: boolean); - procedure setzeNull; + function nichtnegativieren: extended; inline; + function impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; overload; inline; + function impulsIntegral(teilchen: longint): extended; overload; inline; + procedure initialisiereDichte(teilchen: longint; breite,n: extended); inline; + procedure reflektiereTeilchen(rechts: boolean); inline; + procedure setzeNull; inline; + procedure berechnePY; inline; end; { TFelder } tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen private - matAnz: longint; - spezLadungen, - massen: tExtendedArray; - lichters: array[boolean] of tMyStringlist; + 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); public inhalt: array of tRaumPunkt; gitter: tGitter; - constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aPX,aPY: longint; dPX,dPY: extended); + constructor create(groesse: longint; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aP: longint; dP: extended); destructor destroy; override; - procedure berechneAbleitungen(dX,iDX: extended); + procedure berechneAbleitungen(dX,iDX: extended); inline; {$DEFINE LiKotFelderHeader} {$INCLUDE linearkombinationen.inc} {$UNDEF LiKotFelderHeader} - procedure setzeNull; + procedure setzeNull; inline; + procedure dumpErhaltungsgroeszen(prot: tProtokollant); end; { tGitter } @@ -174,9 +178,10 @@ 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; aPX,aPY: longint; dPX,dPY: extended); + 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); destructor destroy; override; procedure iteriereSchritt(dT: extended); + procedure dumpErhaltungsgroeszen; end; { tSimulation } @@ -349,7 +354,7 @@ procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sT, sDX: extended); var i,cnt: longint; sX,cX,val: extended; begin - if (teilchen>=gitter.felders[gitter.aktuelleFelder].matAnz) then 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); pro.destroyall; halt(1); @@ -424,6 +429,8 @@ end; // tTeilchenSpezies ************************************************************ constructor tTeilchenSpezies.create; +var + egr: tTeilchenSpeziesGroeszen; begin inherited create; fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0); @@ -432,7 +439,26 @@ begin setlength(dichteStuecke,0); nMin:=1E-3; nMax:=nMin; - spezifischeLadung:=0; + 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; end; destructor tTeilchenSpezies.destroy; @@ -446,8 +472,11 @@ procedure tTeilchenSpezies.dump(prot: tProtokollant; prefix: string); var i: longint; begin - prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc'); - prot.schreibe(prefix+'q/m = '+floattostr(-spezifischeLadung)+' e/me'); + prot.schreibe(prefix+'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]); @@ -520,7 +549,7 @@ end; // tImpulsPunkt **************************************************************** -constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls,untenImpuls: tImpulsPunkt; anzTeilchen: longint); +constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); var i: longint; b: boolean; @@ -532,13 +561,15 @@ begin for i:=0 to length(werte)-1 do for b:=false to true do werte[i,b]:=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; - nachbarn[2,false]:=untenImpuls; - nachbarn[2,true]:=nil; - for i:=0 to 2 do + for i:=0 to 1 do if assigned(nachbarn[i,false]) then nachbarn[i,false].nachbarn[i,true]:=self; end; @@ -548,33 +579,36 @@ var i: longint; b: boolean; begin - for i:=0 to 3 do + for i:=0 to 1 do for b:=false to true do if assigned(nachbarn[i,b]) then nachbarn[i,b].nachbarn[i,not b]:=nil; setlength(werte,0); + setlength(iMGamma,0); inherited destroy; end; -procedure tImpulsPunkt.nichtnegativieren; // Dichten nicht negativ machen +function tImpulsPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen var - i,j,entfernung: longint; - defizit: extended; - pro: tProtokollant; - jemandErreicht,b: boolean; + i,richtung,entfernung: longint; + defizit: extended; + pro: tProtokollant; + jemandErreicht,vorZurueck: boolean; begin + result:=0; for i:=0 to length(werte)-1 do 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 j:=0 to 3 do - for b:=false to true do - if assigned(nachbarn[j,b]) then - jemandErreicht:=nachbarn[j,b].gibMirWas(defizit,i,entfernung,j,b) or jemandErreicht; + 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; @@ -586,22 +620,24 @@ begin end; end; -procedure tImpulsPunkt.akkumuliereEMQuellen(var emQuellen: tEMQuellen; vx,vy,dVP: extended); +procedure tImpulsPunkt.akkumuliereEMQuellen(var emQuellen: tEMQuellen; pX,dVP: extended); var i: longint; tmp: extended; begin for i:=0 to length(werte)-1 do begin - tmp:=werte[i,false] * raumpunkt.felder.spezLadungen[i] * dVP; - emQuellen[eqRho]:=emQuellen[eqRho] + tmp; - emQuellen[eqJX]:= emQuellen[eqJX] + tmp * vx; - emQuellen[eqJY]:= emQuellen[eqJY] + tmp * vy; + 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]; end; end; -function tImpulsPunkt.gibMirWas(var defizit: extended; entfernung, teilchen,richtung: longint; positiv: boolean): boolean; // ich habe jemanden erreicht +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 @@ -612,17 +648,16 @@ begin end else begin werte[teilchen,false]:=werte[teilchen,false]-defizit; - defizit:=0;; + defizit:=0; end; end; end else - if defizit>0 then begin - result:=false; - for i:=richtung to 4 do - if assigned(nachbarn[i,positiv]) then - result:=nachbarn[i,positiv].gibMirWas(defizit,entfernung-1,teilchen,i,positiv) or result; - end; + 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; end; procedure tImpulsPunkt.setzeNull; @@ -630,41 +665,61 @@ var i: longint; abl: boolean; begin - for i:=0 to length(werte)-1 do + for i:=0 to length(werte)-1 do begin for abl:=false to true do werte[i,abl]:=0; + iMGamma[i]:=0; + end; end; -procedure tImpulsPunkt.berechneAbleitungen(pX,pY,iDX,iDPX,iDPY: extended); +procedure tImpulsPunkt.berechneAbleitungen(pX,iDX,iDPX: extended); var - igam: extended; - i: longint; + i: longint; + tx,tp: extended; begin - igam:=1/sqrt(1+sqr(pX)+sqr(pY)); - pX:=pX*igam; - pY:=pY*igam; - // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f - for i:=0 to length(werte)-1 do + + 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]:= - - pX * (nachbarn[0,true].werte[i,false] - nachbarn[0,false].werte[i,false]) * iDX/2 - - raumpunkt.felder.spezLadungen[i] * ( - (raumpunkt.emFelder[efEX,false] + pY*raumpunkt.emFelder[efBZ,false]) - * (nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false]) * iDPX/2 - -pX*raumpunkt.emFelder[efBZ,false] - * (nachbarn[2,true].werte[i,false] - nachbarn[2,false].werte[i,false]) * iDPY/2 - ); + tX * ( + ( nachbarn[0,true].werte[i,false] - nachbarn[0,false].werte[i,false]) + - sign(tX)*raumpunkt.felder.teilchen[i].diffusion[0] * ( + nachbarn[0,true].werte[i,false] + + nachbarn[0,false].werte[i,false] - + 2*werte[i,false] + ) + ) * iDX/2 + + tP * ( + ( nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false]) + - sign(tP)*raumpunkt.felder.teilchen[i].diffusion[1] * ( + nachbarn[1,true].werte[i,false] + + nachbarn[1,false].werte[i,false] - + 2*werte[i,false] + ) + )* iDPX/2; + end; + + // Der jeweils zweite Summand entspricht (in etwa) einer thermischen + // Verschmierung und soll die numerische Stabilität bis zu + // d(ln n)/d[xp] * d[XP] = "diffusion" gewährleisten. end; // tRaumPunkt ****************************************************************** -constructor tRaumPunkt.create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aPX,_aPY: longint; _dPX,_dPY: extended); +constructor tRaumPunkt.create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aP: longint; _dP: extended); var - emF: tEMFeldGroesze; - emQ: tEMQuellGroesze; - abl: boolean; - i,j: longint; - xN,pxN,pyN: tImpulsPunkt; + emF: tEMFeldGroesze; + emQ: tEMQuellGroesze; + maF: tMaterieFeldGroesze; + abl: boolean; + i: longint; + xN,pN: tImpulsPunkt; begin inherited create; felder:=derBesitzer; @@ -673,32 +728,31 @@ begin 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 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; - aPX:=_aPX+byte(not odd(_aPX)); - aPY:=_aPY+byte(not odd(_aPY)); - dPX:=_dPX; - dPY:=_dPY; + aP:=_aP+byte(not odd(_aP)); + dP:=_dP; fillchar(phasenraum,sizeof(phasenraum),#0); - setlength(phasenraum,aPX*aPY); - for i:=0 to aPX-1 do - for j:=0 to aPY-1 do begin - if assigned(linkerNachbar) then - xN:=linkerNachbar.phasenraum[i+aPX*j] - else - xN:=nil; - if i>0 then - pxN:=phasenraum[i-1+aPX*j] - else - pxN:=nil; - if j>0 then - pyN:=phasenraum[i+aPX*(j-1)] - else - pyN:=nil; - phasenraum[i+aPX*j]:=tImpulsPunkt.create(self,xN,pxN,pyN,teilchen); - end; + 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; destructor tRaumPunkt.destroy; @@ -747,17 +801,19 @@ begin (rN.emFelder[efAY,false]-lN.emFelder[efAY,false])*iDX/2; end; -procedure tRaumPunkt.berechneEMAbleitungen(iDX: extended); // Zeitableitungen der EM-Potentiale berechnen +procedure tRaumPunkt.berechneEMQuellen; var - i,j: longint; + i: longint; emQ: tEMQuellGroesze; begin for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do emQuellen[emQ]:=0; - for i:=0 to aPX-1 do - for j:=0 to aPY-1 do - phasenraum[i+aPX*j].akkumuliereEMQuellen(emQuellen,(i-aPX/2)*dPX,(j-aPY/2)*dPY,dPX*dPY); + for i:=0 to aP-1 do + phasenraum[i].akkumuliereEMQuellen(emQuellen,(i-aP/2)*dP,dP); +end; +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) @@ -773,25 +829,25 @@ end; procedure tRaumPunkt.berechnePhasenraumAbleitungen(iDX: extended); var - i,j: longint; + i: longint; begin // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f - for i:=1 to aPX-2 do - for j:=1 to aPY-2 do - phasenraum[i+aPX*j].berechneAbleitungen((i-(aPX div 2))*dPX,(j-(aPY div 2))*dPY,iDX,1/dPX,1/dPY); + for i:=1 to aP-2 do // Rand wird ausgelassen! + phasenraum[i].berechneAbleitungen((i-(aP div 2))*dP,iDX,1/dP); end; -procedure tRaumPunkt.nichtnegativieren; // Dichten nicht negativ machen +function tRaumPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen var - i: longint; + i: longint; begin + result:=0; for i:=0 to length(phasenraum)-1 do - phasenraum[i].nichtnegativieren; + result:=result+phasenraum[i].nichtnegativieren; end; function tRaumPunkt.impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; var - i,j: longint; + i: longint; begin if teilchen<0 then begin result:=emQuellen[emQ]; // das ist leicht :-) @@ -802,48 +858,71 @@ begin case emQ of eqRho: - for i:=0 to aPX-1 do - for j:=0 to aPY-1 do - result:= - result + - phasenraum[i+aPX*j].werte[teilchen,false]*dPX*dPY; - eqJX: - for i:=0 to aPX-1 do - for j:=0 to aPY-1 do - result:= - result + - phasenraum[i+aPX*j].werte[teilchen,false]*dPX*dPY * (i-aPX/2)*dPX; - eqJY: - for i:=0 to aPX-1 do - for j:=0 to aPY-1 do - result:= - result + - phasenraum[i+aPX*j].werte[teilchen,false]*dPX*dPY * (j-aPY/2)*dPY; + for i:=0 to aP-1 do + result:= + result + + phasenraum[i].werte[teilchen,false]*dP; + eqJX: 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]; + end; + eqJY: 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]; + end; end{of case}; + result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung]; end; -procedure tRaumPunkt.initialisiereDichte(teilchen: longint; n: extended); +function tRaumPunkt.impulsIntegral(teilchen: longint): extended; var - i: longint; + i,j: longint; begin + result:=0; + + if teilchen<0 then begin + for j:=0 to length(matFelder)-1 do + for i:=0 to aP-1 do + result:=result + phasenraum[i].werte[j,false]*dP; + end + else + for i:=0 to aP-1 do + result:=result + phasenraum[i].werte[teilchen,false]*dP; +end; + +procedure tRaumPunkt.initialisiereDichte(teilchen: longint; breite,n: extended); +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 - if (aPX div 2) + aPX*(aPY div 2) = i then - phasenraum[i].werte[teilchen,false]:=n - else - phasenraum[i].werte[teilchen,false]:=0; + phasenraum[i].werte[teilchen,false]:=phasenraum[i].werte[teilchen,false] * ges; end; procedure tRaumPunkt.reflektiereTeilchen(rechts: boolean); var - i,j,k: longint; + i,j: longint; begin - for i:=(aPX div 2 + 1)*byte(rechts) to (aPX div 2)-1 + (apx-(aPX div 2))*byte(rechts) do - for j:=0 to aPY-1 do - for k:=0 to felder.matAnz-1 do - if phasenraum[i+aPX*j].werte[k,false]>0 then begin - phasenraum[aPX-1-i+aPX*j].werte[k,false]:=phasenraum[i+aPX*j].werte[k,false]; - phasenraum[i+aPX*j].werte[k,false]:=0; - end; + 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; @@ -851,6 +930,7 @@ var i: longint; emF: tEMFeldGroesze; emQ: tEMQuellGroesze; + maF: tMaterieFeldGroesze; abl: boolean; begin for i:=0 to length(phasenraum)-1 do @@ -860,6 +940,18 @@ begin 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 @@ -870,7 +962,7 @@ end; // tFelder ********************************************************************* -constructor tFelder.create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aPX,aPY: longint; dPX,dPY: extended); +constructor tFelder.create(groesse: longint; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aP: longint; dP: extended); var i,j: longint; rechts: boolean; @@ -878,25 +970,25 @@ var begin inherited create; gitter:=parent; - matAnz:=length(teilchen); - fillchar(spezLadungen,sizeof(spezLadungen),#0); - setlength(spezLadungen,matAnz); - for i:=0 to length(spezLadungen)-1 do - spezLadungen[i]:=teilchen[i].spezifischeLadung; + 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+4); // zwei Felder links und rechts extra für Randbedingungen - inhalt[0]:=tRaumPunkt.create(nil,self,length(teilchen),aPX,aPY,dPX,dPY); + 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),aPX,aPY,dPX,dPY); + 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,dens); - massen[j]:=massen[j]+dens*gitter.dX; + 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,''); @@ -912,7 +1004,10 @@ destructor tFelder.destroy; var i: longint; begin - setlength(spezLadungen,0); + 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); @@ -948,19 +1043,26 @@ var i: longint; b: boolean; begin - for b:=false to true do + for b:=false to true do // PX ggf. korrigieren inhalt[(length(inhalt)-1)*byte(b)].reflektiereTeilchen(b); - inhalt[0].berechneEMFelder(dX,false); + 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); + setzeRaender(iDX); // Randbedingungen für die Ableitungen des A-Felds setzen - for i:=1 to length(inhalt)-2 do + for i:=1 to length(inhalt)-2 do // sonstige Ableitungen des A-Feldes berechnen inhalt[i].berechneEMAbleitungen(iDX); - for i:=1 to length(inhalt)-2 do + + for i:=1 to length(inhalt)-2 do // Ableitungen der Phasenraumbesetzungsdichten berechnen (Rand wird ausgelassen!) inhalt[i].berechnePhasenraumAbleitungen(iDX); end; @@ -972,9 +1074,27 @@ begin 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); + 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; + // tGitter ********************************************************************* -constructor tGitter.create(derBesitzer: tSimulation; size: longint; deltaX: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string; aPX,aPY: longint; dPX,dPY: extended); +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); var i: longint; begin @@ -1001,7 +1121,7 @@ begin xl:=dX/2; for i:=0 to length(felders)-1 do - felders[i]:=tFelder.create(size,teilchen,lichter,self,aPX,aPY,dPX,dPY); + felders[i]:=tFelder.create(size,teilchen,lichter,self,aP,dP); aktuelleFelder:=0; t:=0; @@ -1111,20 +1231,25 @@ begin t:=t+dT; end; +procedure tGitter.dumpErhaltungsgroeszen; +begin + felders[aktuelleFelder].dumpErhaltungsgroeszen(prot); +end; + // tSimulation ***************************************************************** constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string); var - ifile: tMyStringlist; - zeitverfahren: tZeitverfahren; - s,t,aSuffix,aPrefix: string; - deltaX,breite,dPX,dPY,pXMax,pYMax: extended; - i: longint; - kvs: tKnownValues; - teilchen: array of tTeilchenSpezies; - lichter: tMyStringlist; - pro: tProtokollant; - na: pSigActionRec; + ifile: tMyStringlist; + zeitverfahren: tZeitverfahren; + s,t,aSuffix,aPrefix: string; + dX,breite,dP,pMax: extended; + i: longint; + kvs: tKnownValues; + teilchen: array of tTeilchenSpezies; + lichter: tMyStringlist; + pro: tProtokollant; + na: pSigActionRec; begin inherited create; prot:=tProtokollant.create(protokollant,name); @@ -1140,12 +1265,12 @@ begin // Standardeinstellungen Bereich 'allgemein' zeitverfahren:=zfRungeKuttaVier; - deltaX:=1e-2; + dX:=1e-2; kvs.add('λ',1); kvs.add('T',1); kvs.add('ω',2*pi); - kvs.add('dX',deltaX); - kvs.add('q',1); + kvs.add('dX',dX); + kvs.add('qe',1); kvs.add('me',1); dT:=-1; sDT:=-1; @@ -1154,10 +1279,8 @@ begin breite:=10.0; fortschrittsAnzeige:=false; gotSighup:=false; - dPX:=-1; - dPY:=-1; - pXMax:=3; - pYMax:=3; + dP:=-1; + pMax:=3; // Standardeinstellungen Bereich 'ausgaben' aPrefix:=extractfilepath(paramstr(0)); @@ -1213,8 +1336,8 @@ begin continue; end; if startetMit('ortsschritt ',s) then begin - deltaX:=exprtofloat(false,s,kvs,nil); - kvs.add('dX',deltaX); + dX:=exprtofloat(false,s,kvs,nil); + kvs.add('dX',dX); continue; end; if startetMit('zeitschritt ',s) then begin @@ -1223,29 +1346,11 @@ begin continue; end; if startetMit('impulsschritt ',s) then begin - if startetMit('x ',s) then begin - dPX:=exprtofloat(false,s,kvs,nil); - continue; - end; - if startetMit('y ',s) then begin - dPY:=exprtofloat(false,s,kvs,nil); - continue; - end; - dPX:=exprtofloat(false,s,kvs,nil); - dPY:=dPX; + dP:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('maximalimpuls ',s) then begin - if startetMit('x ',s) then begin - pXMax:=exprtofloat(false,s,kvs,nil); - continue; - end; - if startetMit('y ',s) then begin - pYMax:=exprtofloat(false,s,kvs,nil); - continue; - end; - pXMax:=exprtofloat(false,s,kvs,nil); - pYMax:=pXMax; + pMax:=exprtofloat(false,s,kvs,nil); continue; end; if startetMit('zeit ',s) then begin @@ -1340,9 +1445,15 @@ begin halt(1); end; if s='teilchen'+inttostr(length(teilchen))+'Ende' then break; - if startetMit('spezifische Ladung ',s) then begin - teilchen[length(teilchen)-1].spezifischeLadung:=-exprToFloat(false,s,kvs,nil); - kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung); + 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 @@ -1359,6 +1470,18 @@ 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); prot.destroyall; @@ -1381,18 +1504,22 @@ begin end; if dT<0 then - dT:=deltaX/10; + dT:=dX/10; if sDT<0 then sDT:=dT; sDT:=1/round(1/sDT); sT:=sDT/2; if sDX<0 then - sDX:=deltaX; + 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; - if dPX<0 then - dPX:=deltaX; - if dPY<0 then - dPY:=sqrt(deltaX); pro:=tProtokollant.create(prot,'create'); case zeitverfahren of @@ -1411,13 +1538,11 @@ begin else pro.schreibe('Iteration mittels unbekanntem Verfahren'); end{of case}; - pro.schreibe('dX = '+floattostr(deltaX)); - pro.schreibe('dT = '+floattostr(dT)); - pro.schreibe('breite = '+floattostr(breite)+' (= '+inttostr(round(breite/deltaX)+1)+' Punkte)'); - pro.schreibe('dPX = '+floattostr(dPX)); - pro.schreibe('pXMax = '+floattostr(pXMax)); - pro.schreibe('dPY = '+floattostr(dPY)); - pro.schreibe('pYMax = '+floattostr(pYMax)); + 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 @@ -1456,7 +1581,7 @@ begin setlength(sighupSimulationen,length(sighupSimulationen)+1); sighupSimulationen[length(sighupSimulationen)-1]:=self; - gitter:=tGitter.create(self,round(breite/deltaX)+1,deltaX,kvs,teilchen,lichter,zeitverfahren,prot,'gitter',2*round(pXMax/dPX)+1,2*round(pYMax/dPY)+1,dPX,dPY); + gitter:=tGitter.create(self,round(breite/dX)+1,dX,kvs,teilchen,lichter,zeitverfahren,prot,'gitter',2*round(pMax/dP)+1,dP); for i:=0 to length(teilchen)-1 do teilchen[i].free; setlength(teilchen,0); @@ -1509,6 +1634,7 @@ begin prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true); prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true); prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true); + gitter.dumpErhaltungsgroeszen; end; result:=gitter.t + diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps index adb0f8b..d47f765 100644 --- a/Plasmapropagation.lps +++ b/Plasmapropagation.lps @@ -7,8 +7,7 @@ - - + @@ -17,12 +16,12 @@ - - - - + + + + - + @@ -31,246 +30,244 @@ - + - - - + + + - - - - - + + + + + + + + + + + + - + - - + + - - + + - - - - - - - - + + + + + + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - - - - + - + - - - + - + - - - + + + - + - - + + + - + - - + + - - - - - - - + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + diff --git a/input.epost b/input.epost index c6f8d9a..96ec796 100644 --- a/input.epost +++ b/input.epost @@ -26,8 +26,8 @@ Threadanzahl: 11 parallel lesen -!setze $symFelder: DPSIDX1 DPHIDX AY DAYDT -!setze $asyFelder: N1 +!setze $symFelder: JX1 AY DAYDT # EX +!setze $asyFelder: RHO1 !Schleife: $Feld: $asyFelder $symFelder @@ -41,10 +41,10 @@ Ende !Schleifenende -!setze $symBerFelder: DPSIDX1_N1 +!setze $symBerFelder: VX1 -Multipliziere DPSIDX1 mal N1 - Name: DPSIDX1_N1 +Teile JX1 durch RHO1 + Name: VX1 Ende !Schleife: $symmetrie: symmetrisch asymmetrisch diff --git a/input.plap b/input.plap index 8e6fe6c..a6905dd 100644 --- a/input.plap +++ b/input.plap @@ -7,11 +7,11 @@ allgemein runge-Kutta-3/8 # runge-Kutta-4 # euler-Vorwärts - ortsschritt 10^-2 * λ - zeitschritt 10^-2 * T -# minZeitschritt 10^-2 * T + ortsschritt 5*10^-3 * λ + zeitschritt 5*10^-3 * T + maximalimpuls 10 + impulsschritt 5 * 10^-3 zeit 40 * T -# diffusionsterm 2 # max. p/Lp !setze $breite: (5 * λ) breite $breite mit Fortschrittsanzeige @@ -20,7 +20,7 @@ allgemeinEnde ausgaben prefix /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/ suffix _test.dat - felder AX,AY,dAYDT,EX,BY,BZ,Rho1,Rho2 + felder AX,AY,dAYDT,EX,BY,BZ,Rho1,Rho2,JX1 ausgabenEnde !setze $tFwhm: (2.5 * T) @@ -31,9 +31,12 @@ licht von links 2 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * ω * cos(ω*t) # Zei !setze $IonenMassenFaktor: (1836.15267245 + 1838.68366158) teilchen1 - spezifische Ladung -q/me + ladung -qe + masse me maximaldichte 10 minimaldichte 0 # 10^-2 + maximales dLnN/dX 10 + maximales dLnN/dP 10 !setze $profilbreite: (4 * λ) !setze $randbreite: (0.1 * λ) verteilung stückweise @@ -50,9 +53,12 @@ teilchen1 teilchen1Ende teilchen2 - spezifische Ladung q/$IonenMassenFaktor/me - maximaldichte nmax1*$IonenMassenFaktor - minimaldichte nmin1*$IonenMassenFaktor + ladung qe + masse $IonenMassenFaktor * me + maximaldichte nmax1 + minimaldichte nmin1 + maximales dLnN/dX 10 + maximales dLnN/dP 10 # verteilung wie teilchen1 verteilung stückweise 0 diff --git a/linearkombination.inc b/linearkombination.inc index 0f74497..85ed695 100644 --- a/linearkombination.inc +++ b/linearkombination.inc @@ -320,6 +320,6 @@ begin {$IFDEF lkA32},fak32 {$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}); for i:=0 to length(inhalt)-1 do - inhalt[i].nichtnegativieren; + gesamtDefizit:=gesamtDefizit+inhalt[i].nichtnegativieren; end; diff --git a/linearkombinationen.inc b/linearkombinationen.inc index 724c695..6717141 100644 --- a/linearkombinationen.inc +++ b/linearkombinationen.inc @@ -2,51 +2,51 @@ // interface ******************************************************************* {$IFDEF LiKotImpulsPunktHeader} -procedure liKo(in1,in2: tImpulsPunkt; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt -procedure liKo(in1,in2,in3: tImpulsPunkt; fak2,fak3: extended); overload; -procedure liKo(in1,in2,in3,in4: tImpulsPunkt; fak2,fak3,fak4: extended); overload; -procedure liKo(in1,in2,in3,in4,in5: tImpulsPunkt; fak2,fak3,fak4,fak5: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload; +procedure liKo(in1,in2: tImpulsPunkt; fak2: extended); overload; inline; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt +procedure liKo(in1,in2,in3: tImpulsPunkt; fak2,fak3: extended); overload; inline; +procedure liKo(in1,in2,in3,in4: tImpulsPunkt; fak2,fak3,fak4: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5: tImpulsPunkt; fak2,fak3,fak4,fak5: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tImpulsPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload; inline; {$ENDIF} {$IFDEF LiKotRaumPunktHeader} -procedure liKo(in1,in2: tRaumPunkt; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt -procedure liKo(in1,in2,in3: tRaumPunkt; fak2,fak3: extended); overload; -procedure liKo(in1,in2,in3,in4: tRaumPunkt; fak2,fak3,fak4: extended); overload; -procedure liKo(in1,in2,in3,in4,in5: tRaumPunkt; fak2,fak3,fak4,fak5: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6: tRaumPunkt; fak2,fak3,fak4,fak5,fak6: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload; -procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload; +procedure liKo(in1,in2: tRaumPunkt; fak2: extended); overload; inline; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt +procedure liKo(in1,in2,in3: tRaumPunkt; fak2,fak3: extended); overload; inline; +procedure liKo(in1,in2,in3,in4: tRaumPunkt; fak2,fak3,fak4: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5: tRaumPunkt; fak2,fak3,fak4,fak5: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6: tRaumPunkt; fak2,fak3,fak4,fak5,fak6: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload; inline; +procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tRaumPunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload; inline; {$ENDIF} {$IFDEF LiKotFelderHeader} -- cgit v1.2.3-70-g09d2