From 9e540f01e511e9928460c165dc4a420dd24ba7cb Mon Sep 17 00:00:00 2001 From: Erich Eckner Date: Fri, 24 Jul 2015 15:49:04 +0200 Subject: Tagesendstand, Umstrukturierung angefangen, bei weitem nicht lauffaehig --- Physikunit.pas | 327 +++++++++++++++++++++++++++++--------------------- Plasmapropagation | Bin 1305003 -> 3344888 bytes Plasmapropagation.lpi | 20 +-- Plasmapropagation.lpr | 166 +++---------------------- Plasmapropagation.lps | 183 +++++++++++++++++----------- input.plap | 49 ++++++++ 6 files changed, 380 insertions(+), 365 deletions(-) create mode 100644 input.plap diff --git a/Physikunit.pas b/Physikunit.pas index 0aee452..71e6a79 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -8,7 +8,7 @@ unit Physikunit; interface uses - Classes, SysUtils, Math; + Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit; type TZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier); @@ -21,73 +21,70 @@ type fiGamma,fiiGamma); TAusgabeDatei = record - Ableitung: Boolean; - Inhalt: TFeldInhalt; - Datei: File; + ableitung: boolean; + inhalt: tFeldInhalt; + name: string; + datei: file; end; - TProtokolldatei = record - datei: textfile; - einrueckung: longint; - end; + tFelder = class; - { TProtokollant } + { tWerteGitter } - TProtokollant = class(TObject) - private - sDat: ^TProtokolldatei; - Bes: String; - Kinder: array of TProtokollant; - Elter: TProtokollant; - public - constructor create(Dateiname: string); overload; - constructor create(Elter_: TProtokollant; Besitzer: string); overload; - destructor destroy; override; - procedure schreibe(s: string); - procedure schreibe(s: string; tee: boolean); - procedure destroyall; + tWerteGitter = class(tObject) // repräsentiert ein Gitter von Werten und deren Zeitableitungen + werte: array[boolean] of array of extended; + par: tFelder; end; - { TFeld } + { TFelder } - pTFeld = ^TFeld; - TFeld = class(TObject) + tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen public - Groessen: array[Boolean,TFeldInhalt] of extended; + materieFelder: array of tWerteGitter; // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen + elektroMagnetischeFelder: tWerteGitter; // EM-Felder und deren Zeitableitungen // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier! dX,iDX,x,pDNMax: extended; - lN,rN: pTFeld; + constructor create(deltaX,pDNMa: extended); procedure berechneGammaUndP; procedure berechneNAbleitung; procedure berechneAbleitungen(dT: extended); end; - { TGitter } + { tGitter } - TGitter = class(TObject) + tGitter = class(tObject) private groesse: longint; - Prot: TProtokollant; + prot: tProtokollant; procedure berechneAbleitungen(Felder: longint; dt: extended); overload; procedure berechneAbleitungen(Felder: longint; dt: extended; out dTMax: extended); overload; procedure setzeRaender(Felder: longint); public aktuelleFelder: longint; - Felders: array of array of TFeld; + felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern dTMaximum,dX,iDX,xl,xr,t,pDNMax: extended; pDNMaxDynamisch: boolean; Zeitverfahren: TZeitverfahren; Al,Ar: TVerteilungsfunktion; - constructor create(size: longint; deltaT,deltaX,pDNMa: extended; n0: TVerteilungsfunktion; ZV: TZeitverfahren; Protokollant: TProtokollant); + constructor create(size: longint; deltaT,deltaX,pDNMa: extended; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant); destructor destroy; override; procedure iteriereSchritt(var dT: extended); procedure macheAusgabe(AD: TAusgabedatei; sDX: extended); function gibErhaltungsgroessen: string; end; -function Nullfunktion(x: extended): extended; -function timetostr(t: extended): string; + { tSimulation } + + tSimulation = class(tObject) + private + prot: tProtokollant; + gitter: tGitter; + kvs: tKnownValues; + public + constructor create(inName: string; Protokollant: tProtokollant); + destructor destroy; override; + end; const OptionNamen: array[TFeldInhalt] of string = ( @@ -98,81 +95,6 @@ const implementation -{ TProtokollant } - -constructor TProtokollant.create(Dateiname: string); -begin - inherited create; - getmem(sDat,sizeof(TProtokolldatei)); - sDat^.einrueckung:=10; - assignfile(sDat^.datei,Dateiname); - rewrite(sDat^.datei); - Bes:=''; - setlength(Kinder,0); - Elter:=nil; -end; - -constructor TProtokollant.create(Elter_: TProtokollant; Besitzer: string); -begin - inherited create; - sDat:=Elter_.sDat; - setlength(Elter_.Kinder,length(Elter_.Kinder)+1); - Elter_.Kinder[length(Elter_.Kinder)-1]:=self; - Elter:=Elter_; - setlength(Kinder,0); - Bes:=Elter.Bes+'.'+Besitzer; - if pos('.',Bes)=1 then delete(Bes,1,1); - if length(Bes)+4 > sDat^.einrueckung then - sDat^.einrueckung:=length(Bes)+4; -end; - -destructor TProtokollant.destroy; -var i: longint; -begin - while length(Kinder)>0 do - Kinder[length(Kinder)-1].free; - if assigned(Elter) then begin - flush(sDat^.datei); - i:=0; - while (iself) do - inc(i); - if i>=length(Elter.Kinder) then - schreibe('destroy fehlgeschlagen, ich kann mich nicht sehen.'); - inc(i); - while i=1; - if b then begin - result:=result+inttostr(floor(t))+' Tage '; - t:=t-floor(t); + ifile:=tMyStringlist.create(prot); + ifile.loadfromfile(inName); + if not ifile.unfoldMacros then begin + prot.schreibe('Fehlerhafte Macros in Parameterdatei '''+inName+'''!',true); + halt(1); + end; + + // Standardeinstellungen Bereich 'allgemein' + zeitverfahren:=zfRungeKuttaVier; + deltaX:=1e-2; + kvs.add('λ',1/deltaX); + deltaT:=-1; + sDT:=-1; + endzeit:=100; + breite:=10.0; + pDNMax:=-1; + + // Standardeinstellungen Bereich 'ausgaben' + aPrefix:=extractfilepath(paramstr(0)); + aSuffix:='.dat'; + setlength(ausgabeDateien,0); + + while ifile.readln(s) do begin + if s='allgemein' then begin + repeat + if not ifile.readln(s) then begin + prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich allgemein!',true); + halt(1); + end; + if s='allgemeinEnde' then break; + if s='runge-Kutta-4' then begin + Zeitverfahren:=zfRungeKuttaVier; + continue; + end; + if s='euler-Vorwärts' then begin + Zeitverfahren:=zfEulerVorwaerts; + continue; + end; + if startetMit('ortsschritt',s) then begin + deltaX:=exprtofloat(false,s,kvs,nil); + kvs.add('λ',1/deltaX); + continue; + end; + if startetMit('zeitschritt',s) then begin + deltaT:=exprtofloat(false,s,kvs,nil); + kvs.add('dT',deltaT); + continue; + end; + if startetMit('diffusionsterm',s) then begin + pDNMax:=exprtofloat(false,s,kvs,nil); + continue; + end; + if startetMit('zeit',s) then begin + endzeit:=exprtofloat(false,s,kvs,nil); + continue; + end; + if startetMit('breite',s) then begin + breite:=exprtofloat(false,s,kvs,nil); + continue; + end; + prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich allgemein!',true); + halt(1); + until false; + continue; + end; + + if s='ausgaben' then begin + repeat + if not ifile.readln(s) then begin + prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich ausgaben!',true); + halt(1); + end; + if s='ausgabenEnde' then break; + if startetMit('suffix',s) then begin + aSuffix:=s; + continue; + end; + if startetMit('prefix',s) then begin + aPrefix:=s; + continue; + end; + if startetMit('zeitschritt',s) then begin + sDT:=exprtofloat(false,s,kvs,nil); + continue; + end; + if startetMit('felder',s) then begin + s:=s+','; + while s<>'' do begin + t:=erstesArgument(s,','); + for abl:=false to true do + for fI:=low(tFeldInhalt) to high(tFeldInhalt) do + if not (abl and (fI in [fiGamma,fiP])) then + if t=copy('d',1,byte(abl))+optionNamen[fI] then begin + setlength(ausgabedateien,length(ausgabedateien)+1); + ausgabedateien[length(ausgabedateien)-1].inhalt:=fI; + ausgabedateien[length(ausgabedateien)-1].ableitung:=abl; + ausgabedateien[length(ausgabedateien)-1].name:=aPrefix+t+aSuffix; + assignFile(ausgabedateien[length(ausgabedateien)-1].datei,ausgabedateien[length(ausgabedateien)-1].name); + end; + end; + continue; + end; + prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true); + halt(1); + until false; + continue; + end; + + prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true); + halt(1); end; - t:=t*24; - b:=b or (t>=1); - if b then begin - result:=result+inttostr(floor(t))+'h '; - t:=t-floor(t); + + if length(ausgabedateien)=0 then begin + prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true); + halt(1); end; - t:=t*60; - b:=b or (t>=1); - if b then begin - result:=result+inttostr(floor(t))+'min '; - t:=t-floor(t); + + if deltaT<0 then + deltaT:=deltaX/10; + if sDT<0 then + sDT:=deltaT; + + for i:=0 to length(ausgabedateien)-1 do begin + rewrite(ausgabedateien[i].datei,1); + j:=0; + BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); + j:=floor(Endzeit/sDT+1); + BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); end; - t:=t*60; - result:=result+inttostr(round(t))+'s'; + ifile.free; + + gitter:=tGitter.create(Breite,deltaT,deltaX,pDNMax,@Anfangsdichte,Zeitverfahren,Prot); +end; + +destructor tSimulation.destroy; +begin + gitter.free; + kvs.free; + prot.free; + inherited destroy; end; end. diff --git a/Plasmapropagation b/Plasmapropagation index c461c9f..013329a 100755 Binary files a/Plasmapropagation and b/Plasmapropagation differ diff --git a/Plasmapropagation.lpi b/Plasmapropagation.lpi index 40504f6..f707ff5 100644 --- a/Plasmapropagation.lpi +++ b/Plasmapropagation.lpi @@ -1,4 +1,4 @@ - + @@ -32,17 +32,23 @@ - + - - + + + + + + + + @@ -52,6 +58,7 @@ + @@ -73,13 +80,8 @@ - - - - - diff --git a/Plasmapropagation.lpr b/Plasmapropagation.lpr index 52e1926..b861d2e 100644 --- a/Plasmapropagation.lpr +++ b/Plasmapropagation.lpr @@ -8,7 +8,7 @@ uses {$ENDIF}{$ENDIF} Classes, SysUtils, CustApp, { you can add units after this } - math, Physikunit, crt; + math, Physikunit, crt, protokollunit; var ErrorCode: longint; @@ -22,145 +22,34 @@ type public constructor Create(TheOwner: TComponent); override; destructor Destroy; override; - procedure WriteHelp; virtual; end; -function InputFeld(t: extended): extended; -const a0 = 1; - omega = 2*pi*100; - tfwhm = 10 * 2*pi/omega; - tMitte = 3*tfwhm; -begin -(* result:= - a0 * Exp(-2*Ln(2)*sqr((t-tMitte)/tfwhm)) * - (omega*Cos(t*omega) + 4*Ln(2)*(tMitte - t)*Sin(t*omega)/sqr(tfwhm)); *) - - result:= - a0 * Exp(-2*Ln(2)*sqr((t-tMitte)/tfwhm)) * - Sin(t*omega) -end; - -function Anfangsdichte(x: extended): extended; -const Profilbreite = 0.5; - Randbreite = 0.1; - offset = 0.0001; - Mitte = 0.75; -begin - x:=abs(x-Mitte); - result:=1+offset; - if x<=Profilbreite then exit; - result:=0+offset; - if x>=Profilbreite+Randbreite then exit; - result:=sqr(sin((Profilbreite+Randbreite-x)*0.5*pi/Randbreite))+offset; -end; - { TPlasmapropagation } procedure TPlasmapropagation.DoRun; var - ErrorMsg,s,t,u: String; - Breite,i,j: longint; - Gitter: TGitter; - deltaT,deltaX,start, - sT,sDT,sDX,Endzeit, - zeitDatei,zeitPhysik, - pDNMax: extended; - AusgabeDateien: array of TAusgabeDatei; - FI: TFeldInhalt; - Zeitverfahren: TZeitverfahren; - Abl: Boolean; - Prot: TProtokollant; - c: char; + Breite,i,j: longint; + simulation: tSimulation; + start,zeitDatei,zeitPhysik: extended; + FI: TFeldInhalt; + Abl: Boolean; + Prot: TProtokollant; + c: char; begin - // quick check parameters - start:=now; - s:=''; - for Abl:=False to True do - for FI:=low(TFeldInhalt) to high(TFeldInhalt) do - if not (Abl and (FI in [fiP,fiPX,fiPY,fiPZ,fiGamma,fiiGamma])) then - s:=s+' '+Copy('d',1,Byte(Abl))+OptionNamen[FI]+':'; - ErrorMsg:=CheckOptions('HB:Z:x:t:FEX:T:D:','Hilfe Breite: Zeit: Ortsschritt: Zeitschritt: Fortschrittsanzeige Euler Ortsspeicherschritt: Zeitspeicherschritt: Diffusionsterm:'+s); - if (ErrorMsg='') then begin - setlength(Ausgabedateien,0); - for Abl:=False to True do - for FI:=low(TFeldInhalt) to high(TFeldInhalt) do - if not (Abl and (FI in [fiGamma,fiP])) then - if HasOption(Copy('d',1,Byte(Abl))+OptionNamen[FI]) then begin - setlength(Ausgabedateien,length(Ausgabedateien)+1); - Ausgabedateien[length(Ausgabedateien)-1].Inhalt:=FI; - Ausgabedateien[length(Ausgabedateien)-1].Ableitung:=Abl; - AssignFile(Ausgabedateien[length(Ausgabedateien)-1].Datei,GetOptionValue(Copy('d',1,Byte(Abl))+OptionNamen[FI])); - end; - if length(Ausgabedateien)=0 then - ErrorMsg:='Du solltest irgendetwas abspeichern lassen! ('+AusgabeHilfe+')'; - end; - if ErrorMsg<>'' then begin - ShowException(Exception.Create(ErrorMsg)); - Terminate; - Exit; - end; + prot:=tProtokollant.create('error'); - // parse parameters - if HasOption('H','Hilfe') then begin - WriteHelp; - Terminate; - Exit; + if paramcount<>1 then begin + prot.schreibe('Bitte genau einen Parameter übergeben, nämlich die Parameterdatei!',true); + halt(1); end; - if HasOption('Z','Zeit') then Endzeit:=strtofloat(GetOptionValue('Z','Zeit')) - else Endzeit:=100; - if HasOption('x','Ortsschritt') then deltaX:=strtofloat(GetOptionValue('x','Ortsschritt')) - else deltaX:=1e-2; - if HasOption('B','Breite') then Breite:=round(strtofloat(GetOptionValue('B','Breite'))/deltaX) - else Breite:=1000; - if HasOption('t','Zeitschritt') then deltaT:=strtofloat(GetOptionValue('t','Zeitschritt')) - else deltaT:=deltaX/10; - if HasOption('X','Ortsspeicherschritt') then begin - s:=GetOptionValue('X','Ortsspeicherschritt'); - if pos('#',s)=0 then sDX:=strtofloat(s) - else sDX:=Breite*deltaX/strtofloat(rightstr(s,length(s)-1)); - end - else sDX:=deltaX; - if HasOption('T','Zeitspeicherschritt') then begin - s:=GetOptionValue('T','Zeitspeicherschritt'); - if pos('#',s)=0 then sDT:=strtofloat(s) - else sDT:=Endzeit/strtofloat(rightstr(s,length(s)-1)); - end - else sDT:=sDX; - if HasOption('D','Diffusionsterm') then pDNMax:=strtofloat(GetOptionValue('D','Diffusionsterm')) - else pDNMax:=-1; - if HasOption('E','Euler') then Zeitverfahren:=zfEulerVorwaerts - else Zeitverfahren:=zfRungeKuttaVier; - { add your program here } + start:=now; - Prot:=TProtokollant.create('error'); - - Prot.schreibe('Simulationsbox'); - Prot.schreibe(' '+inttostr(breite)+' Felder breit ('+floattostr(breite*deltaX)+' lambda)'); - Prot.schreibe(' '+floattostr(Endzeit)+' T lang (etwa '+inttostr(round(Endzeit/deltaT))+' Schritte)'); - Prot.schreibe(' '+floattostr(deltaT)+' Zeitschritt und '+floattostr(deltaX)+' Ortsschritt'); - for i:=0 to length(Ausgabedateien)-1 do begin - Prot.schreibe('Mache Ausgabe '+Copy('d',1,Byte(Ausgabedateien[i].Ableitung))+Ausgabenamen[Ausgabedateien[i].Inhalt]+'.'); - Rewrite(Ausgabedateien[i].Datei,1); - j:=0; - BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); - j:=floor(Endzeit/sDT+1); - BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint)); - end; - if pDNMax<0 then Prot.schreibe('Maximales px * dn/dx / n wird automatisch ermittelt.') - else Prot.schreibe('Maximales px * dn/dx / n: '+floattostr(pDNMax/deltaX)+'mc/dx'); - s:='Iterationsverfahren: '; - case Zeitverfahren of - zfEulerVorwaerts: s:=s+'Euler-Vorwärts'; - zfRungeKuttaVier: s:=s+'Runge-Kutta'; - else s:=s+'unbekannt!?'; - end{of case}; - Prot.schreibe(s); zeitDatei:=0; zeitPhysik:=0; sT:=-Min(deltaT,sDT)/2; - Gitter:=TGitter.create(Breite,deltaT,deltaX,pDNMax,@Anfangsdichte,Zeitverfahren,Prot); + simulation:=tSimulation.create(paramstr(1)); // Gitter.Al:=@InputFeld; while Gitter.t + - + - - - - - - + + + - - - - - + + + - + + + + + + + + + + + + + + + + + + - - - - - + + - - + - + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - + - - + + - - + + - - + + - + - + - + - + - - + + - + - - + + - - + + - - + + - - + + - + - + - + - + - - + + - + - + - - + + - + - + - + - + - + - + - - + + - + diff --git a/input.plap b/input.plap new file mode 100644 index 0000000..2e5ecc4 --- /dev/null +++ b/input.plap @@ -0,0 +1,49 @@ +# Parameterdatei für Plasmapropagation + +allgemein + runge-Kutta-4 + ortsschritt 10^-2 λ + zeitschritt 10^-3 T + zeit 100 T + breite 10 λ +allgemeinEnde + +ausgaben + suffix _test.dat + felder AX,AY,dAYDT,N,dPhiDX,dPsiDX,Gamma +ausgabenEnde + +!set $tFwhm: 5 T +!set $tMitte: 20 T + +lichtVonLinks 0.5 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * Sin(ω*t) + +dateiEnde + +teilchen1 + ladung -q + masse me + maximaldichte 0 + beweglich + !setze $profilbreite: 5 λ + !setze $randbreite: .1 λ + verteilung stückweise + 0 + (breite-$profilbreite)/2 - $randbreite + sin((x - (breite-$profilbreite)/2 - $randbreite)*π/2/$randbreite)^2 + (breite-$profilbreite)/2 + 1 + (breite+$profilbreite)/2 + sin((x - (breite+$profilbreite)/2 + $randbreite)*π/2/$randbreite)^2 + (breite+$profilbreite)/2 + $randbreite + 0 + stückweiseEnde +teilchen1Ende + +teilchen2 + ladung q + maximaldichte nmax1 + unbeweglich + verteilung wie teilchen1 +teilchen2Ende + -- cgit v1.2.3-70-g09d2