MATLAB Answers

0

Why do I not get the whole domain plotted?

Asked by Matthias Tiek on 17 Aug 2019
Latest activity Commented on by Walter Roberson
on 18 Aug 2019
I want to get a plot like this:
attenuation_wolf_graph.JPG
but all I get is this:
test_wrong.jpg
Besides the missing resonaces, moreover, I can not plot the whole frequency domain as I want to. I want to plot from1E6[Hz] up to 20E9[Hz]
My code so far:
% Testumgebung
clear all
clc
% Dimensionen
r_0 = 10E-3; % [m] test
d = 1E-3; % [m] test
% Frequenz und Winkelgeschwindigkeit
f1 = 1E6; % [Hz] Frequenz zum Starten
df = 1E6; % [Hz] Schrittweite der Frequenz
fmax = 20E9; % [Hz] letzte Frequenz
f = f1:df:fmax; % [Hz] Frequenz als Vektor ausgeben
omega = 2*pi*f; % [rad/s] Winkelgeschwindigkeit
% mag. Permeabilität
mu0 = 1.256637061E-6; % [N/A^2] mag. Feldkonstante
mu_r = 1; % realtive mag. Permeabilität
mu = mu0*mu_r; % [N/A^2] effektive mag. Permeabilität
% elektrische Leitfähigkeit
sigma_Al = 35.7E6; % [S/m] Leitfähigkeit von Aluminium
sigma_Co = 57.1E6; % [S/m] Leitfähigkeit von Kupfer
% Wirbelstromkonstante
k_w_Al = sqrt(1i*omega*mu*sigma_Al); % Wirbelstromkonstante für Aluminium (3.24)
k_w_Co = sqrt(1i*omega*mu*sigma_Co); % Wirbelstromkonstante für Kupfer (3.24)
c_0 = 299792458; % Lichtgeschwindigkeit im Vakuum [m/s]
lambda_0 = c_0./f; % Wellenlänge [m]
k_0 = 2*pi./lambda_0; % Wellenzahl [rad/m]
eps0 = 8.854187817E-12; % Permitivität des Vakuums [As/Vm]
Z_0 = sqrt(mu0/eps0); % Wellenwiderstand des Vakuums [Ohm]
E_a = 10000; % anregende el. externe Feld [V/m]
H_a = E_a/Z_0; % anregende mag. externe Feld [A/m]
S_a = (E_a^2)/(2*Z_0); % Wirkleistungsdichte [W/m^2] nach Gustrau HF-Technik (2.104)
kr =k_0*r_0; % Abkürzung von k_0*r_0
%Besselfunktionen
J_1 = besselj(1,kr); % Besselfunktionswert der Ordnung 1
Ha_1 = besselh(1,2,kr); % Besselfunktionswert 3. Gattung, 1-ter Ordnung, zweiter Art
% Schirmfaktor
Q_Al = 2./(k_w_Al.*r_0.*sinh(k_w_Al.*d));
% Dämpfung
a_m_Al_wave = -log(abs(Q_Al)) + log(pi.*abs(Ha_1.*J_1)); % E-Feld ist longitudinal zur z-Achse!!!!
semilogx(f,a_m_Al_wave)
xlabel('f [Hz]')
ylabel('a_m [dB]')
legend({'Aluminium ~'})
title('Dämpfung des magnetischen Feldes')
grid on

  0 Comments

Sign in to comment.

2 Answers

Answer by Walter Roberson
on 17 Aug 2019
 Accepted Answer

You do not get a plot up to 20E9 because your calculation overflows to infinity.
The below does not overflow to infinity, but it is noticeably slower than the original code.
% Testumgebung
Q = @(v) sym(v);
Pi = sym('pi');
% Dimensionen
r_0 = Q(10E-3); % [m] test
d = Q(1E-3); % [m] test
% Frequenz und Winkelgeschwindigkeit
f1 = Q(1E6); % [Hz] Frequenz zum Starten
df = Q(1E6); % [Hz] Schrittweite der Frequenz
fmax = Q(20E9); % [Hz] letzte Frequenz
f = f1:df:fmax; % [Hz] Frequenz als Vektor ausgeben
omega = 2*Pi*f; % [rad/s] Winkelgeschwindigkeit
% mag. Permeabilität
mu0 = Q(1.256637061E-6); % [N/A^2] mag. Feldkonstante
mu_r = Q(1); % realtive mag. Permeabilität
mu = mu0*mu_r; % [N/A^2] effektive mag. Permeabilität
% elektrische Leitfähigkeit
sigma_Al = Q(35.7E6); % [S/m] Leitfähigkeit von Aluminium
sigma_Co = Q(57.1E6); % [S/m] Leitfähigkeit von Kupfer
% Wirbelstromkonstante
k_w_Al = sqrt(1i*omega*mu*sigma_Al); % Wirbelstromkonstante für Aluminium (3.24)
k_w_Co = sqrt(1i*omega*mu*sigma_Co); % Wirbelstromkonstante für Kupfer (3.24)
c_0 = Q(299792458); % Lichtgeschwindigkeit im Vakuum [m/s]
lambda_0 = c_0./f; % Wellenlänge [m]
k_0 = 2*Pi./lambda_0; % Wellenzahl [rad/m]
eps0 = Q(8.854187817E-12); % Permitivität des Vakuums [As/Vm]
Z_0 = sqrt(mu0/eps0); % Wellenwiderstand des Vakuums [Ohm]
E_a = Q(10000); % anregende el. externe Feld [V/m]
H_a = E_a/Z_0; % anregende mag. externe Feld [A/m]
S_a = (E_a^2)/(2*Z_0); % Wirkleistungsdichte [W/m^2] nach Gustrau HF-Technik (2.104)
kr =k_0*r_0; % Abkürzung von k_0*r_0
%Besselfunktionen
J_1 = besselj(1,kr); % Besselfunktionswert der Ordnung 1
Ha_1 = besselh(1,2,kr); % Besselfunktionswert 3. Gattung, 1-ter Ordnung, zweiter Art
% Schirmfaktor
Q_Al = 2./(k_w_Al.*r_0.*sinh(k_w_Al.*d));
% Dämpfung
a_m_Al_wave = -log(abs(Q_Al)) + log(Pi.*abs(Ha_1.*J_1)); % E-Feld ist longitudinal zur z-Achse!!!!
semilogx(f,a_m_Al_wave)
xlabel('f [Hz]')
ylabel('a_m [dB]')
legend({'Aluminium ~'})
title('Dämpfung des magnetischen Feldes')
grid on

  3 Comments

Thank you. Unfortunately, your code is not comprehensible and as you have already mantioned, it runs significantly longer. Therefor, it is not applicable. I am not sure, if your suggestion about the overflow ist correct. At least the figure should plot the right range, right?
I took nearly every numeric constant and enclosed it in a call to Q(). is the symbol that is used to indicate the field of rationals. That is, I convert all of the numeric constants to rationals. This is done in the symbolic toolbox: by default when you apply sym() to a floating point number, the result is the same as converting it to rational.
I also changed all pi to Pi, where Pi has been defined as symbolic π. Since π is not rational, it would have been inappropriate for me to have converted pi to rational.
Those are the only changes I made -- hardly incomprehensible.
The effect is to process all of the computations inside the Symbolic Toolbox, using exact solutions where possible and using higher precision calculations (by default, 32 decimal digits) when exact solutions are not possible.
I am not sure, if your suggestion about the overflow ist correct.
In the line after your clc, add the line of code
dbstop if naninf
and run the result. The code will stop on line 37,
Q_Al = 2./(k_w_Al.*r_0.*sinh(k_w_Al.*d));
You will find that sinh(k_w_Al.*d) overflows to double precision infinity at k_w_Al(3582) when it becomes 710520.752519461 + 710520.752519461i which makes sense when you look at asinh(realmax)
Calculating backwards,
syms OM
vpa(solve(sqrt(1i*OM*mu*sigma_Al)==(1+1i)*asinh(realmax)./d)./(2*pi))
This gives you f of 3581547375.1913578738860940474962, which is roughly 3.582E9, which is well inside your frequency range.
Therefore, your calculations can and do overflow double precision.
Your largest frequency generates sinh() values in the range of 1e728 -- more than the square of the maximum double precision number.
At least the figure should plot the right range, right?
No. You do not have any xlim() or equivalent forcing the entire 20E9 to be plotted. The automatic limits are for a "nice" value just beyond the maximum finite value plotted. With everything after 3.581E9 giving inf results, that is as far as it is going to plot.
When you ask to plot with a coordinate that is infinite or nan, MATLAB does not plot that point.
At larger values of X, log(sinh(X)) is approximately X, so with k_w_Al(end).*d approaching 1678*(1+i), then log(sinh()) is about 1678, but roughly 709 is as far as you can go in double precision without overflowing.

Sign in to comment.


Answer by KALYAN ACHARJYA on 17 Aug 2019
Edited by KALYAN ACHARJYA on 17 Aug 2019

Why do I not get the whole domain plotted?
Because: In the first figure, there are three plots, but in your matlab code having single plot.

  3 Comments

Thank you. I see the three plotted lines in the first picture. If you refer to the resonances, two of these lines have resonances and due to the calculation in my code there should also be some resonances. BUT, the more important part ist, why is my plot only from 1E6[Hz] up to 3.5E9[Hz] AND NOT up to 20E9[Hz]?
There are more scale range differences, x scal range started 10^6, original image 10^7??
12356.png
I don't think are any major techical issues, just assign the same numerical data, as per figure 1.
Thank you. I do not want to copy the original image. I need a different range. I want the frequencies from 1E6 [Hz] up to 20E9[Hz]. The figure should look similary but not exactly the same. Look in my code, I specified the frequency:
% Frequenz und Winkelgeschwindigkeit
f1 = 1E6; % [Hz] Frequenz zum Starten
df = 1E6; % [Hz] Schrittweite der Frequenz
fmax = 20E9; % [Hz] letzte Frequenz
f = f1:df:fmax; % [Hz] Frequenz als Vektor ausgeben
And I plot my values over f:
semilogx(f,a_m_Al_wave)
But Matlab does not use all specified frequencies.

Sign in to comment.