path: root/Plasmapropagation.lpr

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.