From 924aa698dcd4fb33c925630aeff7b1320e445ab0 Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Thu, 6 Aug 2015 15:08:42 +0200 Subject: lauffaehig, Materie haut aber ab. --- Physikunit.pas | 233 ++++++++++++++++++++++++++++++++++++++------------ Plasmapropagation.lps | 134 +++++++++++++++-------------- input.epost | 42 ++++++--- input.plap | 33 ++++--- linearkombination.inc | 6 +- 5 files changed, 304 insertions(+), 144 deletions(-) diff --git a/Physikunit.pas b/Physikunit.pas index 142ee4b..ab99170 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -8,6 +8,7 @@ unit Physikunit; {$DEFINE Zeitschrittueberwachung} {$DEFINE Dichteueberwachung} +{ $DEFINE negativeDichteueberwachung} interface @@ -25,11 +26,20 @@ type tMaterieFeldInhalt = ( mfN,mfDPsiDX, // Teilchendichte, Geschwindigkeitsdichte mfP,mfPX,mfPY,mfPZ, // Impuls - mfGamma,mfIGamma // rel. Gamma-Faktor + mfGamma,mfIGamma // rel. Gamma-Faktor ); +const + erstesEMFMitAbleitung = efAX; + letztesEMFMitAbleitung = efDAZDT; + erstesMatFMitAbleitung = mfN; + letztesMatFMitAbleitung = mfDPsiDX; + +type + tGitter = class; tFelder = class; + tSimulation = class; { tAusgabeDatei } @@ -81,15 +91,16 @@ type // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen emWerte: array[tEMFeldInhalt,boolean] of extended; // EM-Felder und deren Zeitableitungen - // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! + // A, dPhiDX, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! lN,rN: tWertePunkt; constructor create(linkerNachbar: tWertePunkt; derBesitzer: tFelder); destructor destroy; override; procedure berechneGammaUndP; + procedure berechneDPhiDX(dX: extended); procedure berechneNAbleitung(iDX,pDNMax: extended); overload; procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload; - procedure berechneAbleitungen(dX,iDX: extended); + procedure berechneAbleitungen(iDX: extended); procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt procedure liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); overload; procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload; @@ -112,6 +123,8 @@ type 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: tWertePunkt; 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: tWertePunkt; 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; function maxDT: extended; + function dump(mitNachbarn: boolean): string; + procedure nichtnegativieren; end; { TFelder } @@ -119,7 +132,8 @@ type tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen private matAnz: longint; - spezLadungen: tExtendedArray; + spezLadungen, + massen: tExtendedArray; lichters: array[boolean] of tMyStringlist; procedure setzeRaender(iDX: extended); public @@ -157,6 +171,7 @@ type tGitter = class(tObject) private + besitzer: tSimulation; prot: tProtokollant; zeitverfahren: tZeitverfahren; pDNMaxDynamisch,abbruch: boolean; @@ -170,7 +185,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(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); + constructor create(derBesitzer: tSimulation; size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); destructor destroy; override; procedure iteriereSchritt(var dT: extended; dTMin: extended); function dumpErhaltungsgroessen: boolean; @@ -429,7 +444,7 @@ var i: longint; begin prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc'); - prot.schreibe(prefix+'q/m = '+floattostr(spezifischeLadung)+' e/me'); + prot.schreibe(prefix+'q/m = '+floattostr(-spezifischeLadung)+' e/me'); prot.schreibe(prefix+'n(x):'); for i:=0 to length(dichteStuecke)-1 do begin prot.schreibe(prefix+' '+dichteStuecke[i]); @@ -557,6 +572,29 @@ begin end; end; +procedure tWertePunkt.berechneDPhiDX(dX: extended); +var + i: longint; +begin + // d2Phi/dx2 = rho ... hier muss integriert werden: + // dPhi/dx = dPhi/dx(x- deltax) + * deltax + + if assigned(lN) then begin + emWerte[efDPhiDX,false]:=lN.emWerte[efDPhiDX,false]; + for i:=0 to length(matWerte)-1 do + emWerte[efDPhiDX,false]:= + emWerte[efDPhiDX,false] + + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,false] + matWerte[i,mfN,false]) * dX/2; + end + else begin + emWerte[efDPhiDX,false]:=0; + for i:=0 to length(matWerte)-1 do + emWerte[efDPhiDX,false]:= + emWerte[efDPhiDX,false] + + besitzer.spezLadungen[i] * matWerte[i,mfN,false] * dX/2; + end; +end; + procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended); // Zeitableitung von n berechnen var i: longint; @@ -602,7 +640,7 @@ begin end end; -procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen +procedure tWertePunkt.berechneAbleitungen(iDX: extended); // Zeitableitungen berechnen var i: longint; begin @@ -610,16 +648,7 @@ begin for i:=0 to length(matWerte)-1 do matWerte[i,mfDPsiDX,true]:= besitzer.spezLadungen[i] * emWerte[efDPhiDX,false] - - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2; - - // d3Phi/dx2dt = dRho/dt ... hier muss integriert werden: - // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + * deltax - emWerte[efDPhiDX,true]:= - lN.emWerte[efDPhiDX,true]; - for i:=0 to length(matWerte)-1 do - emWerte[efDPhiDX,true]:= - emWerte[efDPhiDX,true] - + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,true] + matWerte[i,mfN,true])*dX/2; +;// - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2; // d2A/dt2 = Laplace(A) - Nabla(phi) ... emWerte[efDAXDT,true]:= @@ -659,6 +688,73 @@ begin result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]); end; +function tWertePunkt.dump(mitNachbarn: boolean): string; +var + i: longint; +begin + if mitNachbarn and assigned(lN) then result:=lN.dump(false)+' << ' + else result:=''; + + for i:=0 to length(matWerte)-1 do + result:= + result+ + 'N['+ inttostr(i)+'] = '+floattostr(matWerte[i,mfN,false])+' '+ + 'N''['+inttostr(i)+'] = '+floattostr(matWerte[i,mfN,true])+' '; + + result:=result+'maxDT: '+floattostr(maxDT); + + if mitNachbarn and assigned(rN) then result:=result+' >> '+rN.dump(false); +end; + +procedure tWertePunkt.nichtnegativieren; // Dichten nicht negativ machen +var + i: longint; + defizit: extended; + li,re: tWertePunkt; + pro: tProtokollant; +begin + for i:=0 to length(matWerte)-1 do + if matWerte[i,mfN,false]<0 then begin + defizit:=-matWerte[i,mfN,false]; + matWerte[i,mfN,false]:=0; + li:=self.lN; + re:=self.rN; + repeat + if assigned(li) then begin + if li.matWerte[i,mfN,false]>0 then begin + if defizit>=li.matWerte[i,mfN,false] then begin + defizit:=defizit-li.matWerte[i,mfN,false]; + li.matWerte[i,mfN,false]:=0; + end + else begin + li.matWerte[i,mfN,false]:=li.matWerte[i,mfN,false]-defizit; + defizit:=0;; + end; + end; + li:=li.lN; + end; + if assigned(re) then begin + if re.matWerte[i,mfN,false]>0 then begin + if defizit>=re.matWerte[i,mfN,false] then begin + defizit:=defizit-re.matWerte[i,mfN,false]; + re.matWerte[i,mfN,false]:=0; + end + else begin + re.matWerte[i,mfN,false]:=re.matWerte[i,mfN,false]-defizit; + defizit:=0;; + end; + end; + re:=re.rN; + end; + until (defizit<=0) or not (assigned(li) or assigned(re)); + if defizit>0 then begin + pro:=tProtokollant.create(besitzer.par.prot,'WertePunkt.nichtnegativieren'); + pro.schreibe('Ich konnte ein Teilchendefizit nicht auffüllen!',true); + besitzer.par.abbrechen; + end; + end; +end; + { linearkombinationsspezifische Methoden von tWertePunkt und tFelder } {$INCLUDE linearkombination.inc} @@ -722,10 +818,15 @@ begin inhalt[0]:=tWertePunkt.create(nil,self); for i:=1 to length(inhalt)-1 do inhalt[i]:=tWertePunkt.create(inhalt[i-1],self); + 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 inhalt[i].matWerte[j,mfN,false]:=teilchen[j].gibDichte(par.xl+(i-2)*par.dX,parent.kvs); inhalt[i].matWerte[j,mfDPsiDX,false]:=0; + massen[j]:=massen[j]+inhalt[i].matWerte[j,mfN,false]*par.dX; end; for rechts:=false to true do begin lichters[rechts]:=tMyStringlist.create(nil,''); @@ -778,6 +879,8 @@ var begin for i:=0 to length(inhalt)-1 do inhalt[i].berechneGammaUndP; + for i:=0 to length(inhalt)-1 do + inhalt[i].berechneDPhiDX(dX); for i:=0 to matAnz-1 do begin // Teilchen werden reflektiert inhalt[1].matWerte[i,mfPX,false]:= abs(inhalt[1].matWerte[i,mfPX,false]); @@ -794,7 +897,7 @@ begin setzeRaender(iDX); for i:=1 to length(inhalt)-2 do - inhalt[i].berechneAbleitungen(dX,iDX); + inhalt[i].berechneAbleitungen(iDX); end; function tFelder.maxDT: extended; @@ -808,12 +911,13 @@ end; { tGitter } -constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); +constructor tGitter.create(derBesitzer: tSimulation; size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string); var i: longint; begin inherited create; abbruch:=false; + besitzer:=derBesitzer; zeitverfahren:=zv; kvs:=bekannteWerte; prot:=tProtokollant.create(protokollant,name); @@ -859,17 +963,21 @@ end; procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant); var - i,j: longint; - curMax: extended; + i,j,k: longint; + curMax,tmp: extended; begin curMax:=0; + k:=0; for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do - for j:=0 to felders[aktuelleFelder].matAnz-1 do - curMax:= - max(curMax, - abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]* - (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/ - max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-100))); + for j:=0 to felders[aktuelleFelder].matAnz-1 do begin + tmp:=abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]* + (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/ + max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-100)); + if tmp>curMax then begin + curMax:=tmp; + k:=i; + end; + end; if curMax > pDNMax then begin if pDNMaxDynamisch then begin pDNMax:=2*curMax; @@ -878,10 +986,10 @@ begin else if pDNMax>0 then begin pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+ - 'T), die numerische Stabilität ist nicht mehr gewährleistet!',true); + 'T, x='+floattostr(xL+k*dX)+'), die numerische Stabilität ist nicht mehr gewährleistet, Simulation wird abgebrochen!',true); pro.schreibe(' Lösung: größeren Diffusionsterm wählen',true); pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern',true); - pDNMax:=0; // Term komplett abschalten + abbrechen; end; end; end; @@ -892,18 +1000,17 @@ begin // das maximale dT, welches bei den momentanen Ableitungen nicht zu // unphysikalischen Effekten (z.B. negativen Dichten) führt - result:=dTMax1.1*dTMin); if result then begin // beabsichtigter Zeitschritt ist zu groß, dT:=dTMax/4; // dann machen wir ihn doch kleiner! + if dT '+floattostr(dT)+' (t = '+floattostr(t)+')'); - if dT '+floattostr(dT)+' (t = '+floattostr(t)+')'); + prot.spuelen; + end + {$ENDIF}; end; end; @@ -928,11 +1035,8 @@ begin end; procedure tGitter.iteriereSchritt(var dT: extended; dTMin: extended); -var i: longint; - {$IFDEF Zeitschrittueberwachung} - j: longint; - {$ENDIF} - {$IFDEF Dichteueberwachung} +var {$IFDEF negativeDichteueberwachung} + i,j: longint; pro: tProtokollant; {$ENDIF} dTMax: extended; @@ -973,17 +1077,18 @@ begin break; until abbruch; - prot.schreibe('dT = '+floattostr(dT)); + aktuelleFelder:=1-aktuelleFelder; - {$IFDEF Dichteueberwachung} + {$IFDEF negativeDichteueberwachung} for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do for j:=0 to felders[aktuelleFelder].matAnz-1 do - if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then + if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then begin + pro:=tProtokollant.create(prot,'iteriereSchritt'); + pro.schreibe('n<0 bei:',true); + pro.schreibe(' t = ' +floattostr(t)+',',true); + pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true); + pro.schreibe(' alte Werte:',true); with felders[1-aktuelleFelder] do begin - pro:=tProtokollant.create(prot,'iteriereSchritt'); - pro.schreibe('n<0 bei:',true); - pro.schreibe(' t = ' +floattostr(t)+',',true); - pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true); pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* inhalt[i-1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true); @@ -996,9 +1101,25 @@ begin inhalt[i+1].matWerte[j,mfIGamma,false])+',',true); pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true); pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); - pro.destroyall; - halt(1); end; + pro.schreibe(' neue Werte:',true); + with felders[aktuelleFelder] do begin + pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]* + inhalt[i-1].matWerte[j,mfIGamma,false])+',',true); + pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true); + pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true); + pro.schreibe(' v = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]* + inhalt[i].matWerte[j,mfIGamma,false])+',',true); + pro.schreibe(' n = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true); + pro.schreibe(' n'' = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true); + pro.schreibe(' v_r = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]* + inhalt[i+1].matWerte[j,mfIGamma,false])+',',true); + pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true); + pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true); + end; + pro.free; + abbrechen; + end; {$ENDIF} t:=t+dT; dT:=max(dT,dTMax/100); @@ -1025,9 +1146,9 @@ begin ns[i]:=ns[i]*dX; s:=s+ ' n['+inttostr(i)+']='+floattostr(ns[i]); {$IFDEF Dichteueberwachung} - if ns[i]>1000 then begin + if (not abbruch) and (abs(ns[i]-felders[aktuelleFelder].massen[i])>0.5) then begin result:=false; - pro.schreibe(' n['+inttostr(i)+'] > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true); + pro.schreibe(' n['+inttostr(i)+'] ('+floattostr(ns[i])+') unterscheidet sich stark vom Anfangswert ('+floattostr(felders[aktuelleFelder].massen[i])+'), es scheinen sehr viele Teilchen entstanden oder verloren gegangen zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true); end; {$ENDIF} end; @@ -1276,7 +1397,7 @@ begin 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); + teilchen[length(teilchen)-1].spezifischeLadung:=-exprToFloat(false,s,kvs,nil); kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung); continue; end; @@ -1385,7 +1506,7 @@ begin setlength(sighupSimulationen,length(sighupSimulationen)+1); sighupSimulationen[length(sighupSimulationen)-1]:=self; - gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter'); + gitter:=tGitter.create(self,round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter'); for i:=0 to length(teilchen)-1 do teilchen[i].free; setlength(teilchen,0); diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps index de2ba33..ba80993 100644 --- a/Plasmapropagation.lps +++ b/Plasmapropagation.lps @@ -9,19 +9,20 @@ - + + - - - - + + + + - + @@ -30,16 +31,15 @@ - + - - - - + + + @@ -47,32 +47,34 @@ - - + + + - - + + + - - - - + + + + - - - + + + @@ -80,186 +82,186 @@ - + - + - + - + - + - + - + - + - + - + - - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + - - + + - - + + - - + + - - + + - - + + - + diff --git a/input.epost b/input.epost index 5a448af..39b7444 100644 --- a/input.epost +++ b/input.epost @@ -15,25 +15,47 @@ Palette 000000 Ende -Threadanzahl: 10 +Palette + Name: rotblau + ff0000 + ffffff + 0000ff +Ende -# !Schleife: $Feld: dDPHIDX dDAYDT dAY dAX -!Schleife: $Feld: AY DAYDT N1 N2 +Threadanzahl: 7 + +!Schleife: $symmetrie: symmetrisch asymmetrisch + +?$symmetrie = symmetrisch: !Schleife: $Feld: DPHIDX DPSIDX1 # AY DAYDT +?$symmetrie = asymmetrisch: !Schleife: $Feld: N1 # N2 Daten einlesen Genauigkeit: extended + xmax: 5 + tmax: 40 SpaceTime-Datei: /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld-*_test.dat -# tmin: 0 -# tmax: 5 -# xmin: 0.005 -# xmax: 5 -# Gamma: 1 Ende lineares Bild - Palette: erweiterter Regenbogen - maximale und minimale Dichten bestimmen +?$symmetrie = symmetrisch: Palette: rotblau +?$symmetrie = symmetrisch: maximale und minimale Dichten bestimmen (symmetrisch) +?$symmetrie = asymmetrisch: Palette: erweiterter Regenbogen +?$symmetrie = asymmetrisch: maximale und minimale Dichten bestimmen Datei: /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld_test.bmp + Schriftgröße: 60 + Achse: oben 10+ + Achse: links 20+ + Achse: unten 10+ + Achse: rechts 20+ + Rahmen Ende +externer Befehl: convert /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld_test.bmp /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld_test.png& + +?$symmetrie = symmetrisch: !Schleifenende +?$symmetrie = asymmetrisch: !Schleifenende !Schleifenende + +warte auf externe Befehle + +externer Befehl: rm /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/*_test.bmp diff --git a/input.plap b/input.plap index cad3cd7..94a7c54 100644 --- a/input.plap +++ b/input.plap @@ -7,10 +7,10 @@ allgemein runge-Kutta-3/8 # runge-Kutta-4 # euler-Vorwärts - ortsschritt 10^-3 * λ - zeitschritt 10^-3 * T - minZeitschritt 10^-10 * T - zeit 20 * T + ortsschritt 2*10^-3 * λ + zeitschritt 2*10^-3 * T + minZeitschritt 2*10^-3 * T + zeit 40 * T diffusionsterm 2 # max. p/Lp !setze $breite: (5 * λ) breite $breite @@ -26,12 +26,14 @@ ausgabenEnde !setze $tFwhm: (2.5 * T) !setze $tMitte: (1 * T) -licht von links 0.1 * ω * 2^(-2*((t-$tMitte)/$tFwhm)^2) * cos(ω*t) # Zeitableitung des A-Feldes +# licht von links 0.1 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * ω * cos(ω*t) # Zeitableitung des A-Feldes + +!setze $IonenMassenFaktor: (1836.15267245 + 1838.68366158) teilchen1 spezifische Ladung -q/me maximaldichte 10 -# minimaldichte 1 + minimaldichte 0 # 10^-2 !setze $profilbreite: (4 * λ) !setze $randbreite: (0.1 * λ) verteilung stückweise @@ -48,10 +50,21 @@ teilchen1 teilchen1Ende teilchen2 - spezifische Ladung (q/me)/2 - maximaldichte nmax1*2 -# minimaldichte nmin1*2 - verteilung wie teilchen1 + spezifische Ladung q/$IonenMassenFaktor/me + maximaldichte nmax1*$IonenMassenFaktor + minimaldichte nmin1*$IonenMassenFaktor +# verteilung wie teilchen1 + verteilung stückweise + 0 + (0.1*λ + $breite-$profilbreite)/2 - $randbreite + sin((x - (0.1*λ + $breite-$profilbreite)/2 - $randbreite)*π/2/$randbreite)^2 + (0.1*λ + $breite-$profilbreite)/2 + 1 + (0.1*λ + $breite+$profilbreite)/2 + sin((x - (0.1*λ + $breite+$profilbreite)/2 + $randbreite)*π/2/$randbreite)^2 + (0.1*λ + $breite+$profilbreite)/2 + $randbreite + 0 + stückweiseEnde teilchen2Ende Dateiende diff --git a/linearkombination.inc b/linearkombination.inc index 0f6ceca..6250006 100644 --- a/linearkombination.inc +++ b/linearkombination.inc @@ -55,7 +55,7 @@ begin efDAXDT,efDAYDT,efDAZDT, efDPhiDX ); *) - for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde + for emF:=erstesEMFMitAbleitung to letztesEMFMitAbleitung do emWerte[emF,false]:= in1.emWerte[emF,false] + fak2 * in2.emWerte[emF,true] {$IFDEF lkA3} @@ -97,7 +97,7 @@ begin mfGamma,mfIGamma ); *) for i:=0 to length(matWerte)-1 do // siehe oben - for maF:=mfN to mfDPsiDX do + for maF:=erstesMatFMitAbleitung to letztesMatFMitAbleitung do matWerte[i,maF,false]:= in1.matWerte[i,maF,false] + fak2 * in2.matWerte[i,maF,true] {$IFDEF lkA3} @@ -229,6 +229,8 @@ begin {$IFDEF lkA31},fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31 {$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; end; -- cgit v1.2.3-70-g09d2