Dear all,
I am currently trying to create a custom component block for Simscape Thermal Liquid that models the heat transfer behaviour of a specific heat exchanger between two different thermal liquids. The code creates the block with no errors, nor warnings, and runs fine as long as I do not define thermal properties for both liquids. When I do define thermal properties for both liquids, I get the following error:

I imagine that the error comes from the fact that my code "expects" the 4 nodes (in/out for both liquids) to share properties, but I cannot find where that definition comes from…
Does anyone know how to define nodes that do not share fluid properties?

As a simpler test case, I am trying to make it work with a custom component block based on the Flow Rate Sensor (TL) library block, with 4 ports instead of 2. Same error comes up, but the code is much simpler to work with:
component flow_sensor_two_fluids
nodes
A = foundation.thermal_liquid.thermal_liquid; % A :left
B = foundation.thermal_liquid.thermal_liquid; % B :right
C = foundation.thermal_liquid.thermal_liquid; % C :left
D = foundation.thermal_liquid.thermal_liquid; % D :right
end
outputs
M1 = {0, ‘kg/s’};
M2 = {0, ‘kg/s’};
end
annotations
[M1] : Side=top
[M2] : Side=bottom
end
variables (Access = protected, ExternalAccess = none)
mdot1 = {0, ‘kg/s’}; % Mass flow rate from port A to port B
Phi1 = {0, ‘kW’ }; % Energy flow rate from port A to port B
mdot2 = {0, ‘kg/s’}; % Mass flow rate from port C to port D
Phi2 = {0, ‘kW’ }; % Energy flow rate from port C to port D
end
branches
mdot1 : A.mdot -> B.mdot;
Phi1 : A.Phi -> B.Phi;
mdot2 : C.mdot -> D.mdot;
Phi2 : C.Phi -> D.Phi;
end
equations
M1 == mdot1;
A.p == B.p;
A.T == B.T;
M2 == mdot2;
C.p == D.p;
C.T == D.T;
end
end

Thank you for your help!
PS: (I am aware that there already exists a heat exchanger (TL-TL) block that MathWorks provide, but would prefer to create my own)Dear all,
I am currently trying to create a custom component block for Simscape Thermal Liquid that models the heat transfer behaviour of a specific heat exchanger between two different thermal liquids. The code creates the block with no errors, nor warnings, and runs fine as long as I do not define thermal properties for both liquids. When I do define thermal properties for both liquids, I get the following error:

I imagine that the error comes from the fact that my code "expects" the 4 nodes (in/out for both liquids) to share properties, but I cannot find where that definition comes from…
Does anyone know how to define nodes that do not share fluid properties?

As a simpler test case, I am trying to make it work with a custom component block based on the Flow Rate Sensor (TL) library block, with 4 ports instead of 2. Same error comes up, but the code is much simpler to work with:
component flow_sensor_two_fluids
nodes
A = foundation.thermal_liquid.thermal_liquid; % A :left
B = foundation.thermal_liquid.thermal_liquid; % B :right
C = foundation.thermal_liquid.thermal_liquid; % C :left
D = foundation.thermal_liquid.thermal_liquid; % D :right
end
outputs
M1 = {0, ‘kg/s’};
M2 = {0, ‘kg/s’};
end
annotations
[M1] : Side=top
[M2] : Side=bottom
end
variables (Access = protected, ExternalAccess = none)
mdot1 = {0, ‘kg/s’}; % Mass flow rate from port A to port B
Phi1 = {0, ‘kW’ }; % Energy flow rate from port A to port B
mdot2 = {0, ‘kg/s’}; % Mass flow rate from port C to port D
Phi2 = {0, ‘kW’ }; % Energy flow rate from port C to port D
end
branches
mdot1 : A.mdot -> B.mdot;
Phi1 : A.Phi -> B.Phi;
mdot2 : C.mdot -> D.mdot;
Phi2 : C.Phi -> D.Phi;
end
equations
M1 == mdot1;
A.p == B.p;
A.T == B.T;
M2 == mdot2;
C.p == D.p;
C.T == D.T;
end
end

Thank you for your help!
PS: (I am aware that there already exists a heat exchanger (TL-TL) block that MathWorks provide, but would prefer to create my own) Dear all,
I am currently trying to create a custom component block for Simscape Thermal Liquid that models the heat transfer behaviour of a specific heat exchanger between two different thermal liquids. The code creates the block with no errors, nor warnings, and runs fine as long as I do not define thermal properties for both liquids. When I do define thermal properties for both liquids, I get the following error:

I imagine that the error comes from the fact that my code "expects" the 4 nodes (in/out for both liquids) to share properties, but I cannot find where that definition comes from…
Does anyone know how to define nodes that do not share fluid properties?

As a simpler test case, I am trying to make it work with a custom component block based on the Flow Rate Sensor (TL) library block, with 4 ports instead of 2. Same error comes up, but the code is much simpler to work with:
component flow_sensor_two_fluids
nodes
A = foundation.thermal_liquid.thermal_liquid; % A :left
B = foundation.thermal_liquid.thermal_liquid; % B :right
C = foundation.thermal_liquid.thermal_liquid; % C :left
D = foundation.thermal_liquid.thermal_liquid; % D :right
end
outputs
M1 = {0, ‘kg/s’};
M2 = {0, ‘kg/s’};
end
annotations
[M1] : Side=top
[M2] : Side=bottom
end
variables (Access = protected, ExternalAccess = none)
mdot1 = {0, ‘kg/s’}; % Mass flow rate from port A to port B
Phi1 = {0, ‘kW’ }; % Energy flow rate from port A to port B
mdot2 = {0, ‘kg/s’}; % Mass flow rate from port C to port D
Phi2 = {0, ‘kW’ }; % Energy flow rate from port C to port D
end
branches
mdot1 : A.mdot -> B.mdot;
Phi1 : A.Phi -> B.Phi;
mdot2 : C.mdot -> D.mdot;
Phi2 : C.Phi -> D.Phi;
end
equations
M1 == mdot1;
A.p == B.p;
A.T == B.T;
M2 == mdot2;
C.p == D.p;
C.T == D.T;
end
end

Thank you for your help!
PS: (I am aware that there already exists a heat exchanger (TL-TL) block that MathWorks provide, but would prefer to create my own) simscape, custom components, themal liquid MATLAB Answers — New Questions

​

I want to make curve fitting to these points but when i made it the blue straight line appeared. What does it mean ?I want to make curve fitting to these points but when i made it the blue straight line appeared. What does it mean ? I want to make curve fitting to these points but when i made it the blue straight line appeared. What does it mean ? matlab, curve fitting MATLAB Answers — New Questions

​

Hello,
We have Matlab R2023a installed on Oracle Linux8. Below is the message user gets. He reports that he would get a pop up Window that says "Oops Were sorry, it looks like a problem occurred." He says that when this happens matlab doesn’t crash, it marks the beginning of the GUI glitches. MatLab turns all white, And I get a ldd error in the console that spawned matlab, it was related to loading/drawing from what he remembers
We had MATE desktop and that was not working, it was worse with glitches according to the user.

id 6f789c523fa312a7f9199f39a72694e68987d191
reason: traps: MATLABWindow[3569572] general protection fault ip:7f76175fabba sp:7ffe77621cf0 error:0 in libgtk-x11-2.0.so.0.2400.31[7f7617419000+260000]

Can anyone tell me if they have seen this or can assist on how to resolve?
thanks,
DennisHello,
We have Matlab R2023a installed on Oracle Linux8. Below is the message user gets. He reports that he would get a pop up Window that says "Oops Were sorry, it looks like a problem occurred." He says that when this happens matlab doesn’t crash, it marks the beginning of the GUI glitches. MatLab turns all white, And I get a ldd error in the console that spawned matlab, it was related to loading/drawing from what he remembers
We had MATE desktop and that was not working, it was worse with glitches according to the user.

id 6f789c523fa312a7f9199f39a72694e68987d191
reason: traps: MATLABWindow[3569572] general protection fault ip:7f76175fabba sp:7ffe77621cf0 error:0 in libgtk-x11-2.0.so.0.2400.31[7f7617419000+260000]

Can anyone tell me if they have seen this or can assist on how to resolve?
thanks,
Dennis Hello,
We have Matlab R2023a installed on Oracle Linux8. Below is the message user gets. He reports that he would get a pop up Window that says "Oops Were sorry, it looks like a problem occurred." He says that when this happens matlab doesn’t crash, it marks the beginning of the GUI glitches. MatLab turns all white, And I get a ldd error in the console that spawned matlab, it was related to loading/drawing from what he remembers
We had MATE desktop and that was not working, it was worse with glitches according to the user.

id 6f789c523fa312a7f9199f39a72694e68987d191
reason: traps: MATLABWindow[3569572] general protection fault ip:7f76175fabba sp:7ffe77621cf0 error:0 in libgtk-x11-2.0.so.0.2400.31[7f7617419000+260000]

Can anyone tell me if they have seen this or can assist on how to resolve?
thanks,
Dennis general protection fault, ibgtk-x11-2.0.so.0.2400.31 MATLAB Answers — New Questions

​

I am attempting to create a whirl flutter diagram based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system. As my coefficients are periodic
Is my approach for determining the Stiffnesses and angles correct, and how can I change the process using the Floquet technique and the Mathieu equation? Because ultimately, the boundries are not shown correctly according to the reference.
clc;
clear all;

% Define parameters
N = 2; % Number of blades

I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b

C_thetaoverI_b = 0.00; % Damping coefficient over I_b
C_psioverI_b = 0.00; % Damping coefficient over I_b

h = 0.3; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 1000; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

% Frequency ranges
f_pitch= 0.001:0.5:5;
f_yaw= 0.001:0.5:5;

% Time periods for pitch and yaw
T_pitch = 1 ./ (2 * pi * f_pitch);
T_yaw = 1 ./ (2 * pi * f_yaw);

divergence_map = [];
Rdivergence_map = [];
unstable = [];

% Modify the loop to iterate over time points
for i = 1:length(T_pitch)
for j = 1:length(T_yaw)

T = max(T_pitch(i), T_yaw(j)); % Use the maximum period to cover all dynamics
t_steps = linspace(0, T, 100); % Time steps within one period

for t = t_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi / T_pitch(i);
w_omega_yaw = 2 * pi / T_yaw(j);

K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;

% Calculate matrices at time t using harmonic motion expressions
phi = 2 * pi * t / T; % Phase variation over the period

% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];

% Define damping matrix [D]
D11 = h^2*gamma*H_u*(1 – cos(2*phi)) – gamma*M_b*(1 + cos(2*phi)) – (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) – 2*(1 + cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 – cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D22 = h^2*gamma*H_u*(1 + cos(2*phi)) – gamma*M_b*(1 – cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);

D_matrix = [D11, D12;
D21, D22];

% Define stiffness matrix [K]
K11 = K_theta – h*gamma*V*H_u*(1 – cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) – gamma*V*M_u*(1 – cos(2*phi));
K22 = K_psi – h*gamma*V*H_u*(1 + cos(2*phi)) – gamma*V*M_u*sin(2*phi);

K_matrix = [K11, K12;
K21, K22];

% Compute the system matrices

M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);

P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;

P = roots([P0, P1, P2, P3, P4]);

r = 1 * P;

%Flutter
for k = 1:length(r)

if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];

% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 1e-5
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];

end
end
end

end

%Divergence

if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; t, K_psi, K_theta];
end

end
end
end

% Plotting section

figure;
hold on;

scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);

xlabel(‘K_psi’);
ylabel(‘K_theta’);

title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);

hold off;I am attempting to create a whirl flutter diagram based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system. As my coefficients are periodic
Is my approach for determining the Stiffnesses and angles correct, and how can I change the process using the Floquet technique and the Mathieu equation? Because ultimately, the boundries are not shown correctly according to the reference.
clc;
clear all;

% Define parameters
N = 2; % Number of blades

I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b

C_thetaoverI_b = 0.00; % Damping coefficient over I_b
C_psioverI_b = 0.00; % Damping coefficient over I_b

h = 0.3; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 1000; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

% Frequency ranges
f_pitch= 0.001:0.5:5;
f_yaw= 0.001:0.5:5;

% Time periods for pitch and yaw
T_pitch = 1 ./ (2 * pi * f_pitch);
T_yaw = 1 ./ (2 * pi * f_yaw);

divergence_map = [];
Rdivergence_map = [];
unstable = [];

% Modify the loop to iterate over time points
for i = 1:length(T_pitch)
for j = 1:length(T_yaw)

T = max(T_pitch(i), T_yaw(j)); % Use the maximum period to cover all dynamics
t_steps = linspace(0, T, 100); % Time steps within one period

for t = t_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi / T_pitch(i);
w_omega_yaw = 2 * pi / T_yaw(j);

K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;

% Calculate matrices at time t using harmonic motion expressions
phi = 2 * pi * t / T; % Phase variation over the period

% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];

% Define damping matrix [D]
D11 = h^2*gamma*H_u*(1 – cos(2*phi)) – gamma*M_b*(1 + cos(2*phi)) – (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) – 2*(1 + cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 – cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D22 = h^2*gamma*H_u*(1 + cos(2*phi)) – gamma*M_b*(1 – cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);

D_matrix = [D11, D12;
D21, D22];

% Define stiffness matrix [K]
K11 = K_theta – h*gamma*V*H_u*(1 – cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) – gamma*V*M_u*(1 – cos(2*phi));
K22 = K_psi – h*gamma*V*H_u*(1 + cos(2*phi)) – gamma*V*M_u*sin(2*phi);

K_matrix = [K11, K12;
K21, K22];

% Compute the system matrices

M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);

P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;

P = roots([P0, P1, P2, P3, P4]);

r = 1 * P;

%Flutter
for k = 1:length(r)

if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];

% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 1e-5
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];

end
end
end

end

%Divergence

if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; t, K_psi, K_theta];
end

end
end
end

% Plotting section

figure;
hold on;

scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);

xlabel(‘K_psi’);
ylabel(‘K_theta’);

title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);

hold off; I am attempting to create a whirl flutter diagram based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system. As my coefficients are periodic
Is my approach for determining the Stiffnesses and angles correct, and how can I change the process using the Floquet technique and the Mathieu equation? Because ultimately, the boundries are not shown correctly according to the reference.
clc;
clear all;

% Define parameters
N = 2; % Number of blades

I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b

C_thetaoverI_b = 0.00; % Damping coefficient over I_b
C_psioverI_b = 0.00; % Damping coefficient over I_b

h = 0.3; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;

gamma = 4; % lock number

V = 1000; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]

freq_1_over_Omega = 1 / Omega;

%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;

%the propeller aerodynamic coefficients
H_u = V/2;

% Frequency ranges
f_pitch= 0.001:0.5:5;
f_yaw= 0.001:0.5:5;

% Time periods for pitch and yaw
T_pitch = 1 ./ (2 * pi * f_pitch);
T_yaw = 1 ./ (2 * pi * f_yaw);

divergence_map = [];
Rdivergence_map = [];
unstable = [];

% Modify the loop to iterate over time points
for i = 1:length(T_pitch)
for j = 1:length(T_yaw)

T = max(T_pitch(i), T_yaw(j)); % Use the maximum period to cover all dynamics
t_steps = linspace(0, T, 100); % Time steps within one period

for t = t_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi / T_pitch(i);
w_omega_yaw = 2 * pi / T_yaw(j);

K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;

% Calculate matrices at time t using harmonic motion expressions
phi = 2 * pi * t / T; % Phase variation over the period

% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];

% Define damping matrix [D]
D11 = h^2*gamma*H_u*(1 – cos(2*phi)) – gamma*M_b*(1 + cos(2*phi)) – (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) – 2*(1 + cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 – cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D22 = h^2*gamma*H_u*(1 + cos(2*phi)) – gamma*M_b*(1 – cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);

D_matrix = [D11, D12;
D21, D22];

% Define stiffness matrix [K]
K11 = K_theta – h*gamma*V*H_u*(1 – cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) – gamma*V*M_u*(1 – cos(2*phi));
K22 = K_psi – h*gamma*V*H_u*(1 + cos(2*phi)) – gamma*V*M_u*sin(2*phi);

K_matrix = [K11, K12;
K21, K22];

% Compute the system matrices

M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);

P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;

P = roots([P0, P1, P2, P3, P4]);

r = 1 * P;

%Flutter
for k = 1:length(r)

if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];

% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 1e-5
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];

end
end
end

end

%Divergence

if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; t, K_psi, K_theta];
end

end
end
end

% Plotting section

figure;
hold on;

scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);

xlabel(‘K_psi’);
ylabel(‘K_theta’);

title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);

hold off; floquet technic, dynamical system MATLAB Answers — New Questions

​

Hi Everyone,
I am currently working with a very long equation. And when I output it in the livescript, MATLAB finds several sigmas (see screenshot). I would like to have this equation with the sigmas in order to process it in another function.
So far I have tried [r,sigma] = subexpr(expr). I cancelled this command after more than 1 hour of calculation. In addition, subexpr, if I have understood it correctly, only looks for one sigma.
I hope someone here can help me.

I am working with Matla 2017b and the symbolic toolboxHi Everyone,
I am currently working with a very long equation. And when I output it in the livescript, MATLAB finds several sigmas (see screenshot). I would like to have this equation with the sigmas in order to process it in another function.
So far I have tried [r,sigma] = subexpr(expr). I cancelled this command after more than 1 hour of calculation. In addition, subexpr, if I have understood it correctly, only looks for one sigma.
I hope someone here can help me.

I am working with Matla 2017b and the symbolic toolbox Hi Everyone,
I am currently working with a very long equation. And when I output it in the livescript, MATLAB finds several sigmas (see screenshot). I would like to have this equation with the sigmas in order to process it in another function.
So far I have tried [r,sigma] = subexpr(expr). I cancelled this command after more than 1 hour of calculation. In addition, subexpr, if I have understood it correctly, only looks for one sigma.
I hope someone here can help me.

I am working with Matla 2017b and the symbolic toolbox subexpr, matlab, simplify MATLAB Answers — New Questions

​

I am working on simulating heat transfer in a component that consists of multiple solid layers and an enclosed air cavity using the PDE Toolbox in MATLAB. My goal is to account for heat transfer mechanisms such as radiation and convection between the air in the cavity and the surrounding solid layers. is this possible with the pde toolbox? Is there a solutuion without modelling the air as a solid material?I am working on simulating heat transfer in a component that consists of multiple solid layers and an enclosed air cavity using the PDE Toolbox in MATLAB. My goal is to account for heat transfer mechanisms such as radiation and convection between the air in the cavity and the surrounding solid layers. is this possible with the pde toolbox? Is there a solutuion without modelling the air as a solid material? I am working on simulating heat transfer in a component that consists of multiple solid layers and an enclosed air cavity using the PDE Toolbox in MATLAB. My goal is to account for heat transfer mechanisms such as radiation and convection between the air in the cavity and the surrounding solid layers. is this possible with the pde toolbox? Is there a solutuion without modelling the air as a solid material? heat transfer, air cavity, pde MATLAB Answers — New Questions

​

Hi, I have been running a matlab script on a daily basis using windows task scheduler for the past few weeks and have noticed that the output image plots appear compressed (the arrows and contours are squished together, implying some sort of scale issue(?)) when compared with running the same script with matlab already running on my desktop. My script does specify the dimensions of the figure, but these seem to be ignored. Is there a way to correct for this? Related, if I can resolve the former issues, is there a way to force Task Scheduler to run matlab so that it is maximised (not behind the scene).
I’m using the following in Task scheduler to initiate my script:
matlab -r "cd(‘C:MatlabScripts’); run scriptname.m; pause(10)"
Thank you for any thoughts on this,
JonHi, I have been running a matlab script on a daily basis using windows task scheduler for the past few weeks and have noticed that the output image plots appear compressed (the arrows and contours are squished together, implying some sort of scale issue(?)) when compared with running the same script with matlab already running on my desktop. My script does specify the dimensions of the figure, but these seem to be ignored. Is there a way to correct for this? Related, if I can resolve the former issues, is there a way to force Task Scheduler to run matlab so that it is maximised (not behind the scene).
I’m using the following in Task scheduler to initiate my script:
matlab -r "cd(‘C:MatlabScripts’); run scriptname.m; pause(10)"
Thank you for any thoughts on this,
Jon Hi, I have been running a matlab script on a daily basis using windows task scheduler for the past few weeks and have noticed that the output image plots appear compressed (the arrows and contours are squished together, implying some sort of scale issue(?)) when compared with running the same script with matlab already running on my desktop. My script does specify the dimensions of the figure, but these seem to be ignored. Is there a way to correct for this? Related, if I can resolve the former issues, is there a way to force Task Scheduler to run matlab so that it is maximised (not behind the scene).
I’m using the following in Task scheduler to initiate my script:
matlab -r "cd(‘C:MatlabScripts’); run scriptname.m; pause(10)"
Thank you for any thoughts on this,
Jon automatically run script problem with plots MATLAB Answers — New Questions

​

How do I generate a sphere with center at [x = 4.26, y = 1.76, z = 1.62] with a radius of 0.25 units? The sphere should have 101 points. The output should look like following figure. In the figure, the blue markers are the first 51 points of the 101 points and the red markers are the remaining 50 points. Attached also find, .fig format of the required output. Below this figure, is my code. However, I am unable to reproduce the required figure. Can anyone please help?
I = imread(‘Sphere.png’);
imshow(I);
%Center of sphere
recPosition = [4.26 1.76 1.62];

rng(0,’twister’);
radius = 0.25;

%Number of points
nPoints = 100;

%Create array
test_xyz = zeros(101,3);

test_xyz(1,1) = recPosition(1);
test_xyz(1,2) = recPosition(2);
test_xyz(1,3) = recPosition(3);

radius = zeros(nPoints,1)+radius ;
azimuth = rand(size(radius)) * 2*pi;
elevation = (rand(size(radius)) .* pi)-pi/2;

[x,y,z] = sph2cart(azimuth,elevation,radius);

test_xyz(2:end,1) = test_xyz(1,1)+x;
test_xyz(2:end,2) = test_xyz(1,2)+y;
test_xyz(2:end,3) = test_xyz(1,3)+z;

figure
scatter3(test_xyz(:,1),test_xyz(:,2),test_xyz(:,3));How do I generate a sphere with center at [x = 4.26, y = 1.76, z = 1.62] with a radius of 0.25 units? The sphere should have 101 points. The output should look like following figure. In the figure, the blue markers are the first 51 points of the 101 points and the red markers are the remaining 50 points. Attached also find, .fig format of the required output. Below this figure, is my code. However, I am unable to reproduce the required figure. Can anyone please help?
I = imread(‘Sphere.png’);
imshow(I);
%Center of sphere
recPosition = [4.26 1.76 1.62];

rng(0,’twister’);
radius = 0.25;

%Number of points
nPoints = 100;

%Create array
test_xyz = zeros(101,3);

test_xyz(1,1) = recPosition(1);
test_xyz(1,2) = recPosition(2);
test_xyz(1,3) = recPosition(3);

radius = zeros(nPoints,1)+radius ;
azimuth = rand(size(radius)) * 2*pi;
elevation = (rand(size(radius)) .* pi)-pi/2;

[x,y,z] = sph2cart(azimuth,elevation,radius);

test_xyz(2:end,1) = test_xyz(1,1)+x;
test_xyz(2:end,2) = test_xyz(1,2)+y;
test_xyz(2:end,3) = test_xyz(1,3)+z;

figure
scatter3(test_xyz(:,1),test_xyz(:,2),test_xyz(:,3)); How do I generate a sphere with center at [x = 4.26, y = 1.76, z = 1.62] with a radius of 0.25 units? The sphere should have 101 points. The output should look like following figure. In the figure, the blue markers are the first 51 points of the 101 points and the red markers are the remaining 50 points. Attached also find, .fig format of the required output. Below this figure, is my code. However, I am unable to reproduce the required figure. Can anyone please help?
I = imread(‘Sphere.png’);
imshow(I);
%Center of sphere
recPosition = [4.26 1.76 1.62];

rng(0,’twister’);
radius = 0.25;

%Number of points
nPoints = 100;

%Create array
test_xyz = zeros(101,3);

test_xyz(1,1) = recPosition(1);
test_xyz(1,2) = recPosition(2);
test_xyz(1,3) = recPosition(3);

radius = zeros(nPoints,1)+radius ;
azimuth = rand(size(radius)) * 2*pi;
elevation = (rand(size(radius)) .* pi)-pi/2;

[x,y,z] = sph2cart(azimuth,elevation,radius);

test_xyz(2:end,1) = test_xyz(1,1)+x;
test_xyz(2:end,2) = test_xyz(1,2)+y;
test_xyz(2:end,3) = test_xyz(1,3)+z;

figure
scatter3(test_xyz(:,1),test_xyz(:,2),test_xyz(:,3)); sphere, sph2cart, radius MATLAB Answers — New Questions

​

Hi

I want to create a for-loop that calculates the weights of portfolios using a moving window for the period I am investigating. The moving window should move one day at a time and there are 1000 days in the window.

For example, compute the optimal weights on day1001 based on observations for period 1-1000days. Then you move on to day 1002, and re-calculate the weights, based on observations for period 2-1001days, etc.

I have a matrix of returns (Rets) that is 3740×6. The first column has the dates and the rest of the columns have daily returns for five different asset classes.

Thanks!Hi

I want to create a for-loop that calculates the weights of portfolios using a moving window for the period I am investigating. The moving window should move one day at a time and there are 1000 days in the window.

For example, compute the optimal weights on day1001 based on observations for period 1-1000days. Then you move on to day 1002, and re-calculate the weights, based on observations for period 2-1001days, etc.

I have a matrix of returns (Rets) that is 3740×6. The first column has the dates and the rest of the columns have daily returns for five different asset classes.

Thanks! Hi

I want to create a for-loop that calculates the weights of portfolios using a moving window for the period I am investigating. The moving window should move one day at a time and there are 1000 days in the window.

For example, compute the optimal weights on day1001 based on observations for period 1-1000days. Then you move on to day 1002, and re-calculate the weights, based on observations for period 2-1001days, etc.

I have a matrix of returns (Rets) that is 3740×6. The first column has the dates and the rest of the columns have daily returns for five different asset classes.

Thanks! for loop, moving window, portfolio optimization, rebalancing MATLAB Answers — New Questions

​

%Load the Aircraft Trajectory
aircraft = load("aircraftTT.mat");
aircraft1 = load("aircraftTT1.mat");
aircraft2 = load("aircraftTT2.mat");
aircraft3 = load("aircraftTT3.mat");
aircraft4 = load("aircraftTT4.mat");
geoplot(aircraft.aircraftTT.AircraftLLA(:,1), aircraft.aircraftTT.AircraftLLA(:,2), "b–",…
aircraft1.aircraftTT1.AircraftLLA1(:,1), aircraft1.aircraftTT1.AircraftLLA1(:,2), "r–",…
aircraft2.aircraftTT2.AircraftLLA2(:,1), aircraft2.aircraftTT2.AircraftLLA2(:,2), "k–",…
aircraft3.aircraftTT3.AircraftLLA3(:,1), aircraft3.aircraftTT3.AircraftLLA3(:,2), "g–",…
aircraft4.aircraftTT4.AircraftLLA4(:,1), aircraft4.aircraftTT4.AircraftLLA4(:,2), "c–"…
)
title("Aircrafts Trajectories")
%legend(‘show’);
geobasemap streets
geolimits([38.0887 38.1099],[140.6330 140.6660])
%Load the Vehicles path
gx = geoaxes;
geolimits(gx,[38.0887,38.1099],[140.6330,140.6660]);
title("Vehicles Path")
geobasemap(gx,"streets");
useManualInput = true;
if useManualInput
waypoints = zeros(15,3); %#ok<*UNRCH>
numWaypoints = size(waypoints,1);
[lat,lon] = ginput(numWaypoints);
waypoints = [lat lon ones(numWaypoints,1)];
else
waypoints = […
42.265387594440746, -71.566584651695194, 1; …
42.263702708358274, -71.569614091115028, 1; …
42.262575723421115, -71.573098554161504, 1; …
42.260479485544565, -71.575579488911472, 1; …
42.258318030927740, -71.578013745650068, 1; …
42.256140084405871, -71.580417734892066, 1; …
42.253864954613888, -71.582655636609843, 1; …
42.251556291741792, -71.584827877421873, 1; …
42.249201618438420, -71.586910025508118, 1; …
42.246834381495290, -71.588965321234468, 1; …
42.244433745609328, -71.590949233154006, 1; …
42.242264173597064, -71.593357390593781, 1; …
42.240253137371795, -71.596002172990822, 1; …
42.238444180204169, -71.598914079082988, 1; …
42.236661606027504, -71.601851940116248, 1; …
42.235163293921218, -71.605069451220544, 1; …
42.233664981814933, -71.608286962324840, 1; …
42.232221634012099, -71.611549664950445, 1; …
42.230785995268782, -71.614818705946718, 1; …
42.229491349524331, -71.618178093198480, 1; …
42.228523450114722, -71.621746854731867, 1; …
42.227591252458026, -71.625329339092573, 1; …
42.227021657836168, -71.629051198649620, 1; …
42.226488781649429, -71.632774933477108, 1; …
42.226848543969822, -71.636544256810438, 1; …
42.227208306290223, -71.640313580143768, 1; …
42.227475917842995, -71.644096605079170, 1; …
42.227730436282329, -71.647881576787640, 1; …
42.227949419806208, -71.651670397035346, 1; …
42.228142838835581, -71.655461985995458, 1; …
42.228336257864946, -71.659253574955557, 1];
lat = waypoints(:,1);
lon = waypoints(:,2);
end
geoplot(gx,lat,lon,’r-‘)
title("Vehicles Path")
% Find the distance in statute miles between each waypoint
distInDeg = distance(waypoints(1:end-1,1:2),waypoints(2:end,1:2)); % deg
distInSM = deg2sm(distInDeg); % statute miles
carSpeed = 70 / 3600; % miles per second
timeOfArrival = [0; cumsum(distInSM) / carSpeed]; % seconds
carPath = geoTrajectory(waypoints,timeOfArrival);
%Define Mission Start Date and Duration"Creating the scenario"
mission.StartDate = datetime(2024,7,24,02,42,30,TimeZone="UTC");
mission.StopTime = mission.StartDate + hours(2) + minutes(2);
mission.Duration = hours(2); %should remove, I’ve already defined above
sampleTime = 10;
mission.scenario = satelliteScenario(mission.StartDate, mission.StopTime,sampleTime);
mission.viewer = satelliteScenarioViewer(mission.scenario);
%Add Vehicles to a Satellite Scenario
car = platform(mission.scenario,carPath, Name = "Car");
%Add Aircrafts to a Satellite Scenario
aircraft = platform(mission.scenario,aircraft.aircraftTT,Name="Aircraft1", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft);
aircraft1 = platform(mission.scenario,aircraft1.aircraftTT1,Name="Aircraft2", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft1);
aircraft2 = platform(mission.scenario,aircraft2.aircraftTT2,Name="Aircraft3", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft2);
aircraft3 = platform(mission.scenario,aircraft3.aircraftTT3,Name="Aircraft4", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft3);
aircraft4 = platform(mission.scenario,aircraft4.aircraftTT4,Name="Aircraft5", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft4);
%Add a Starlink Constellation to the Satellite Scenario-TLE files
sat = satellite(mission.scenario,"2024-001B.tle", Visual3DModel="SmallSat.glb");
show(sat)
groundTrack(sat,"LeadTime",1200);
%Aircraft-to-Vehicles Access Analysis
Preq = -35.5; % Required signal power in dBm
Ptx = 17.5; % Transmitted power in dBm
% Configure the ground station, satellites, and link characteristics
% Set the ground station characteristics with parabolic telescope
%gs = struct;
car.Height = 1; % Height above the mean sea level in km
car.OpticsEfficiency = 0.8; % Optical antenna efficiency
car.ApertureDiameter = 1; % Antenna aperture diameter in m
car.PointingError = 1e-6; % Pointing error in rad
% Set the satellite A characteristics with parabolic telescope
aircraft = struct;
aircraft.Height = 550; % Height above the mean sea level in km
aircraft.OpticsEfficiency = 0.8; % Optical antenna efficiency
aircraft.ApertureDiameter = 0.07; % Antenna aperture diameter in m
aircraft.PointingError = 1e-6; % Pointing error in rad
% Set the link characteristics
link = struct;
link.Wavelength = 1550e-9; % m
link.TroposphereHeight = 20; % km (Typically ranges from 6-20 km)
link.ElevationAngle = 50; % degrees
%link.Type = "downlink"; % "downlink"|"inter-satellite"|"uplink"
% When the Type field is set to "uplink" or "downlink", you must specify
% the CloudType field, as defined in [5] table 1
link.CloudType = "Thin cirrus";
tx = transmitter(car);
rx = receiver(aircraft,MountingAngles=[0; 180; 0]);
% Calculate transmitter and receiver gain
txGain = (pi*tx.ApertureDiameter/link.Wavelength)^2;
Gtx = 10*log10(txGain); % in dB
rxGain = (pi*rx.ApertureDiameter/link.Wavelength)^2;
Grx = 10*log10(rxGain); % in dB
% Calculate transmitter and receiver pointing loss in dB
txPointingLoss = 4.3429*(txGain*(tx.PointingError)^2);
rxPointingLoss = 4.3429*(rxGain*(rx.PointingError)^2);
absorptionLoss = 0.01; % Absorption loss in dB
% Calculate the distance of the optical beam that propagates through
% the troposphere layer of the atmosphere in km
dT = (link.TroposphereHeight – gs.Height).*cscd(link.ElevationAngle);

% Calculate the slant distance for uplink and downlink between
% satellite A and the ground station for circular orbit in m
dGS = slantRangeCircularOrbit(link.ElevationAngle,satA.Height*1e3,gs.Height*1e3);

% Calculate free-space path loss between the ground station and
% satellite in dB
pathLoss = fspl(dGS,link.Wavelength);

% Calculate loss due to geometrical scattering
% cnc – cloud number concentration in cm-3
% lwc – Liquid water content in g/m-3
[cnc,lwc] = getCloudParameters(link.CloudType);
visibility = 1.002/((lwc*cnc)^0.6473); % Calculate visibility in km
% Get particle size related coefficient
if visibility<=0.5
delta = 0;
elseif visibility>0.5 && visibility<=1
delta = visibility – 0.5;
elseif visibility>1 && visibility<=6
delta = 0.16*visibility + 0.34;
elseif visibility>=6 && visibility<=50
delta = 1.3;
else % visibility>50
delta = 1.6;
end
geoCoeff = (3.91/visibility)* …
((link.Wavelength*1e9/550)^-delta); % Extinction coefficient
geoScaLoss = 4.3429*geoCoeff*dT; % Geometrical scattering loss in dB

% Calculate loss due to Mie scattering
lambda_mu = link.Wavelength*1e6; % Wavelength in microns
% Calculate empirical coefficients
a = (0.000487*(lambda_mu^3)) – (0.002237*(lambda_mu^2)) + …
(0.003864*lambda_mu) – 0.004442;
b = (-0.00573*(lambda_mu^3)) + (0.02639*(lambda_mu^2)) – …
(0.04552*lambda_mu) + 0.05164;
c = (0.02565*(lambda_mu^3)) – (0.1191*(lambda_mu^2)) + …
(0.20385*lambda_mu) – 0.216;
d = (-0.0638*(lambda_mu^3)) + (0.3034*(lambda_mu^2)) – …
(0.5083*lambda_mu) + 0.425;
mieER = a*(gs.Height^3) + b*(gs.Height^2) + …
c*(gs.Height) + d; % Extinction ratio
mieScaLoss = (4.3429*mieER)./sind(link.ElevationAngle); % Mie scattering loss in dB

% Calculate link margin for uplink or downlink in dB
linkMargin = Ptx + 10*log10(tx.OpticsEfficiency) + …
10*log10(rx.OpticsEfficiency) + Gtx + Grx – …
txPointingLoss – rxPointingLoss – pathLoss – …
absorptionLoss – geoScaLoss – mieScaLoss – Preq;
disp("Link margin for "+num2str(link.Type)+" is "+num2str(linkMargin)+" dB")%Load the Aircraft Trajectory
aircraft = load("aircraftTT.mat");
aircraft1 = load("aircraftTT1.mat");
aircraft2 = load("aircraftTT2.mat");
aircraft3 = load("aircraftTT3.mat");
aircraft4 = load("aircraftTT4.mat");
geoplot(aircraft.aircraftTT.AircraftLLA(:,1), aircraft.aircraftTT.AircraftLLA(:,2), "b–",…
aircraft1.aircraftTT1.AircraftLLA1(:,1), aircraft1.aircraftTT1.AircraftLLA1(:,2), "r–",…
aircraft2.aircraftTT2.AircraftLLA2(:,1), aircraft2.aircraftTT2.AircraftLLA2(:,2), "k–",…
aircraft3.aircraftTT3.AircraftLLA3(:,1), aircraft3.aircraftTT3.AircraftLLA3(:,2), "g–",…
aircraft4.aircraftTT4.AircraftLLA4(:,1), aircraft4.aircraftTT4.AircraftLLA4(:,2), "c–"…
)
title("Aircrafts Trajectories")
%legend(‘show’);
geobasemap streets
geolimits([38.0887 38.1099],[140.6330 140.6660])
%Load the Vehicles path
gx = geoaxes;
geolimits(gx,[38.0887,38.1099],[140.6330,140.6660]);
title("Vehicles Path")
geobasemap(gx,"streets");
useManualInput = true;
if useManualInput
waypoints = zeros(15,3); %#ok<*UNRCH>
numWaypoints = size(waypoints,1);
[lat,lon] = ginput(numWaypoints);
waypoints = [lat lon ones(numWaypoints,1)];
else
waypoints = […
42.265387594440746, -71.566584651695194, 1; …
42.263702708358274, -71.569614091115028, 1; …
42.262575723421115, -71.573098554161504, 1; …
42.260479485544565, -71.575579488911472, 1; …
42.258318030927740, -71.578013745650068, 1; …
42.256140084405871, -71.580417734892066, 1; …
42.253864954613888, -71.582655636609843, 1; …
42.251556291741792, -71.584827877421873, 1; …
42.249201618438420, -71.586910025508118, 1; …
42.246834381495290, -71.588965321234468, 1; …
42.244433745609328, -71.590949233154006, 1; …
42.242264173597064, -71.593357390593781, 1; …
42.240253137371795, -71.596002172990822, 1; …
42.238444180204169, -71.598914079082988, 1; …
42.236661606027504, -71.601851940116248, 1; …
42.235163293921218, -71.605069451220544, 1; …
42.233664981814933, -71.608286962324840, 1; …
42.232221634012099, -71.611549664950445, 1; …
42.230785995268782, -71.614818705946718, 1; …
42.229491349524331, -71.618178093198480, 1; …
42.228523450114722, -71.621746854731867, 1; …
42.227591252458026, -71.625329339092573, 1; …
42.227021657836168, -71.629051198649620, 1; …
42.226488781649429, -71.632774933477108, 1; …
42.226848543969822, -71.636544256810438, 1; …
42.227208306290223, -71.640313580143768, 1; …
42.227475917842995, -71.644096605079170, 1; …
42.227730436282329, -71.647881576787640, 1; …
42.227949419806208, -71.651670397035346, 1; …
42.228142838835581, -71.655461985995458, 1; …
42.228336257864946, -71.659253574955557, 1];
lat = waypoints(:,1);
lon = waypoints(:,2);
end
geoplot(gx,lat,lon,’r-‘)
title("Vehicles Path")
% Find the distance in statute miles between each waypoint
distInDeg = distance(waypoints(1:end-1,1:2),waypoints(2:end,1:2)); % deg
distInSM = deg2sm(distInDeg); % statute miles
carSpeed = 70 / 3600; % miles per second
timeOfArrival = [0; cumsum(distInSM) / carSpeed]; % seconds
carPath = geoTrajectory(waypoints,timeOfArrival);
%Define Mission Start Date and Duration"Creating the scenario"
mission.StartDate = datetime(2024,7,24,02,42,30,TimeZone="UTC");
mission.StopTime = mission.StartDate + hours(2) + minutes(2);
mission.Duration = hours(2); %should remove, I’ve already defined above
sampleTime = 10;
mission.scenario = satelliteScenario(mission.StartDate, mission.StopTime,sampleTime);
mission.viewer = satelliteScenarioViewer(mission.scenario);
%Add Vehicles to a Satellite Scenario
car = platform(mission.scenario,carPath, Name = "Car");
%Add Aircrafts to a Satellite Scenario
aircraft = platform(mission.scenario,aircraft.aircraftTT,Name="Aircraft1", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft);
aircraft1 = platform(mission.scenario,aircraft1.aircraftTT1,Name="Aircraft2", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft1);
aircraft2 = platform(mission.scenario,aircraft2.aircraftTT2,Name="Aircraft3", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft2);
aircraft3 = platform(mission.scenario,aircraft3.aircraftTT3,Name="Aircraft4", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft3);
aircraft4 = platform(mission.scenario,aircraft4.aircraftTT4,Name="Aircraft5", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft4);
%Add a Starlink Constellation to the Satellite Scenario-TLE files
sat = satellite(mission.scenario,"2024-001B.tle", Visual3DModel="SmallSat.glb");
show(sat)
groundTrack(sat,"LeadTime",1200);
%Aircraft-to-Vehicles Access Analysis
Preq = -35.5; % Required signal power in dBm
Ptx = 17.5; % Transmitted power in dBm
% Configure the ground station, satellites, and link characteristics
% Set the ground station characteristics with parabolic telescope
%gs = struct;
car.Height = 1; % Height above the mean sea level in km
car.OpticsEfficiency = 0.8; % Optical antenna efficiency
car.ApertureDiameter = 1; % Antenna aperture diameter in m
car.PointingError = 1e-6; % Pointing error in rad
% Set the satellite A characteristics with parabolic telescope
aircraft = struct;
aircraft.Height = 550; % Height above the mean sea level in km
aircraft.OpticsEfficiency = 0.8; % Optical antenna efficiency
aircraft.ApertureDiameter = 0.07; % Antenna aperture diameter in m
aircraft.PointingError = 1e-6; % Pointing error in rad
% Set the link characteristics
link = struct;
link.Wavelength = 1550e-9; % m
link.TroposphereHeight = 20; % km (Typically ranges from 6-20 km)
link.ElevationAngle = 50; % degrees
%link.Type = "downlink"; % "downlink"|"inter-satellite"|"uplink"
% When the Type field is set to "uplink" or "downlink", you must specify
% the CloudType field, as defined in [5] table 1
link.CloudType = "Thin cirrus";
tx = transmitter(car);
rx = receiver(aircraft,MountingAngles=[0; 180; 0]);
% Calculate transmitter and receiver gain
txGain = (pi*tx.ApertureDiameter/link.Wavelength)^2;
Gtx = 10*log10(txGain); % in dB
rxGain = (pi*rx.ApertureDiameter/link.Wavelength)^2;
Grx = 10*log10(rxGain); % in dB
% Calculate transmitter and receiver pointing loss in dB
txPointingLoss = 4.3429*(txGain*(tx.PointingError)^2);
rxPointingLoss = 4.3429*(rxGain*(rx.PointingError)^2);
absorptionLoss = 0.01; % Absorption loss in dB
% Calculate the distance of the optical beam that propagates through
% the troposphere layer of the atmosphere in km
dT = (link.TroposphereHeight – gs.Height).*cscd(link.ElevationAngle);

% Calculate the slant distance for uplink and downlink between
% satellite A and the ground station for circular orbit in m
dGS = slantRangeCircularOrbit(link.ElevationAngle,satA.Height*1e3,gs.Height*1e3);

% Calculate free-space path loss between the ground station and
% satellite in dB
pathLoss = fspl(dGS,link.Wavelength);

% Calculate loss due to geometrical scattering
% cnc – cloud number concentration in cm-3
% lwc – Liquid water content in g/m-3
[cnc,lwc] = getCloudParameters(link.CloudType);
visibility = 1.002/((lwc*cnc)^0.6473); % Calculate visibility in km
% Get particle size related coefficient
if visibility<=0.5
delta = 0;
elseif visibility>0.5 && visibility<=1
delta = visibility – 0.5;
elseif visibility>1 && visibility<=6
delta = 0.16*visibility + 0.34;
elseif visibility>=6 && visibility<=50
delta = 1.3;
else % visibility>50
delta = 1.6;
end
geoCoeff = (3.91/visibility)* …
((link.Wavelength*1e9/550)^-delta); % Extinction coefficient
geoScaLoss = 4.3429*geoCoeff*dT; % Geometrical scattering loss in dB

% Calculate loss due to Mie scattering
lambda_mu = link.Wavelength*1e6; % Wavelength in microns
% Calculate empirical coefficients
a = (0.000487*(lambda_mu^3)) – (0.002237*(lambda_mu^2)) + …
(0.003864*lambda_mu) – 0.004442;
b = (-0.00573*(lambda_mu^3)) + (0.02639*(lambda_mu^2)) – …
(0.04552*lambda_mu) + 0.05164;
c = (0.02565*(lambda_mu^3)) – (0.1191*(lambda_mu^2)) + …
(0.20385*lambda_mu) – 0.216;
d = (-0.0638*(lambda_mu^3)) + (0.3034*(lambda_mu^2)) – …
(0.5083*lambda_mu) + 0.425;
mieER = a*(gs.Height^3) + b*(gs.Height^2) + …
c*(gs.Height) + d; % Extinction ratio
mieScaLoss = (4.3429*mieER)./sind(link.ElevationAngle); % Mie scattering loss in dB

% Calculate link margin for uplink or downlink in dB
linkMargin = Ptx + 10*log10(tx.OpticsEfficiency) + …
10*log10(rx.OpticsEfficiency) + Gtx + Grx – …
txPointingLoss – rxPointingLoss – pathLoss – …
absorptionLoss – geoScaLoss – mieScaLoss – Preq;
disp("Link margin for "+num2str(link.Type)+" is "+num2str(linkMargin)+" dB") %Load the Aircraft Trajectory
aircraft = load("aircraftTT.mat");
aircraft1 = load("aircraftTT1.mat");
aircraft2 = load("aircraftTT2.mat");
aircraft3 = load("aircraftTT3.mat");
aircraft4 = load("aircraftTT4.mat");
geoplot(aircraft.aircraftTT.AircraftLLA(:,1), aircraft.aircraftTT.AircraftLLA(:,2), "b–",…
aircraft1.aircraftTT1.AircraftLLA1(:,1), aircraft1.aircraftTT1.AircraftLLA1(:,2), "r–",…
aircraft2.aircraftTT2.AircraftLLA2(:,1), aircraft2.aircraftTT2.AircraftLLA2(:,2), "k–",…
aircraft3.aircraftTT3.AircraftLLA3(:,1), aircraft3.aircraftTT3.AircraftLLA3(:,2), "g–",…
aircraft4.aircraftTT4.AircraftLLA4(:,1), aircraft4.aircraftTT4.AircraftLLA4(:,2), "c–"…
)
title("Aircrafts Trajectories")
%legend(‘show’);
geobasemap streets
geolimits([38.0887 38.1099],[140.6330 140.6660])
%Load the Vehicles path
gx = geoaxes;
geolimits(gx,[38.0887,38.1099],[140.6330,140.6660]);
title("Vehicles Path")
geobasemap(gx,"streets");
useManualInput = true;
if useManualInput
waypoints = zeros(15,3); %#ok<*UNRCH>
numWaypoints = size(waypoints,1);
[lat,lon] = ginput(numWaypoints);
waypoints = [lat lon ones(numWaypoints,1)];
else
waypoints = […
42.265387594440746, -71.566584651695194, 1; …
42.263702708358274, -71.569614091115028, 1; …
42.262575723421115, -71.573098554161504, 1; …
42.260479485544565, -71.575579488911472, 1; …
42.258318030927740, -71.578013745650068, 1; …
42.256140084405871, -71.580417734892066, 1; …
42.253864954613888, -71.582655636609843, 1; …
42.251556291741792, -71.584827877421873, 1; …
42.249201618438420, -71.586910025508118, 1; …
42.246834381495290, -71.588965321234468, 1; …
42.244433745609328, -71.590949233154006, 1; …
42.242264173597064, -71.593357390593781, 1; …
42.240253137371795, -71.596002172990822, 1; …
42.238444180204169, -71.598914079082988, 1; …
42.236661606027504, -71.601851940116248, 1; …
42.235163293921218, -71.605069451220544, 1; …
42.233664981814933, -71.608286962324840, 1; …
42.232221634012099, -71.611549664950445, 1; …
42.230785995268782, -71.614818705946718, 1; …
42.229491349524331, -71.618178093198480, 1; …
42.228523450114722, -71.621746854731867, 1; …
42.227591252458026, -71.625329339092573, 1; …
42.227021657836168, -71.629051198649620, 1; …
42.226488781649429, -71.632774933477108, 1; …
42.226848543969822, -71.636544256810438, 1; …
42.227208306290223, -71.640313580143768, 1; …
42.227475917842995, -71.644096605079170, 1; …
42.227730436282329, -71.647881576787640, 1; …
42.227949419806208, -71.651670397035346, 1; …
42.228142838835581, -71.655461985995458, 1; …
42.228336257864946, -71.659253574955557, 1];
lat = waypoints(:,1);
lon = waypoints(:,2);
end
geoplot(gx,lat,lon,’r-‘)
title("Vehicles Path")
% Find the distance in statute miles between each waypoint
distInDeg = distance(waypoints(1:end-1,1:2),waypoints(2:end,1:2)); % deg
distInSM = deg2sm(distInDeg); % statute miles
carSpeed = 70 / 3600; % miles per second
timeOfArrival = [0; cumsum(distInSM) / carSpeed]; % seconds
carPath = geoTrajectory(waypoints,timeOfArrival);
%Define Mission Start Date and Duration"Creating the scenario"
mission.StartDate = datetime(2024,7,24,02,42,30,TimeZone="UTC");
mission.StopTime = mission.StartDate + hours(2) + minutes(2);
mission.Duration = hours(2); %should remove, I’ve already defined above
sampleTime = 10;
mission.scenario = satelliteScenario(mission.StartDate, mission.StopTime,sampleTime);
mission.viewer = satelliteScenarioViewer(mission.scenario);
%Add Vehicles to a Satellite Scenario
car = platform(mission.scenario,carPath, Name = "Car");
%Add Aircrafts to a Satellite Scenario
aircraft = platform(mission.scenario,aircraft.aircraftTT,Name="Aircraft1", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft);
aircraft1 = platform(mission.scenario,aircraft1.aircraftTT1,Name="Aircraft2", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft1);
aircraft2 = platform(mission.scenario,aircraft2.aircraftTT2,Name="Aircraft3", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft2);
aircraft3 = platform(mission.scenario,aircraft3.aircraftTT3,Name="Aircraft4", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft3);
aircraft4 = platform(mission.scenario,aircraft4.aircraftTT4,Name="Aircraft5", Visual3DModel="NarrowBodyAirliner.glb")
%camtarget(mission.viewer,aircraft4);
%Add a Starlink Constellation to the Satellite Scenario-TLE files
sat = satellite(mission.scenario,"2024-001B.tle", Visual3DModel="SmallSat.glb");
show(sat)
groundTrack(sat,"LeadTime",1200);
%Aircraft-to-Vehicles Access Analysis
Preq = -35.5; % Required signal power in dBm
Ptx = 17.5; % Transmitted power in dBm
% Configure the ground station, satellites, and link characteristics
% Set the ground station characteristics with parabolic telescope
%gs = struct;
car.Height = 1; % Height above the mean sea level in km
car.OpticsEfficiency = 0.8; % Optical antenna efficiency
car.ApertureDiameter = 1; % Antenna aperture diameter in m
car.PointingError = 1e-6; % Pointing error in rad
% Set the satellite A characteristics with parabolic telescope
aircraft = struct;
aircraft.Height = 550; % Height above the mean sea level in km
aircraft.OpticsEfficiency = 0.8; % Optical antenna efficiency
aircraft.ApertureDiameter = 0.07; % Antenna aperture diameter in m
aircraft.PointingError = 1e-6; % Pointing error in rad
% Set the link characteristics
link = struct;
link.Wavelength = 1550e-9; % m
link.TroposphereHeight = 20; % km (Typically ranges from 6-20 km)
link.ElevationAngle = 50; % degrees
%link.Type = "downlink"; % "downlink"|"inter-satellite"|"uplink"
% When the Type field is set to "uplink" or "downlink", you must specify
% the CloudType field, as defined in [5] table 1
link.CloudType = "Thin cirrus";
tx = transmitter(car);
rx = receiver(aircraft,MountingAngles=[0; 180; 0]);
% Calculate transmitter and receiver gain
txGain = (pi*tx.ApertureDiameter/link.Wavelength)^2;
Gtx = 10*log10(txGain); % in dB
rxGain = (pi*rx.ApertureDiameter/link.Wavelength)^2;
Grx = 10*log10(rxGain); % in dB
% Calculate transmitter and receiver pointing loss in dB
txPointingLoss = 4.3429*(txGain*(tx.PointingError)^2);
rxPointingLoss = 4.3429*(rxGain*(rx.PointingError)^2);
absorptionLoss = 0.01; % Absorption loss in dB
% Calculate the distance of the optical beam that propagates through
% the troposphere layer of the atmosphere in km
dT = (link.TroposphereHeight – gs.Height).*cscd(link.ElevationAngle);

% Calculate the slant distance for uplink and downlink between
% satellite A and the ground station for circular orbit in m
dGS = slantRangeCircularOrbit(link.ElevationAngle,satA.Height*1e3,gs.Height*1e3);

% Calculate free-space path loss between the ground station and
% satellite in dB
pathLoss = fspl(dGS,link.Wavelength);

% Calculate loss due to geometrical scattering
% cnc – cloud number concentration in cm-3
% lwc – Liquid water content in g/m-3
[cnc,lwc] = getCloudParameters(link.CloudType);
visibility = 1.002/((lwc*cnc)^0.6473); % Calculate visibility in km
% Get particle size related coefficient
if visibility<=0.5
delta = 0;
elseif visibility>0.5 && visibility<=1
delta = visibility – 0.5;
elseif visibility>1 && visibility<=6
delta = 0.16*visibility + 0.34;
elseif visibility>=6 && visibility<=50
delta = 1.3;
else % visibility>50
delta = 1.6;
end
geoCoeff = (3.91/visibility)* …
((link.Wavelength*1e9/550)^-delta); % Extinction coefficient
geoScaLoss = 4.3429*geoCoeff*dT; % Geometrical scattering loss in dB

% Calculate loss due to Mie scattering
lambda_mu = link.Wavelength*1e6; % Wavelength in microns
% Calculate empirical coefficients
a = (0.000487*(lambda_mu^3)) – (0.002237*(lambda_mu^2)) + …
(0.003864*lambda_mu) – 0.004442;
b = (-0.00573*(lambda_mu^3)) + (0.02639*(lambda_mu^2)) – …
(0.04552*lambda_mu) + 0.05164;
c = (0.02565*(lambda_mu^3)) – (0.1191*(lambda_mu^2)) + …
(0.20385*lambda_mu) – 0.216;
d = (-0.0638*(lambda_mu^3)) + (0.3034*(lambda_mu^2)) – …
(0.5083*lambda_mu) + 0.425;
mieER = a*(gs.Height^3) + b*(gs.Height^2) + …
c*(gs.Height) + d; % Extinction ratio
mieScaLoss = (4.3429*mieER)./sind(link.ElevationAngle); % Mie scattering loss in dB

% Calculate link margin for uplink or downlink in dB
linkMargin = Ptx + 10*log10(tx.OpticsEfficiency) + …
10*log10(rx.OpticsEfficiency) + Gtx + Grx – …
txPointingLoss – rxPointingLoss – pathLoss – …
absorptionLoss – geoScaLoss – mieScaLoss – Preq;
disp("Link margin for "+num2str(link.Type)+" is "+num2str(linkMargin)+" dB") satellite, car MATLAB Answers — New Questions

​