From 5eaf0e1928f4024ac5c1935dfecfe7d29f93417c Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Tue, 11 Aug 2015 14:17:49 +0200 Subject: lauffaehig, es entstehen aber noch "Ripple" im Phasenraum ... --- Physikunit.pas | 228 ++++++++++++++++++++++++++++++++------------------ Plasmapropagation.lps | 141 ++++++++++++------------------- input.epost | 51 +++++++---- input.plap | 32 +++---- 4 files changed, 247 insertions(+), 205 deletions(-) diff --git a/Physikunit.pas b/Physikunit.pas index 4861a14..5d8881a 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -27,6 +27,7 @@ type eqRho,eqJX,eqJY ); tMaterieFeldGroesze = (mfPY); + tMaterieSpeicherGroesze = (msPX,msPXSqr,msPY,msVX,msVY,msN,msPXRipple); tTeilchenSpeziesGroeszen = (tsgMasse,tsgIQdrMasse,tsgLadung); tEMQuellen = array[tEMQuellGroesze] of extended; @@ -46,9 +47,11 @@ type tAusgabeDatei = class private ableitung: boolean; - emFelder: tEMFeldGroesze; - emQuellen: tEMQuellGroesze; - teilchen, // -2: EM-Feld; -1: Gesamt-EM-Quelle; sonst: EM-Quelle der entsprechenden Teilchen + emFeld: tEMFeldGroesze; + emQuelle: tEMQuellGroesze; + matFeld: tMaterieSpeicherGroesze; + wasSpeichern, // 0: EM-Feld; 1: EM-Quelle; 2: Mat-Feld + teilchen, // -1: Gesamt-Feld; sonst: Feld des entsprechenden Teilchen nNum,tCnt: longint; pre,suf: string; datei: file; @@ -132,7 +135,7 @@ type {$UNDEF LiKotRaumPunktHeader} function nichtnegativieren: extended; inline; function impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; overload; inline; - function impulsIntegral(teilchen: longint): extended; overload; inline; + function impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended; overload; inline; procedure initialisiereDichte(teilchen: longint; breite,n: extended); inline; procedure reflektiereTeilchen(rechts: boolean); inline; procedure setzeNull; inline; @@ -211,6 +214,9 @@ const emQuellNamen: array[tEMQuellGroesze] of string = ( 'RHO','JX','JY' ); + matSpeicherNamen: array[tMaterieSpeicherGroesze] of string = ( + 'PX','PXSQR','PY','VX','VY','N','PXRIPPLE' + ); minAusgabeBuffer = 1024*1024; var @@ -222,6 +228,7 @@ constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokol var emF: tEMFeldGroesze; emQ: tEMQuellGroesze; + maF: tMaterieSpeicherGroesze; num: longint; abl,gef: boolean; begin @@ -238,8 +245,10 @@ begin else teilchen:=-1; - emFelder:=efAX; - emQuellen:=eqRho; + emFeld:=efAX; + emQuelle:=eqRho; + matFeld:=msN; + wasSpeichern:=0; feldName:=ansiUpperCase(feldName); nNum:=0; // Header 0 wurde also noch nicht geschrieben zeitAnz:=-1; @@ -250,7 +259,8 @@ begin if not gef then for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin - emFelder:=emF; + emFeld:=emF; + wasSpeichern:=0; ableitung:=abl; teilchen:=-2; gef:=true; @@ -261,24 +271,45 @@ begin for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do if emQuellNamen[emQ] = feldName then begin ableitung:=false; - emQuellen:=emQ; + emQuelle:=emQ; + wasSpeichern:=1; + gef:=true; + break; + end; + if not gef then + for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do + if matSpeicherNamen[maF] = feldName then begin + ableitung:=false; + matFeld:=maF; + wasSpeichern:=2; gef:=true; break; end; if not gef then begin pro.schreibe('tAusgabeDatei.create kennt Größe '''+feldName+''' nicht!',true); + nam:=''; + for abl:=false to true do + for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do + nam:=nam+', '+copy('D',1,byte(abl))+emFeldNamen[emF]; + for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do + nam:=nam+', '+emQuellNamen[emQ]; + for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do + nam:=nam+', '+matSpeicherNamen[maF]; + delete(nam,1,2); + pro.schreibe('Ich kenne nur: '+nam,true); + pro.destroyAll; halt(1); end; - if teilchen=-2 then - pre:=emFeldNamen[emF] - else begin - pre:=emQuellNamen[emQ]; - if teilchen>=0 then - pre:=pre+inttostr(teilchen+1); + case wasSpeichern of + 0: pre:=emFeldNamen[emFeld]; + 1: pre:=emQuellNamen[emQuelle]; + 2: pre:=matSpeicherNamen[matFeld]; end; + if teilchen>=0 then + pre:=pre+inttostr(teilchen+1); if ableitung then pre:='d'+pre; @@ -320,13 +351,13 @@ end; function tAusgabeDatei.dump: string; begin - if teilchen=-2 then - result:=emFeldNamen[emFelder] - else begin - result:=emQuellNamen[emQuellen]; - if teilchen>=0 then - result:=result+inttostr(teilchen); + case wasSpeichern of + 0: result:=emFeldNamen[emFeld]; + 1: result:=emQuellNamen[emQuelle]; + 2: result:=matSpeicherNamen[matFeld]; end; + if teilchen>=0 then + result:=result+'['+inttostr(teilchen+1)+']'; if ableitung then result:=result+''''; result:=result+' -> '+pre+'-$i'+suf; @@ -385,40 +416,40 @@ begin inc(bufPos,sizeof(integer)); sX:=cX-Min(gitter.dX,sDX)/2; - if teilchen=-2 then begin - for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin - while cX>=sX do begin - dec(cnt); - move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFelder,ableitung],(buf+bufPos)^,sizeof(extended)); - inc(bufPos,sizeof(extended)); - sX:=sX+sDX; + case wasSpeichern of + 0: // em-Feld speichern + for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin + while cX>=sX do begin + dec(cnt); + move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFeld,ableitung],(buf+bufPos)^,sizeof(extended)); + inc(bufPos,sizeof(extended)); + sX:=sX+sDX; + end; + cX:=cX+gitter.dX; end; - cX:=cX+gitter.dX; - end; - end - else if teilchen=-1 then begin - for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin - while cX>=sX do begin - dec(cnt); - move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emQuellen[emQuellen],(buf+bufPos)^,sizeof(extended)); - inc(bufPos,sizeof(extended)); - sX:=sX+sDX; + 1: // em-Quelle + for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin + while cX>=sX do begin + dec(cnt); + val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuelle); + move(val,(buf+bufPos)^,sizeof(extended)); + inc(bufPos,sizeof(extended)); + sX:=sX+sDX; + end; + cX:=cX+gitter.dX; end; - cX:=cX+gitter.dX; - end; - end - else begin - for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin - while cX>=sX do begin - dec(cnt); - val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuellen); - move(val,(buf+bufPos)^,sizeof(extended)); - inc(bufPos,sizeof(extended)); - sX:=sX+sDX; + 2: // Materiefeld + for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin + while cX>=sX do begin + dec(cnt); + val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,matFeld); + move(val,(buf+bufPos)^,sizeof(extended)); + inc(bufPos,sizeof(extended)); + sX:=sX+sDX; + end; + cX:=cX+gitter.dX; end; - cX:=cX+gitter.dX; - end; - end; + end{of case}; if cnt<>0 then begin pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true); pro.destroyall; @@ -590,17 +621,18 @@ end; function tImpulsPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen var - i,richtung,entfernung: longint; defizit: extended; + i: longint; pro: tProtokollant; - jemandErreicht,vorZurueck: boolean; +(* i,richtung,entfernung: longint; + jemandErreicht,vorZurueck: boolean; *) begin result:=0; - for i:=0 to length(werte)-1 do + for i:=0 to length(werte)-1 do begin if werte[i,false]<0 then begin defizit:=-werte[i,false]; result:=result+defizit; - werte[i,false]:=0; + werte[i,false]:=0; (* entfernung:=1; jemandErreicht:=true; while (defizit>0) and jemandErreicht do begin @@ -616,8 +648,15 @@ begin pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'WertePunkt.nichtnegativieren'); pro.schreibe('Ich konnte ein Teilchendefizit der Sorte '+inttostr(i+1)+' nicht auffüllen!',true); raumpunkt.felder.gitter.abbrechen; - end; + end; *) end; + if werte[i,false]>1E10 then begin + pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'impulsPunkt.nichtnegativieren'); + pro.schreibe('An einer Stelle im Phasenraum sind bereits mehr als 10^100 Teilchen, ich breche ab!',true); + pro.free; + raumpunkt.felder.gitter.abbrechen; + end; + end; end; procedure tImpulsPunkt.akkumuliereEMQuellen(var emQuellen: tEMQuellen; pX,dVP: extended); @@ -688,16 +727,16 @@ begin + raumpunkt.matFelder[i,mfPY,false] * iMGamma[i] * raumpunkt.emFelder[efBZ,false]); werte[i,true]:= 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]) + + 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]) + + sign(tP)*raumpunkt.felder.teilchen[i].diffusion[1] * ( nachbarn[1,true].werte[i,false] + nachbarn[1,false].werte[i,false] - 2*werte[i,false] @@ -846,8 +885,6 @@ begin end; function tRaumPunkt.impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; -var - i: longint; begin if teilchen<0 then begin result:=emQuellen[emQ]; // das ist leicht :-) @@ -857,12 +894,50 @@ begin result:=0; case emQ of - eqRho: + eqRho: result:=impulsIntegral(teilchen,msN); + eqJX: result:=impulsIntegral(teilchen,msVX); + eqJY: result:=impulsIntegral(teilchen,msVY); + end{of case}; + result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung]; +end; + +function tRaumPunkt.impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended; +var + i: longint; +begin + if teilchen<0 then begin + result:=0; + for i:=0 to length(felder.teilchen)-1 do + result:=result+impulsIntegral(i,maF); + exit; + end; + + result:=0; + + case maF of + msN: for i:=0 to aP-1 do result:= result + phasenraum[i].werte[teilchen,false]*dP; - eqJX: begin + msPX: + for i:=0 to aP-1 do + result:= + result + + phasenraum[i].werte[teilchen,false]*dP * (i-aP/2)*dP; + msPXSqr: + for i:=0 to aP-1 do + result:= + result + + phasenraum[i].werte[teilchen,false]*dP * sqr((i-aP/2)*dP); + msPY: begin + for i:=0 to aP-1 do + result:= + result + + phasenraum[i].werte[teilchen,false]*dP; + result:=result * matFelder[teilchen,mfPY,false]; + end; + msVX: begin for i:=0 to aP-1 do result:= result + @@ -870,7 +945,7 @@ begin / 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 + msVY: begin for i:=0 to aP-1 do result:= result + @@ -878,24 +953,11 @@ begin / 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; + msPXRipple: + for i:=1 to aP-1 do + result:= + result + abs(phasenraum[i].werte[teilchen,false]-phasenraum[i-1].werte[teilchen,false])*dP; end{of case}; - result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung]; -end; - -function tRaumPunkt.impulsIntegral(teilchen: longint): extended; -var - 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); @@ -1086,7 +1148,7 @@ begin 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+inhalt[j].impulsIntegral(i,msN); dens:=dens*teilchen[i].eigenschaften[tsgMasse]*gitter.dX; prot.schreibe('n['+inttostr(i+1)+'] = '+floattostr(dens)+' (relative Abweichung: '+floattostr(dens/massen[i]-1)+')',true); end; diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps index d47f765..9a7b886 100644 --- a/Plasmapropagation.lps +++ b/Plasmapropagation.lps @@ -3,25 +3,25 @@ - + - + + - - - - - + + + + - + @@ -29,15 +29,17 @@ - - + + + - - + + + @@ -45,24 +47,23 @@ - - - + + + - - - + + - - - + + + @@ -77,9 +78,9 @@ - + - + @@ -144,131 +145,93 @@ - - + + + + + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/input.epost b/input.epost index 96ec796..6558880 100644 --- a/input.epost +++ b/input.epost @@ -26,50 +26,65 @@ Threadanzahl: 11 parallel lesen -!setze $symFelder: JX1 AY DAYDT # EX -!setze $asyFelder: RHO1 +!setze $symFelder: VX1 # AY DAYDT # EX +!setze $asyFelder: N1 PXRIPPLE1 +!setze $reineInputFelder: PX1 PXSQR1 -!Schleife: $Feld: $asyFelder $symFelder +!Schleife: $Feld: $reineInputFelder $asyFelder $symFelder Daten einlesen Name: $Feld Genauigkeit: extended xmax: 5 - tmax: 40 + tmax: 2 SpaceTime-Datei: /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld-*_test.dat Ende !Schleifenende -!setze $symBerFelder: VX1 +!setze $asyBerFelder: PXBreite -Teile JX1 durch RHO1 - Name: VX1 +Teile VX1 durch N1 zu VX1 Ende -!Schleife: $symmetrie: symmetrisch asymmetrisch +Teile PX1 durch N1 zu PX1 +Ende + +Teile PXRIPPLE1 durch N1 zu PXRIPPLE1 +Ende + +Teile PXSQR1 durch N1 zu PXSQR1 +Ende + +Multipliziere PX1 mal PX1 + Name: PX1SQR +Ende + +Linearkombination + Name: PXBreite + PXSQR1 1 + PX1SQR -1 +Ende -?$symmetrie = symmetrisch: !Schleife: $Feld: $symFelder $symBerFelder -?$symmetrie = asymmetrisch: !Schleife: $Feld: $asyFelder +!Schleife: $Feld: $symFelder $asyFelder $asyBerFelder lineares Bild $Feld -?$symmetrie = symmetrisch: Palette: rotblau -?$symmetrie = symmetrisch: maximale und minimale Dichten bestimmen (symmetrisch) -?$symmetrie = asymmetrisch: Palette: erweiterter Regenbogen -?$symmetrie = asymmetrisch: maximale und minimale Dichten bestimmen +?$Feld in $symFelder $symBerFelder: Palette: rotblau +?$Feld in $symFelder $symBerFelder: maximale und minimale Dichten bestimmen (symmetrisch) +?$Feld in $asyFelder $asyBerFelder: Palette: erweiterter Regenbogen +?$Feld in $asyFelder $asyBerFelder: maximale und minimale Dichten bestimmen +?$Feld = N1: Nachbearbeitung: Log: 10^-8 Datei: /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld_test.bmp Schriftgröße: 60 Achse: oben 10+ - Achse: links 20+ + Achse: links 2+ Achse: unten 10+ - Achse: rechts 20+ + Achse: rechts 2+ 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 diff --git a/input.plap b/input.plap index a6905dd..f83acba 100644 --- a/input.plap +++ b/input.plap @@ -7,11 +7,11 @@ allgemein runge-Kutta-3/8 # runge-Kutta-4 # euler-Vorwärts - ortsschritt 5*10^-3 * λ - zeitschritt 5*10^-3 * T - maximalimpuls 10 - impulsschritt 5 * 10^-3 - zeit 40 * T + ortsschritt 10^-2 * λ + zeitschritt 10^-2 * T + maximalimpuls 5 + impulsschritt 10^-2 + zeit 2 * T !setze $breite: (5 * λ) breite $breite mit Fortschrittsanzeige @@ -20,13 +20,13 @@ allgemeinEnde ausgaben prefix /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/ suffix _test.dat - felder AX,AY,dAYDT,EX,BY,BZ,Rho1,Rho2,JX1 + felder AX,AY,dAYDT,EX,BY,BZ,Rho1,Rho2,VX1,PX1,PXSqr1,N1,PXRipple1 ausgabenEnde !setze $tFwhm: (2.5 * T) !setze $tMitte: (1 * T) -licht von links 2 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * ω * cos(ω*t) # Zeitableitung des A-Feldes +# licht von links 2 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * ω * cos(ω*t) # Zeitableitung des A-Feldes !setze $IonenMassenFaktor: (1836.15267245 + 1838.68366158) @@ -36,7 +36,8 @@ teilchen1 maximaldichte 10 minimaldichte 0 # 10^-2 maximales dLnN/dX 10 - maximales dLnN/dP 10 + maximales dLnN/dP 2 + impulsbreite 1/20 !setze $profilbreite: (4 * λ) !setze $randbreite: (0.1 * λ) verteilung stückweise @@ -58,17 +59,18 @@ teilchen2 maximaldichte nmax1 minimaldichte nmin1 maximales dLnN/dX 10 - maximales dLnN/dP 10 + maximales dLnN/dP 2 + impulsbreite 1/20 # verteilung wie teilchen1 verteilung stückweise 0 - ($breite-$profilbreite)/2 - $randbreite - sin((x - ($breite-$profilbreite)/2 - $randbreite)*π/2/$randbreite)^2 - ($breite-$profilbreite)/2 + (0.1 + $breite-$profilbreite)/2 - $randbreite + sin((x - (0.1 + $breite-$profilbreite)/2 - $randbreite)*π/2/$randbreite)^2 + (0.1 + $breite-$profilbreite)/2 1 - ($breite+$profilbreite)/2 - sin((x - ($breite+$profilbreite)/2 + $randbreite)*π/2/$randbreite)^2 - ($breite+$profilbreite)/2 + $randbreite + (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 -- cgit v1.2.3-70-g09d2