Forum: HF, Funk und Felder Simulation EM-Feld (Matlab)


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
von M. M. (blackcow)


Bewertung
0 lesenswert
nicht lesenswert
Ich wollte mit matlab eine 1-Dimensionale EM-Welle Simulieren. Mein 
Ausgangspunkt sind die Maxwellschen Gleichungen (zur Vereinfachung habe 
ich E=D und B=H gesetzt:
rot(E)=dB/dt
rot(B)=dE/dt

Bei einer Dimension habe ich nun in Matlab sechs 1-Dimensionale Arrays 
erstellt:
buf_rotE = aktuelle Rotation von E-Feld
buf_rotB = aktuelle Rotation von B-Feld
buf_E =    vergangene&aktuelle E-Feld-Vektoren
buf_B =    vergangene&aktuelle B-Feld-Vektoren

Ein Simulationsschritt sieht bei mir folgendermaßen aus:
Aus buf_rotE und buf_rotB werden jeweils das aktuelle E- und B-Feld 
berechnet. Mit diesen Werten und den des Schrittes davor werden dann die 
Differenzen (sprich zeitliche Ableitung) berechnet und als neue 
Rotationswerte in buf_rotE, bzw. buf_rotB eingetragen.

Das Problem ist, es funktioniert nicht. Das Ergebnis sieht überhaupt 
nicht nach Welle aus. Kann mir jemand sagen wo der Haken liegt?
1
% Initialisiere Variablen
2
    % Einstellungen
3
        % Feldgröße
4
        X_width = 200;
5
    % Simulationsgeschwindigkeit
6
    d_t  = 0.01;
7
    % Verknüpfung rotE <- dB/dt
8
    k_rotE_dB = 1;
9
    % Verknüpfung rotB <- dE/dt
10
    k_rotB_dE = 1;
11
12
  buf_rotE = zeros(X_width,1);
13
  buf_rotB = zeros(X_width,1);
14
  buf_E = zeros(X_width,2);    % X, alt/neu
15
  buf_B = zeros(X_width,2);       % X, alt/neu
16
  i_alt = 1;
17
  i_neu = 2;
18
  gw_rotE  = zeros(2, 1);          % min/max
19
  gw_offset = 0;
20
  gw_skng = 0;
21
  X = 0;
22
    rain = 0;
23
    
24
buf_rotE(X_width/2) = 250;
25
% Führe Simulationsschritt durch
26
for step = 1:50
27
% buf_rotE(X_width/2) = 250*sin(step*0.4);
28
    figure(1);
29
    subplot(4,1,1);
30
        bar(buf_rotE);
31
    title(['step: ' num2str(step)]);
32
        axis([1 X_width -250 250]);
33
    subplot(4,1,3);
34
        bar(buf_rotB);
35
        axis([1 X_width -250 250]);
36
        title('rotB');
37
    pause(d_t);
38
  % Rotiere Array Index: alt <=> neu
39
        if i_alt==1
40
            i_alt = 2;
41
      i_neu = 1;
42
        else
43
            i_alt = 1;
44
            i_neu = 2;
45
        end
46
  % Berechne Felder und schreibe min/max Werte mit
47
        gw_rotE(1) = 0;
48
        gw_rotE(2) = 0;
49
        gw_B(1) = 0;
50
        gw_B(2) = 0;
51
        for X = 1:X_width
52
            % aktuelles E-Feld
53
                buf_E(X, i_neu) = 0;
54
                if X>=X_width
55
                    buf_E(X, i_neu) = 0; %buf_rotE(X)-buf_rotE(X-1);
56
                elseif X==1
57
                    buf_E(X, i_neu) = 0; %buf_rotE(X+1)-buf_rotE(X);
58
                else
59
                    buf_E(X, i_neu) = buf_E(X, i_alt) + buf_rotE(X)-buf_rotE(X-1);
60
                    buf_E(X, i_neu) = buf_rotE(X+1)-buf_rotE(X) + buf_E(X, i_neu);
61
                end
62
                % aktuelles B-Feld
63
                    buf_B(X, i_neu) = 0;
64
                    if X>=X_width
65
                        buf_B(X, i_neu) = 0; %buf_rotB(X)-buf_rotB(X-1);
66
                    elseif X==1
67
                        buf_B(X, i_neu) = 0; %buf_rotB(X+1)-buf_rotB(X);
68
                    else
69
                        buf_B(X, i_neu) = buf_B(X, i_neu) + buf_rotB(X+1)-buf_rotB(X);
70
                        buf_B(X, i_neu) = buf_rotB(X)-buf_rotB(X-1) + buf_B(X, i_neu);
71
                    end
72
        end
73
        for X = 1:X_width
74
                % aktuelles E-Rotation
75
                    buf_rotE(X) = ( buf_B(X, i_neu)-buf_B(X, i_alt) ) * k_rotE_dB;
76
                    % schreibe min/max Werte mit
77
                    if buf_rotE(X)<gw_rotE(1)
78
                        gw_rotE(1) = buf_rotE(X);
79
                    end
80
                    if buf_rotE(X)>gw_rotE(2)
81
                        gw_rotE(2) = buf_rotE(X);
82
                    end
83
                % aktuelles B-Rotation
84
                    buf_rotB(X) = ( buf_E(X, i_neu)-buf_E(X, i_alt) ) * k_rotB_dE;
85
        end
86
87
    subplot(4,1,2);
88
        bar(buf_E(:,i_neu));
89
        axis([1 X_width -250 250]);
90
        title('E');
91
    subplot(4,1,4);
92
        bar(buf_B(:,i_neu));
93
        axis([1 X_width -250 250]);
94
        title('B');
95
    hold off;
96
    % pause;
97
end

von Rolf R. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Ich verstehe nichts von Matlab.

Aber warum nimmst du nicht die Wellengleichung, die man aus den 
Maxwellschen Gleichungen erhält, und löst diese zuerst mathematisch auf 
dem Papier für deinen Anwendungsfall. Mit dem Ergebnis hast du dann dein 
E- und dein B-Feld in Abhängigkeit von Zeit und Ort.

Damit würde ich dann in Matlab weiter machen.

Das Einsetzten der Maxwellschen Gleichungen ineinander und die Umformung 
zur Wellengleichung und die Lösung dieser ist bei mir aber schon 
Jahrzehnte her. Ich müsste mich da auch wieder einarbeiten.

von Rolf R. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Oder willst du FDTD wie hier

https://www.youtube.com/watch?v=RoUrqR71UWo

machen?

von Bolle (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Matlab ist dafür nicht so geeignet. Ich empfehle Comsol.

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.