diff options
Diffstat (limited to 'Plasmapropagation.lpr')
| -rw-r--r-- | Plasmapropagation.lpr | 279 |
1 files changed, 279 insertions, 0 deletions
diff --git a/Plasmapropagation.lpr b/Plasmapropagation.lpr new file mode 100644 index 0000000..52e1926 --- /dev/null +++ b/Plasmapropagation.lpr @@ -0,0 +1,279 @@ +program Plasmapropagation; + +{$mode objfpc}{$H+} + +uses + {$IFDEF UNIX}{$IFDEF UseCThreads} + cthreads, + {$ENDIF}{$ENDIF} + Classes, SysUtils, CustApp, + { you can add units after this } + math, Physikunit, crt; + +var ErrorCode: longint; + +type + + { TPlasmapropagation } + + TPlasmapropagation = class(TCustomApplication) + protected + procedure DoRun; override; + 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; +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; + + // parse parameters + if HasOption('H','Hilfe') then begin + WriteHelp; + Terminate; + Exit; + 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 } + + 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); +// Gitter.Al:=@InputFeld; + while Gitter.t<Endzeit do begin + if HasOption('F','Fortschrittsanzeige') then begin + if errorCode<2 then + s:=Gitter.gibErhaltungsgroessen; + if (floor(100*Gitter.t/Endzeit) < floor(100*(Gitter.t+deltaT)/Endzeit)) or keypressed then begin + if keypressed then c:=readkey + else c:=#0; + case c of + #27,'q': begin + errorCode:=3; + break; + end; + ' ': begin + writeln(' ... Pause (beliebige Taste drücken um fortzufahren) ...'); + readkey; + writeln(' ... weiter geht''s ...'); + end; + else begin + Prot.schreibe(inttostr(round(100*Gitter.t/Endzeit))+'% (t='+floattostr(Gitter.t)+'T)',true); + Prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true); + Prot.schreibe('ETA: '+timetostr((now-start)*(Endzeit-Gitter.t)/max(Gitter.t,deltaT)),true); + Prot.schreibe('aktueller Zeitschritt: '+floattostr(deltaT)+'T',true); + Prot.schreibe(s); + end; + end{of case}; + end; + end; + + zeitPhysik:=zeitPhysik-now; + if errorCode<2 then + Gitter.iteriereSchritt(deltaT); + zeitPhysik:=zeitPhysik+now; + zeitDatei:=zeitDatei-now; + while Gitter.t>=sT do begin + sT:=sT+sDT; + for j:=0 to length(Ausgabedateien)-1 do + Gitter.macheAusgabe(Ausgabedateien[j],sDX); + end; + zeitDatei:=zeitDatei+now; + end; + + case errorCode of + 2: Prot.schreibe('Simulation wurde auf halbem Wege abgebrochen wegen Überlaufs.',true); + 3: Prot.schreibe('Simulation wurde auf halbem Wege abgebrochen wegen Benutzereingriffs.',true); + end{of case}; + + for i:=0 to length(Ausgabedateien)-1 do + CloseFile(Ausgabedateien[i].Datei); + + Prot.schreibe('fertig!',true); + s:=timetostr(now-start); + t:=timetostr(zeitDatei); + u:=timetostr(zeitPhysik); + while length(s)<max(length(t),length(u)) do s:=' '+s; + while length(t)<max(length(s),length(u)) do t:=' '+t; + while length(u)<max(length(s),length(t)) do u:=' '+u; + Prot.schreibe('Das hat '+s+' gedauert,',true); + Prot.schreibe(' davon '+t+' für Dateizugriffe',true); + Prot.schreibe('und nur '+u+' für die eigentliche Physik!',true); + + Gitter.Free; + Prot.Free; + // stop program loop + Terminate; + if errorCode=1 then errorCode:=0; +end; + +constructor TPlasmapropagation.Create(TheOwner: TComponent); +begin + inherited Create(TheOwner); + errorCode:=1; + StopOnException:=True; +end; + +destructor TPlasmapropagation.Destroy; +begin + inherited Destroy; +end; + +procedure TPlasmapropagation.WriteHelp; +begin + { add your help code here } + writeln('Verwendung:'); + writeln(ExeName,' -HBZxtFEP '+AusgabeHilfe); + writeln(' -H/--Hilfe Hilfe anzeigen'); + writeln(' -B/--Breite=x'); + writeln(' Breite der Simulationsbox in Plasmawellenlängen'); + writeln(' -Z/Zeit=x'); + writeln(' Länge der Simulationsbox in Plasmaoszillationsperioden'); + writeln(' -x/Ortsschritt=x'); + writeln(' Ortsdiskretisierung in Plasmaoszillationsvakuumwellenlängen'); + writeln(' -t/Zeitschritt=x'); + writeln(' Zeitdiskretisierung in Plasmaoszillationsperioden'); + writeln(' -X/Ortsspeicherschritt=x'); + writeln(' Ortsdiskretisierung der gespeicherten Werte in Plasmaoszillationsvakuumwellenlängen'); + writeln(' -T/Zeitspeicherschritt=x'); + writeln(' Zeitdiskretisierung der gespeicherten Werte in Plasmaoszillationsperioden'); + writeln(' -F/Fortschrittsanzeige'); + writeln(' -E/Euler'); + writeln(' verwende Eulerverfahren statt Runge-Kutta Stufe 4'); + writeln(' -D/Diffusionsterm=x'); + writeln(' p_x * dn/dx / n, bis zu dem numerische Stabilität gewähleistet ist.'); + writeln(' Ein hoher Wert bringt in Verbindung mit einer schlechten Ortsauflösung unphysikalische Ergebnisse!'); + writeln(' Ein negativer Wert sowie keine Angabe aktivieren die automatische Bestimmung. Das garantiert numerische Stabilität, führt jedoch möglicherweise unbemerkt zu unphysikalischen Ergebnissen.'); +end; + +var + Application: TPlasmapropagation; +begin + Application:=TPlasmapropagation.Create(nil); + Application.Run; + Application.Free; + Halt(errorCode); +end. + |