diff options
| -rw-r--r-- | Physikunit.pas | 233 | ||||
| -rw-r--r-- | Plasmapropagation.lps | 134 | ||||
| -rw-r--r-- | input.epost | 42 | ||||
| -rw-r--r-- | input.plap | 33 | ||||
| -rw-r--r-- | 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) + <rho> * 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) + <dRho/dt> * 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:=dTMax<dT*2; + result:=(dTMax<dT) and (dT>1.1*dTMin); if result then begin // beabsichtigter Zeitschritt ist zu groß, dT:=dTMax/4; // dann machen wir ihn doch kleiner! + if dT<dTMin then + dT:=dTMin {$IFDEF Zeitschrittueberwachung} - prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT)+' (t = '+floattostr(t)+')'); - if dT<dTMin then begin - prot.schreibe('pruefeMaxDT: Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < '+floattostr(dTMin)+') - irgendwas scheint grundsätzlich kaputt zu sein!',true); - abbrechen; - end; - prot.spuelen; - {$ENDIF} - exit; + else begin + prot.schreibe('pruefeMaxDT: 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 @@ <IsPartOfProject Value="True"/> <TopLine Value="54"/> <CursorPos X="37" Y="34"/> - <UsageCount Value="188"/> + <UsageCount Value="199"/> <Loaded Value="True"/> </Unit0> <Unit1> <Filename Value="Physikunit.pas"/> <IsPartOfProject Value="True"/> + <IsVisibleTab Value="True"/> <EditorIndex Value="1"/> - <TopLine Value="943"/> - <CursorPos Y="1413"/> - <FoldState Value=" T3lf0,5 pjSlM0A[94BjZ084]9Ija09]9ifkK07[I45jO0I013 ppasC0G2 T0=HQ/0F17121b"/> - <UsageCount Value="129"/> + <TopLine Value="619"/> + <CursorPos X="13" Y="648"/> + <FoldState Value=" T3mD0,5 pjSlM0A[94BjZ0812114]DLnd0B[944jR08013]Bdl20U2 pp6q10G2 T0\9Q/0F171A"/> + <UsageCount Value="140"/> <Bookmarks Count="1"> - <Item0 Y="974"/> + <Item0 Y="1078"/> </Bookmarks> <Loaded Value="True"/> </Unit1> @@ -30,16 +31,15 @@ <IsPartOfProject Value="True"/> <EditorIndex Value="-1"/> <CursorPos X="15" Y="46"/> - <UsageCount Value="91"/> + <UsageCount Value="102"/> </Unit2> <Unit3> <Filename Value="input.plap"/> <IsPartOfProject Value="True"/> - <IsVisibleTab Value="True"/> <EditorIndex Value="3"/> - <TopLine Value="6"/> - <CursorPos X="19" Y="14"/> - <UsageCount Value="90"/> + <TopLine Value="31"/> + <CursorPos X="26" Y="36"/> + <UsageCount Value="101"/> <Loaded Value="True"/> <DefaultSyntaxHighlighter Value="None"/> </Unit3> @@ -47,32 +47,34 @@ <Filename Value="linearkombination.inc"/> <IsPartOfProject Value="True"/> <EditorIndex Value="2"/> - <CursorPos X="48" Y="220"/> - <UsageCount Value="38"/> + <TopLine Value="88"/> + <CursorPos X="50" Y="100"/> + <UsageCount Value="49"/> <Loaded Value="True"/> </Unit4> <Unit5> <Filename Value="input.epost"/> <EditorIndex Value="4"/> - <CursorPos X="37" Y="22"/> - <UsageCount Value="78"/> + <TopLine Value="13"/> + <CursorPos X="5" Y="22"/> + <UsageCount Value="84"/> <Loaded Value="True"/> <DefaultSyntaxHighlighter Value="None"/> </Unit5> <Unit6> <Filename Value="../units/matheunit.pas"/> <EditorIndex Value="-1"/> - <TopLine Value="53"/> - <CursorPos Y="53"/> - <FoldState Value=" T3i905B pj0jV034 piaj60U2-"/> - <UsageCount Value="17"/> + <TopLine Value="366"/> + <CursorPos X="27" Y="390"/> + <FoldState Value=" T3i905B pj0jV034 piaj60U511+"/> + <UsageCount Value="16"/> </Unit6> <Unit7> <Filename Value="../units/lowlevelunit.pas"/> <EditorIndex Value="-1"/> - <TopLine Value="4"/> - <CursorPos X="86" Y="23"/> - <UsageCount Value="7"/> + <TopLine Value="46"/> + <CursorPos X="3" Y="12"/> + <UsageCount Value="9"/> </Unit7> <Unit8> <Filename Value="../units/mystringlistunit.pas"/> @@ -80,186 +82,186 @@ <TopLine Value="367"/> <CursorPos X="17" Y="390"/> <FoldState Value=" T3i3075 piZjD0WQ"/> - <UsageCount Value="8"/> + <UsageCount Value="7"/> </Unit8> <Unit9> <Filename Value="../epost/werteunit.pas"/> <EditorIndex Value="-1"/> <TopLine Value="950"/> <CursorPos X="30" Y="1054"/> - <UsageCount Value="5"/> + <UsageCount Value="4"/> </Unit9> <Unit10> <Filename Value="../epost/typenunit.pas"/> <EditorIndex Value="-1"/> <TopLine Value="347"/> <CursorPos X="62" Y="358"/> - <UsageCount Value="5"/> + <UsageCount Value="4"/> </Unit10> <Unit11> <Filename Value="../units/systemunit.pas"/> <EditorIndex Value="-1"/> <CursorPos X="3" Y="79"/> - <UsageCount Value="5"/> + <UsageCount Value="4"/> </Unit11> <Unit12> <Filename Value="/usr/lib/fpc/src/rtl/inc/objpash.inc"/> <EditorIndex Value="-1"/> <TopLine Value="232"/> <CursorPos X="23" Y="192"/> - <UsageCount Value="5"/> + <UsageCount Value="4"/> </Unit12> <Unit13> <Filename Value="rk14.inc"/> <EditorIndex Value="-1"/> - <UsageCount Value="8"/> + <UsageCount Value="7"/> </Unit13> <Unit14> <Filename Value="rk1210.inc"/> <EditorIndex Value="-1"/> <TopLine Value="492"/> <CursorPos X="28" Y="565"/> - <UsageCount Value="9"/> + <UsageCount Value="8"/> </Unit14> <Unit15> <Filename Value="rk1412.inc"/> <EditorIndex Value="-1"/> <TopLine Value="945"/> <CursorPos X="5" Y="985"/> - <UsageCount Value="9"/> + <UsageCount Value="8"/> </Unit15> <Unit16> <Filename Value="rk108.inc"/> <EditorIndex Value="-1"/> <CursorPos X="28" Y="8"/> - <UsageCount Value="9"/> + <UsageCount Value="8"/> </Unit16> <Unit17> <Filename Value="rk3_8.inc"/> <EditorIndex Value="-1"/> <CursorPos X="75" Y="23"/> - <UsageCount Value="10"/> + <UsageCount Value="9"/> </Unit17> </Units> <JumpHistory Count="30" HistoryIndex="29"> <Position1> <Filename Value="Physikunit.pas"/> - <Caret Line="341" Column="8" TopLine="320"/> </Position1> <Position2> <Filename Value="Physikunit.pas"/> - <Caret Line="1430" Column="42" TopLine="1419"/> + <Caret Line="126" Column="17" TopLine="93"/> </Position2> <Position3> <Filename Value="Physikunit.pas"/> - <Caret Line="648" Column="104" TopLine="619"/> + <Caret Line="521" Column="50" TopLine="433"/> </Position3> <Position4> <Filename Value="Physikunit.pas"/> - <Caret Line="608" Column="52" TopLine="213"/> + <Caret Line="550" Column="51" TopLine="509"/> </Position4> <Position5> <Filename Value="Physikunit.pas"/> - <Caret Line="977" Column="63" TopLine="953"/> + <Caret Line="551" Column="51" TopLine="510"/> </Position5> <Position6> <Filename Value="Physikunit.pas"/> - <Caret Line="92" Column="38" TopLine="72"/> + <Caret Line="665" Column="28" TopLine="326"/> </Position6> <Position7> <Filename Value="Physikunit.pas"/> - <Caret Line="947" Column="12" TopLine="910"/> + <Caret Line="798" Column="44" TopLine="765"/> </Position7> <Position8> <Filename Value="Physikunit.pas"/> - <Caret Line="930" Column="64" TopLine="707"/> + <Caret Line="799" Column="25" TopLine="766"/> </Position8> <Position9> <Filename Value="Physikunit.pas"/> - <Caret Line="166" Column="72" TopLine="145"/> + <Caret Line="800" Column="34" TopLine="767"/> </Position9> <Position10> <Filename Value="Physikunit.pas"/> - <Caret Line="175" Column="60" TopLine="165"/> + <Caret Line="801" Column="51" TopLine="774"/> </Position10> <Position11> <Filename Value="Physikunit.pas"/> - <Caret Line="946" Column="84" TopLine="929"/> + <Caret Line="67" Column="22" TopLine="47"/> </Position11> <Position12> <Filename Value="Physikunit.pas"/> - <Caret Line="945" Column="20" TopLine="909"/> + <Caret Line="425" Column="4" TopLine="210"/> </Position12> <Position13> <Filename Value="Physikunit.pas"/> - <Caret Line="1423" Column="40" TopLine="1394"/> + <Caret Line="438" Column="63" TopLine="216"/> </Position13> <Position14> <Filename Value="Physikunit.pas"/> - <Caret Line="92" Column="38" TopLine="72"/> + <Caret Line="801" Column="51" TopLine="768"/> </Position14> <Position15> <Filename Value="Physikunit.pas"/> - <Caret Line="605" Column="46" TopLine="585"/> + <Caret Line="1385" Column="60" TopLine="1232"/> </Position15> <Position16> <Filename Value="Physikunit.pas"/> - <Caret Line="775" Column="39" TopLine="704"/> + <Caret Line="789" Column="63" TopLine="773"/> </Position16> <Position17> <Filename Value="Physikunit.pas"/> - <Caret Line="797" Column="38" TopLine="778"/> + <Caret Line="10" Column="28"/> </Position17> <Position18> <Filename Value="Physikunit.pas"/> - <Caret Line="131" Column="35" TopLine="111"/> + <Caret Line="1024" Column="31" TopLine="784"/> </Position18> <Position19> <Filename Value="Physikunit.pas"/> - <Caret Line="775" Column="39" TopLine="753"/> + <Caret Line="1067" Column="29" TopLine="1056"/> </Position19> <Position20> <Filename Value="Physikunit.pas"/> - <Caret Line="946" Column="38" TopLine="910"/> + <Caret Line="1019" TopLine="861"/> </Position20> <Position21> <Filename Value="Physikunit.pas"/> - <Caret Line="166" Column="25" TopLine="158"/> + <Caret Line="11" Column="3"/> </Position21> <Position22> <Filename Value="Physikunit.pas"/> - <Caret Line="901" Column="17" TopLine="708"/> + <Caret Line="1023" Column="5" TopLine="1000"/> </Position22> <Position23> <Filename Value="Physikunit.pas"/> - <Caret Line="889" Column="18" TopLine="860"/> + <Caret Line="9" Column="33"/> </Position23> <Position24> <Filename Value="Physikunit.pas"/> + <Caret Line="991" Column="36" TopLine="787"/> </Position24> <Position25> - <Filename Value="Physikunit.pas"/> - <Caret Line="1423" Column="37" TopLine="1394"/> + <Filename Value="input.epost"/> + <Caret Line="27" TopLine="5"/> </Position25> <Position26> - <Filename Value="Physikunit.pas"/> - <Caret TopLine="145"/> + <Filename Value="input.epost"/> + <Caret Line="44" Column="15" TopLine="5"/> </Position26> <Position27> - <Filename Value="Physikunit.pas"/> - <Caret Line="1412" TopLine="774"/> + <Filename Value="input.epost"/> + <Caret Line="56" Column="27" TopLine="19"/> </Position27> <Position28> - <Filename Value="Physikunit.pas"/> - <Caret Line="175" Column="54" TopLine="161"/> + <Filename Value="input.epost"/> + <Caret Line="47" TopLine="25"/> </Position28> <Position29> - <Filename Value="Physikunit.pas"/> - <Caret Line="946" Column="12" TopLine="774"/> + <Filename Value="input.epost"/> + <Caret Line="57" TopLine="25"/> </Position29> <Position30> <Filename Value="Physikunit.pas"/> - <Caret Line="1424" Column="36" TopLine="988"/> + <Caret Line="989" Column="35" TopLine="903"/> </Position30> </JumpHistory> </ProjectSession> 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 @@ -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; |