% coupler.m
% Andreas Franzen, Dl1YDR, 2019-11-26
%
% Das folgende Programm simuliert einen Parallelleitungsrichtkoppler. Zum
% Start unter Debian GNU/Linux das Paket "openems" installiert. Eine
% ausführliche Beschreibung ist unter "http://openems.de" zu finden. Wichtig
% ist, die Pfade in der ".octaverc" zu setzen. In der hier verwendeten
% Konfiguration sind dies die folgenden zwei Einträge:
%
% addpath('/usr/share/octave/packages/openems-0.0.35');
% addpath('/usr/share/octave/packages/csxcad-0.0.35');
%
% Unter dem oberen Pfad steht z.B. "InitFDTD.m". In dem unteren Verzeichnis
% findet man z.B. "AddBox.m".
%
% Das Programm wird gestartet mit dem Befehl "octave coupler.m".
%
% Um später die Ausgangsdaten darzustellen wird das Programm ParaView ver-
% wendet. Dazu installiert man das Paket "paraview". Um die Grafiken zu
% erzeugen wird das Kommando "graph" verwendet. Dazu installiert man das
% Paket "plotutils".

% InitFDTD initialisiert das Objekt FDTD, das die Simulationsrechnung
% beschreibt. Der Parameter NrTS (number of time steps) wird auf 3000 gesetzt.
% Dies ist die maximale Anzahl an Iterationsschritten. EndCriteria ist die
% Energie, die im simulierten Gesamtsystem unterschritten werden muss, um
% die Simulation zu beenden. 1e-12 ist sehr klein und bewirkt, dass die 3000
% Schritte tatsächlich abgearbeitet werden. OverSampling bezieht sich auf
% die Anzahl an Bildern des elektrischen Feldes, die pro Halbwelle
% abgespeichert werden.
FDTD = InitFDTD('NrTS',3000,'EndCriteria',1e-12,'OverSampling',10);
% SetGaussExcite setzt die Anregung auf ein Wellenpaket, das zunächst
% ansteigt und dann wieder abfällt. Der erste Parameter ist die
% Mittenfrequenz von 2 GHz. Der zweite Parameter beschreibt den
% Frequenzabstand, bei dem die Intensität auf -20dB gefallen ist. Praktisch
% beschreibt der Impuls den Frequenzbereich von sehr niedrigen Frequenzen
% bis 4 GHz. Aus der Impulsantwort des Systems kann dann das
% Frequenzverhalten in diesem Frequenzbereich berechnet werden.
FDTD = SetGaussExcite(FDTD,2e9,2e9);
% Mit SetBoundaryCond werden die Eigenschaften der Wände des quaderförmigen
% Simulationsraums gesetzt. PEC bedeute "perfect electrical conductor", also
% eine allseits geschlossene Kiste aus gut leitendem Metall.
FDTD = SetBoundaryCond(FDTD,{'PEC' 'PEC' 'PEC' 'PEC' 'PEC' 'PEC'});
% InitCSX initialisiert das Objekt CSX, welches die Geometrie beschreibt.
CSX = InitCSX();
% Die Struktur "mesh" wird mit den Komponenten x, y und z definiert. Diese
% enthalten Arrays, die die Position der Gitterlinien der Simulation
% beschreiben. In der Simulation gibt es eine Metallkiste von 12+224=236 mm
% Tiefe, 18 mm Höhe und 200+500=700 mm Breite. Die x-Achse zeigt in die
% Tiefe, die y-Achse in die Höhe und die z-Achse nach rechtes. "-12:224"
% bedeutet Arrayelement "-12, -11, -10, ..., 223, 224".
mesh.x = -12:224;
mesh.y = 0:18;
mesh.z = -200:500;
% DefineRectGrid bestimmt das rechtwinklige Gitter mit den Gitterlinien der
% Struktur "mesh". Der zweite Parameter legt fest, dass ein Abstand von 1
% in der Struktur "mesh" 1 mm bedeutet.
CSX = DefineRectGrid(CSX,1e-3,mesh);
% AddMetal definiert ein gut elektrisch Leitfähiges Objekt mit der
% Bezeichnung "line1". Dies ist eine Vierkantstange von links nach rechts
% von 6 mm Kantenlänge, die 2 mm vor der linken Wand und 2 mm vor der
% rechten Wand der quaderförmigen Kiste endet.
CSX = AddMetal(CSX,'line1');
% AddBox fügt dem Objekt "line1" eine quaderförmige Kiste hinzu. Der
% Parameter 10 bezieht sich auf die Reihenfolge, in der das Objekt aufgebaut
% wird, spielt hier keine Rolle, da nur elektrisch leitfähige Objekte
% additiv aufgebaut werden. Die folgenden Parameter beschreiben die
% Eckpunkte des Quaders in den Komponenten x, y und z.
CSX = AddBox(CSX,'line1',10,[6 6 -198],[12 12 498]);
% AddMetal definiert das Objekt "line2" als gut elektrisch leitfähiges
% Objekt. Dies ist die Koppelleitung, die von hinten kommt, parallel zur
% ersten Leitung verläuft und dann nach hinten zurückgeht. Die Ecken sind
% angeschnitten, um die Reflexion einer Welle an der Ecke weitgehend zu
% unterdrücken. Das Objekt wird definiert als ein Prisma mit Grundfläche in
% der x-z-Ebene.
CSX = AddMetal(CSX,'line2');
% In der Struktur "p" werden die Koordinaten der Punkte des Prismas nach z
% und x abgelegt. Der Startpunkt ist also 30 mm nach rechts und 222 mm nach
% hinten vom Koordinatenursprung aus gemessen. Dieser Punkt befindet sich
% also 2 mm entfernt von der hinteren Wand.
p(1,1)=30;p(2,1)=222;
p(1,2)=30;p(2,2)=28;
p(1,3)=40;p(2,3)=18;
p(1,4)=260;p(2,4)=18;
p(1,5)=270;p(2,5)=28;
p(1,6)=270;p(2,6)=222;
p(1,7)=264;p(2,7)=222;
p(1,8)=264;p(2,8)=24;
p(1,9)=36;p(2,9)=24;
p(1,10)=36;p(2,10)=222;
% AddLinPoly legt "line2" als geschlossenes Polygon in der z-x-Ebene fest,
% welches zu einem Prisma in y-Richtung auseinandergezogen wird. Die "10"
% ist wieder ein hier unwichtiger Parameter der Reihenfolge des Aufbaus. Die
% "1" ist die Richtung, in welche auseinandergezogen wird. "1" bedeutet
% y-Richtung. Die erste "6" ist die Position des Basispolygons in y-Richtung.
% Es befindet sich also in einer Höhe von 6 mm von der Bodenplatte. Die
% zweite "6" bedeutet, dass das Polygon um 6 mm in y-Richtung gezogen wird.
% Es endet also 6 mm unterhalb der Deckelplatte.
CSX = AddLinPoly(CSX,'line2',10,1,6,p,6);
% Mit AddLumpedPort werden zweipolige Anschlüsse an der Struktur angebracht.
% Der Parameter "5" ist wieder eine hier unwichtige Aufbaupriorität. Der
% folgende Parameter gibt die laufende Nummer des Anschlusses an. "50" ist
% der elektrische Widerstand der angeschlossen ist in Ohm. Danach folgt der
% erste Anschlusspunkt an der Wand und dann der zweite Anschlusspunkt an der
% linken Stirnseite des Leiters "line1". "[0 0 1] bedeute, dass der
% Anschluss in z-Richtung verläuft. "true" bedeutet, dass der Anschluss über
% eine Spannungsquelle angeschlossen ist, die mit dem oben definierten
% Impuls angesteuert wird. Anschluss 1 ist der Eingang.
[CSX port] = AddLumpedPort(CSX,5,1,50,[9 9 -200],[9 9 -198],[0 0 1],true);
% Anschluss 2 ist passiv. Es ist der Ausgang der Leitung "line1" auf der
% rechten Seite der Kiste. Anschluss 2 ist der direkte Ausgang.
[CSX port] = AddLumpedPort(CSX,5,2,50,[9 9 500],[9 9 498],[0 0 1],false);
% Anschluss 3 ist der Anschluss der u-förmigen Koppelleitung links oben,
% also bei einem Leitungsrichtkoppler der gekoppelte Ausgang, an welchem die
% nach rechts laufende Welle ausgekoppelt wird.
[CSX port] = AddLumpedPort(CSX,5,3,50,[224 9 33],[222 9 33],[1 0 0],false);
% Anschluss 4 ist der Anschluss rechts oben, also der ungekoppelte Ausgang,
% an dem die hier nicht absichtlich angeregte nach links laufende Welle
% ausgekoppelt wird.
[CSX port] = AddLumpedPort(CSX,5,4,50,[224 9 267],[222 9 267],[1 0 0],false);
% Um das Wandern der Impulse auf den Leitern beobachten zu können, wird das
% elektrische Feld in der Fläche, die 1 mm oberhalb der beiden Leiter liegt
% aufgezeichnet. Dazu wird zunächst mit AddDump eine solche Aufzeichnung
% definiert. "E" ist der Name der Aufzeichnung, "DumpMode" gibt an, welche
% Größe aufgezeichnet werden soll. "0" bedeutet, dass das elektrische Feld
% aufgezeichnet werden soll. 
CSX = AddDump(CSX,'E','DumpMode',0);
% Jetzt wird für die Aufzeichnung mit dem Namen "E" eine Ebene 1 mm oberhalb
% der Leiter festgelegt, die die gesamte Fläche der Metallkiste ausmacht.
CSX = AddBox(CSX,'E',0,[-12 13 -200],[224 13 500]);
% Die Simulationsdaten werden im Verzeichnis "tmp" im Heimatverzeichnis des
% Benutzers abgelegt. Dieses Verzeichnis wird gegebenenfalls erzeugt.
mkdir('~/tmp');
% Die Daten liegen hier im Unterverzeichnis "openEMS". Dies wird zunächst
% entfernt. "rmdir" ist ein Kommando von "Octave". Eine ausführliche
% Beschreibung findet man unter "https://octave.org/octave.pdf". Der
% Parameter "s" bedeutet, dass auch Unterverzeichnisse rekursiv entfernt
% werden.
rmdir('~/tmp/openEMS','s');
% Danach wird "openEMS" leer neu angelegt.
mkdir('~/tmp/openEMS');
% WriteOpenEMS erzeugt eine Beschreibung der Simulation als XML-File.
WriteOpenEMS('~/tmp/openEMS/coupler.xml',FDTD,CSX);
% Mit CSXGeomPlot wird die Struktur angezeigt. Dies dient der Kontrolle. Das
% Programm wird mit "File->exit" geschlossen.
CSXGeomPlot('~/tmp/openEMS/coupler.xml');
% Mit RunOpenEMS wird die Simulation durchgeführt.
RunOpenEMS('~/tmp/openEMS','coupler.xml');
% Die Schrittweite der Simulation ist automatisch bestimmt 1,92583 ps. Der
% symmetrische Ansteuerungsimpuls hat eine Länge von 744 Schritten. Bei 3000
% Schritten insgesamt hat das Maximum des Impulses also Zeit, um in der
% Struktur (3000-744/2)x1,92x0,3 = 1518 mm zurückzulegen. Da die Kiste eine
% Breite von 700 mm hat, kann der Reflex des Impulsmaximums am
% Ausgangsanschluss, welcher durch den nicht exakten Abschluss der Leitung
% bewirkt wird, die Kiste als am Eingang gerade wieder verlassen.

% Nach Ablauf der Simulation können die Daten mit dem Programm "ParaView"
% angezeigt werden.

f = linspace(0,4e9,501);
%port = calcPort(port,path,f)
U = ReadUI({'port_ut1','port_ut2','port_ut3','port_ut4'},'~/tmp/openEMS',f);

t = U.TD{1}.t./1e-9;
ut1 = U.TD{1}.val./1e-3;
ut2 = U.TD{2}.val./1e-3;
ut3 = U.TD{3}.val./1e-3;
ut4 = U.TD{4}.val./1e-3;
AUS = fopen('plot_t.dat','w');
fprintf(AUS,'%e %e\n',[t;ut1]);
fprintf(AUS,'\n');
fprintf(AUS,'%e %e\n',[t;ut2]);
fprintf(AUS,'\n');
fprintf(AUS,'%e %e\n',[t;ut3]);
fprintf(AUS,'\n');
fprintf(AUS,'%e %e\n',[t;ut4]);
fclose(AUS);

uf1 = abs(U.FD{1}.val);
uf2 = abs(U.FD{2}.val);
uf3 = abs(U.FD{3}.val);
uf4 = abs(U.FD{4}.val);
V21 = uf2./uf1;
V31 = uf3./uf1;
V41 = uf4./uf1;
f2 = f./1e9;
AUS = fopen('plot_f.dat','w');
fprintf(AUS,'%e %e\n',[f2;V21]);
fprintf(AUS,'\n');
fprintf(AUS,'%e %e\n',[f2;V31]);
fprintf(AUS,'\n');
fprintf(AUS,'%e %e\n',[f2;V41]);
fclose(AUS);
