Post Content Post Content license MATLAB Answers — New Questions

​

Problem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.

the code: –
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)

% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)

% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)

% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)

% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)

% Distance between center of mass and center of pressure
d_cp_cm = cp_z – cm_z;

% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)

% Calculate initial velocity components
v0x = v0_aircraft – v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);

% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];

% Time settings
tspan = [0 300]; % Time span for integration

% ODE solver options
options = odeset(‘Events’, @ground_impact, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);

% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);

% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);

% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);

% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel(‘Horizontal distance (km)’);
ylabel(‘Altitude (km)’);
title(‘2D Trajectory of Body (Advanced 6DOF Simulation)’);
grid on;

% Plot velocity
figure;
plot(t, v);
xlabel(‘Time (s)’);
ylabel(‘Velocity (m/s)’);
title(‘Velocity vs Time’);
grid on;

% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel(‘Angular Position (degrees)’);
title(‘Angular Motion’);
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel(‘Time (s)’);
ylabel(‘Angular Velocity (deg/s)’);
grid on;

% Display results
disp([‘Time of flight: ‘, num2str(t(end)), ‘ s’]);
disp([‘Horizontal distance traveled: ‘, num2str(x(end)/1000), ‘ km’]);
disp([‘Impact velocity: ‘, num2str(v(end)), ‘ m/s’]);
disp([‘Final angle: ‘, num2str(rad2deg(theta(end))), ‘ degrees’]);
disp([‘Max altitude: ‘, num2str(max(h)/1000), ‘ km’]);
disp([‘Max velocity: ‘, num2str(max(v)), ‘ m/s’]);

% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp([‘State at t = ‘, num2str(t(idx)), ‘ s:’]);
disp([‘ Position: (‘, num2str(x(idx)/1000), ‘ km, ‘, num2str(h(idx)/1000), ‘ km)’]);
disp([‘ Velocity: ‘, num2str(v(idx)), ‘ m/s’]);
disp([‘ Angle: ‘, num2str(rad2deg(theta(idx))), ‘ degrees’]);
end
end
end

function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);

% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) – theta;

% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);

% Calculate Mach number
M = v / a;

% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);

% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment

% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;

% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end

% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) – Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration

% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end

function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air

% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)

if h < 11000 % Troposphere
T = T0 – 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h – 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h – 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end

rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end

function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models

% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);

% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;

% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 – M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility

% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;

% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end

function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
endProblem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.

the code: –
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)

% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)

% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)

% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)

% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)

% Distance between center of mass and center of pressure
d_cp_cm = cp_z – cm_z;

% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)

% Calculate initial velocity components
v0x = v0_aircraft – v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);

% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];

% Time settings
tspan = [0 300]; % Time span for integration

% ODE solver options
options = odeset(‘Events’, @ground_impact, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);

% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);

% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);

% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);

% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel(‘Horizontal distance (km)’);
ylabel(‘Altitude (km)’);
title(‘2D Trajectory of Body (Advanced 6DOF Simulation)’);
grid on;

% Plot velocity
figure;
plot(t, v);
xlabel(‘Time (s)’);
ylabel(‘Velocity (m/s)’);
title(‘Velocity vs Time’);
grid on;

% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel(‘Angular Position (degrees)’);
title(‘Angular Motion’);
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel(‘Time (s)’);
ylabel(‘Angular Velocity (deg/s)’);
grid on;

% Display results
disp([‘Time of flight: ‘, num2str(t(end)), ‘ s’]);
disp([‘Horizontal distance traveled: ‘, num2str(x(end)/1000), ‘ km’]);
disp([‘Impact velocity: ‘, num2str(v(end)), ‘ m/s’]);
disp([‘Final angle: ‘, num2str(rad2deg(theta(end))), ‘ degrees’]);
disp([‘Max altitude: ‘, num2str(max(h)/1000), ‘ km’]);
disp([‘Max velocity: ‘, num2str(max(v)), ‘ m/s’]);

% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp([‘State at t = ‘, num2str(t(idx)), ‘ s:’]);
disp([‘ Position: (‘, num2str(x(idx)/1000), ‘ km, ‘, num2str(h(idx)/1000), ‘ km)’]);
disp([‘ Velocity: ‘, num2str(v(idx)), ‘ m/s’]);
disp([‘ Angle: ‘, num2str(rad2deg(theta(idx))), ‘ degrees’]);
end
end
end

function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);

% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) – theta;

% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);

% Calculate Mach number
M = v / a;

% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);

% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment

% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;

% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end

% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) – Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration

% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end

function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air

% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)

if h < 11000 % Troposphere
T = T0 – 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h – 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h – 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end

rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end

function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models

% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);

% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;

% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 – M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility

% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;

% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end

function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
end Problem Statement: –
I have a 3d body which i have launched from cruising aircraft travelling at 200m/sec at 10km altitude. The launch speed is 50m/sec and it is launched at 45degrees angle pointing downwards. I want to create the whole trajectory of the body using everything like change in pressure, density, drag, lift etc. I created a code but it is not giving the accurate answer.

the code: –
function aircraft_launch_trajectory_6DOF_advanced()
% Constants
g = 9.81; % Acceleration due to gravity (m/s^2)

% Convert given units to SI
m = 49.964 * 0.0283495; % Mass (kg) (1 ouncemass = 0.0283495 kg)
A = 68.166 * 0.00064516; % Area (m^2) (1 in^2 = 0.00064516 m^2)
Izz = 19.58 * 0.0283495 * 0.0254^2; % Moment of inertia (kg*m^2)

% Body dimensions
length = 1.80 * 0.0254; % Length (m)
width = 3.30 * 0.0254; % Width (m)
height = 8.00 * 0.0254; % Height (m)

% Center of mass location (assuming from the base of the body)
cm_z = 1.708 * 0.0254; % (m)

% Assumed center of pressure location (you may need to adjust this)
cp_z = 4.5 * 0.0254; % (m)

% Distance between center of mass and center of pressure
d_cp_cm = cp_z – cm_z;

% Initial conditions
h0 = 10000; % Initial height (m)
v0_aircraft = 200; % Aircraft velocity (m/s)
v0_launch = 50; % Launch velocity (m/s)
angle = -45; % Launch angle (degrees)
omega0 = 0; % Initial angular velocity (rad/s)

% Calculate initial velocity components
v0x = v0_aircraft – v0_launch * cosd(angle);
v0y = -v0_launch * sind(angle);

% Initial state vector [x, y, vx, vy, theta, omega]
y0 = [0; h0; v0x; v0y; deg2rad(angle); omega0];

% Time settings
tspan = [0 300]; % Time span for integration

% ODE solver options
options = odeset(‘Events’, @ground_impact, ‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);

% Solve ODE using ode45 (4th order Runge-Kutta)
[t, y] = ode45(@(t, y) equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length), tspan, y0, options);

% Extract results
x = y(:, 1);
h = y(:, 2);
vx = y(:, 3);
vy = y(:, 4);
theta = y(:, 5);
omega = y(:, 6);

% Calculate velocity magnitude
v = sqrt(vx.^2 + vy.^2);

% Plot trajectory
figure;
plot(x/1000, h/1000);
xlabel(‘Horizontal distance (km)’);
ylabel(‘Altitude (km)’);
title(‘2D Trajectory of Body (Advanced 6DOF Simulation)’);
grid on;

% Plot velocity
figure;
plot(t, v);
xlabel(‘Time (s)’);
ylabel(‘Velocity (m/s)’);
title(‘Velocity vs Time’);
grid on;

% Plot angular position and velocity
figure;
subplot(2,1,1);
plot(t, rad2deg(theta));
ylabel(‘Angular Position (degrees)’);
title(‘Angular Motion’);
grid on;
subplot(2,1,2);
plot(t, rad2deg(omega));
xlabel(‘Time (s)’);
ylabel(‘Angular Velocity (deg/s)’);
grid on;

% Display results
disp([‘Time of flight: ‘, num2str(t(end)), ‘ s’]);
disp([‘Horizontal distance traveled: ‘, num2str(x(end)/1000), ‘ km’]);
disp([‘Impact velocity: ‘, num2str(v(end)), ‘ m/s’]);
disp([‘Final angle: ‘, num2str(rad2deg(theta(end))), ‘ degrees’]);
disp([‘Max altitude: ‘, num2str(max(h)/1000), ‘ km’]);
disp([‘Max velocity: ‘, num2str(max(v)), ‘ m/s’]);

% Output state at specific times for debugging
debug_times = [0, 5, 10, 30, 60, 120];
for debug_time = debug_times
idx = find(t >= debug_time, 1);
if ~isempty(idx)
disp([‘State at t = ‘, num2str(t(idx)), ‘ s:’]);
disp([‘ Position: (‘, num2str(x(idx)/1000), ‘ km, ‘, num2str(h(idx)/1000), ‘ km)’]);
disp([‘ Velocity: ‘, num2str(v(idx)), ‘ m/s’]);
disp([‘ Angle: ‘, num2str(rad2deg(theta(idx))), ‘ degrees’]);
end
end
end

function dydt = equations_of_motion_6DOF_advanced(t, y, m, Izz, A, d_cp_cm, length)
% Unpack state vector
x = y(1);
h = y(2);
vx = y(3);
vy = y(4);
theta = y(5);
omega = y(6);

% Calculate velocity magnitude and angle of attack
v = sqrt(vx^2 + vy^2);
alpha = atan2(vy, vx) – theta;

% Calculate air density and speed of sound using improved atmospheric model
[rho, a] = atmosphere_model(h);

% Calculate Mach number
M = v / a;

% Get aerodynamic coefficients based on angle of attack and Mach number
[Cd, Cl, Cm] = aero_coefficients(alpha, M);

% Calculate aerodynamic forces and moment
q = 0.5 * rho * v^2; % Dynamic pressure
Fd = q * Cd * A; % Drag force
Fl = q * Cl * A; % Lift force
M_aero = q * Cm * A * length; % Aerodynamic moment

% Add moment due to center of pressure offset
M_total = M_aero + Fl * d_cp_cm;

% Apply initial boost for first 7 seconds
boost_duration = 7; % seconds
if t <= boost_duration
boost_force = m * 200 / boost_duration; % Force to achieve 200 m/s in 7 seconds
boost_fx = boost_force * cos(theta);
boost_fy = boost_force * sin(theta);
else
boost_fx = 0;
boost_fy = 0;
end

% Calculate accelerations
ax = (-Fd * cos(theta) + Fl * sin(theta) + boost_fx) / m;
ay = -9.81 + (-Fd * sin(theta) – Fl * cos(theta) + boost_fy) / m;
alpha_dot = M_total / Izz; % Angular acceleration

% Return state derivatives
dydt = [vx; vy; ax; ay; omega; alpha_dot];
end

function [rho, a] = atmosphere_model(h)
% Constants
R = 287.05; % Specific gas constant for air (J/(kg*K))
g0 = 9.80665; % Standard gravity (m/s^2)
gamma = 1.4; % Ratio of specific heats for air

% Sea level conditions
T0 = 288.15; % Temperature at sea level (K)
P0 = 101325; % Pressure at sea level (Pa)

if h < 11000 % Troposphere
T = T0 – 0.0065 * h;
P = P0 * (T / T0)^(g0 / (R * 0.0065));
elseif h < 20000 % Lower Stratosphere
T = 216.65;
P = 22632.1 * exp(-g0 * (h – 11000) / (R * T));
else % Upper Stratosphere
T = 216.65 + 0.001 * (h – 20000);
P = 5474.89 * (216.65 / T)^(34.1632);
end

rho = P / (R * T);
a = sqrt(gamma * R * T); % Speed of sound
end

function [Cd, Cl, Cm] = aero_coefficients(alpha, M)
% This function returns aerodynamic coefficients based on angle of attack and Mach number
% Note: These are simplified models and should be replaced with actual data or more complex models

% Convert alpha to degrees for easier calculations
alpha_deg = rad2deg(alpha);

% Drag coefficient
Cd0 = 0.1 + 0.1 * M; % Base drag increases with Mach number
Cdi = 0.01 * alpha_deg^2; % Induced drag
Cd = Cd0 + Cdi;

% Lift coefficient
Cl_slope = 0.1; % Lift curve slope
Cl = Cl_slope * alpha_deg * (1 – M^2/4)^(-0.5); % Prandtl-Glauert correction for compressibility

% Moment coefficient
Cm0 = -0.05; % Zero-lift moment coefficient
Cm_alpha = -0.01; % Moment curve slope
Cm = Cm0 + Cm_alpha * alpha_deg;

% Limit coefficients to realistic ranges
Cd = max(0, min(Cd, 2));
Cl = max(-2, min(Cl, 2));
Cm = max(-1, min(Cm, 1));
end

function [position, isterminal, direction] = ground_impact(t, y)
position = y(2); % Height
isterminal = 1; % Stop integration when ground is reached
direction = -1; % Negative direction (approaching ground)
end trajectory, aerospace MATLAB Answers — New Questions

​

% I want to do feature extraction for signals and then apply ANN classification on them,, can you help please?

Fs = 50;
sf = waveletScattering(‘SignalLength’, 4950, ‘SamplingFrequency’, Fs) ;

s = cell(120,1);
parfor j = 1:120;
r = featureMatrix(sf, f(j, :)) ;
s{j, 1} = r;
end

h = cell2mat(s) ;
b = cell2table(s) ;
c = table2dataset(b);

%the single signal was in one row, but after this it becomes 202rows×10columns, and I can’t apply ANN like this
Any ideas?% I want to do feature extraction for signals and then apply ANN classification on them,, can you help please?

Fs = 50;
sf = waveletScattering(‘SignalLength’, 4950, ‘SamplingFrequency’, Fs) ;

s = cell(120,1);
parfor j = 1:120;
r = featureMatrix(sf, f(j, :)) ;
s{j, 1} = r;
end

h = cell2mat(s) ;
b = cell2table(s) ;
c = table2dataset(b);

%the single signal was in one row, but after this it becomes 202rows×10columns, and I can’t apply ANN like this
Any ideas? % I want to do feature extraction for signals and then apply ANN classification on them,, can you help please?

Fs = 50;
sf = waveletScattering(‘SignalLength’, 4950, ‘SamplingFrequency’, Fs) ;

s = cell(120,1);
parfor j = 1:120;
r = featureMatrix(sf, f(j, :)) ;
s{j, 1} = r;
end

h = cell2mat(s) ;
b = cell2table(s) ;
c = table2dataset(b);

%the single signal was in one row, but after this it becomes 202rows×10columns, and I can’t apply ANN like this
Any ideas? feature extraction, neural network, classification, wavelet, wavelet scattering, ecg signals, ecg recognition, ecg as signature, deep learning, machine learning MATLAB Answers — New Questions

​

Hi,

I want to find the sensitivity of model output which is a tissue mass for models parameters. I received an red flag that I can’t do that for repeaded variable.
How can I handle that?

Thank you.Hi,

I want to find the sensitivity of model output which is a tissue mass for models parameters. I received an red flag that I can’t do that for repeaded variable.
How can I handle that?

Thank you. Hi,

I want to find the sensitivity of model output which is a tissue mass for models parameters. I received an red flag that I can’t do that for repeaded variable.
How can I handle that?

Thank you. sensitivity MATLAB Answers — New Questions

​

How can I set the format of the ticks of the y axis?
I want to change the upper one to the lower one, as the former is too wide in space.How can I set the format of the ticks of the y axis?
I want to change the upper one to the lower one, as the former is too wide in space. How can I set the format of the ticks of the y axis?
I want to change the upper one to the lower one, as the former is too wide in space. plotting MATLAB Answers — New Questions

​

I’m trying to analyze the frequency accuracy of the DFT as a function of DFT size assuming a rectangular
window. Using single sinusoidal waveform sampled at 5 kHz vary its frequency in small
increments to show the how the frequency accuracy changes for at least 4 different DFT length,
256, 512, 1024, and 4096.

Here is my code:

%Part 1
% MATLAB code to analyze the frequency accuracy of the DFT manually
fs = 5000; % Sampling frequency 5 kHz
f_start = 1000; % Start frequency of sinusoid (1 kHz)
f_end = 1050; % End frequency of sinusoid (1.05 kHz)
f_step = 1; % Frequency step size (small increments)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Prepare figure
figure;
hold on;
for N = N_dft
freq_accuracy = []; % Store accuracy results for each DFT size
for f = f_start:f_step:f_end
% Generate sinusoidal waveform
x = sin(2*pi*f*t);

rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;

% Manually compute the DFT
X = x_windowed*dftmtx(N);

% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);

% Find the index of the peak in the magnitude of the DFT
[~, peak_idx] = max(abs(X));

% Estimated frequency from the DFT
f_estimated = freq_axis(peak_idx);

% Store the frequency error (absolute difference between true and estimated frequency)
freq_error = abs(f_estimated – f);
freq_accuracy = [freq_accuracy, freq_error];
end

% Plot the frequency accuracy for the current DFT size
plot(f_start:f_step:f_end, freq_accuracy, ‘DisplayName’, [‘N = ‘, num2str(N)]);
end
% Add labels and legend
xlabel(‘True Frequency (Hz)’);
ylabel(‘Frequency Error (Hz)’);
title(‘Frequency Accuracy of the DFT for Different DFT Sizes (Manual DFT)’);
legend(‘show’);
grid on;
hold off;

However, the result I got is really weird.

At some frequency the estimate error between true frequency can be 3000, while most of the time it’s under 1.

However, it I substitue this line: X = x_windowed*dftmtx(N);
with this line: X = fft(x_windowed, N);

Then I got this image

It seems that this is the correct figure.
I read from this website that fft and dftmtx function can be interchangeble
https://www.mathworks.com/help/signal/ref/dftmtx.html
But my results doesn’t approve it. I want to know why.

Besides, am I applying windowing function in a correct way?
The code is as following:
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;I’m trying to analyze the frequency accuracy of the DFT as a function of DFT size assuming a rectangular
window. Using single sinusoidal waveform sampled at 5 kHz vary its frequency in small
increments to show the how the frequency accuracy changes for at least 4 different DFT length,
256, 512, 1024, and 4096.

Here is my code:

%Part 1
% MATLAB code to analyze the frequency accuracy of the DFT manually
fs = 5000; % Sampling frequency 5 kHz
f_start = 1000; % Start frequency of sinusoid (1 kHz)
f_end = 1050; % End frequency of sinusoid (1.05 kHz)
f_step = 1; % Frequency step size (small increments)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Prepare figure
figure;
hold on;
for N = N_dft
freq_accuracy = []; % Store accuracy results for each DFT size
for f = f_start:f_step:f_end
% Generate sinusoidal waveform
x = sin(2*pi*f*t);

rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;

% Manually compute the DFT
X = x_windowed*dftmtx(N);

% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);

% Find the index of the peak in the magnitude of the DFT
[~, peak_idx] = max(abs(X));

% Estimated frequency from the DFT
f_estimated = freq_axis(peak_idx);

% Store the frequency error (absolute difference between true and estimated frequency)
freq_error = abs(f_estimated – f);
freq_accuracy = [freq_accuracy, freq_error];
end

% Plot the frequency accuracy for the current DFT size
plot(f_start:f_step:f_end, freq_accuracy, ‘DisplayName’, [‘N = ‘, num2str(N)]);
end
% Add labels and legend
xlabel(‘True Frequency (Hz)’);
ylabel(‘Frequency Error (Hz)’);
title(‘Frequency Accuracy of the DFT for Different DFT Sizes (Manual DFT)’);
legend(‘show’);
grid on;
hold off;

However, the result I got is really weird.

At some frequency the estimate error between true frequency can be 3000, while most of the time it’s under 1.

However, it I substitue this line: X = x_windowed*dftmtx(N);
with this line: X = fft(x_windowed, N);

Then I got this image

It seems that this is the correct figure.
I read from this website that fft and dftmtx function can be interchangeble
https://www.mathworks.com/help/signal/ref/dftmtx.html
But my results doesn’t approve it. I want to know why.

Besides, am I applying windowing function in a correct way?
The code is as following:
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window; I’m trying to analyze the frequency accuracy of the DFT as a function of DFT size assuming a rectangular
window. Using single sinusoidal waveform sampled at 5 kHz vary its frequency in small
increments to show the how the frequency accuracy changes for at least 4 different DFT length,
256, 512, 1024, and 4096.

Here is my code:

%Part 1
% MATLAB code to analyze the frequency accuracy of the DFT manually
fs = 5000; % Sampling frequency 5 kHz
f_start = 1000; % Start frequency of sinusoid (1 kHz)
f_end = 1050; % End frequency of sinusoid (1.05 kHz)
f_step = 1; % Frequency step size (small increments)
t_duration = 1; % Signal duration in seconds
t = 0:1/fs:t_duration-1/fs; % Time vector
% DFT sizes to analyze
N_dft = [256, 512, 1024, 4096];
% Prepare figure
figure;
hold on;
for N = N_dft
freq_accuracy = []; % Store accuracy results for each DFT size
for f = f_start:f_step:f_end
% Generate sinusoidal waveform
x = sin(2*pi*f*t);

rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window;

% Manually compute the DFT
X = x_windowed*dftmtx(N);

% Compute the frequency axis for the current DFT
freq_axis = (0:N-1)*(fs/N);

% Find the index of the peak in the magnitude of the DFT
[~, peak_idx] = max(abs(X));

% Estimated frequency from the DFT
f_estimated = freq_axis(peak_idx);

% Store the frequency error (absolute difference between true and estimated frequency)
freq_error = abs(f_estimated – f);
freq_accuracy = [freq_accuracy, freq_error];
end

% Plot the frequency accuracy for the current DFT size
plot(f_start:f_step:f_end, freq_accuracy, ‘DisplayName’, [‘N = ‘, num2str(N)]);
end
% Add labels and legend
xlabel(‘True Frequency (Hz)’);
ylabel(‘Frequency Error (Hz)’);
title(‘Frequency Accuracy of the DFT for Different DFT Sizes (Manual DFT)’);
legend(‘show’);
grid on;
hold off;

However, the result I got is really weird.

At some frequency the estimate error between true frequency can be 3000, while most of the time it’s under 1.

However, it I substitue this line: X = x_windowed*dftmtx(N);
with this line: X = fft(x_windowed, N);

Then I got this image

It seems that this is the correct figure.
I read from this website that fft and dftmtx function can be interchangeble
https://www.mathworks.com/help/signal/ref/dftmtx.html
But my results doesn’t approve it. I want to know why.

Besides, am I applying windowing function in a correct way?
The code is as following:
rectangular_window = rectwin(N)’;
% Apply rectangular window (implicitly applied since no windowing function)
% Select only the first N samples for each DFT calculation
x_windowed = x(1:N).*rectangular_window; dft, discrete fourier transform, matlab, digital signal processing MATLAB Answers — New Questions

​

% Please note that this section requires the toolbox m_map
% Now we would like to know the mean states and annual trends of MHW
% frequency, i.e. how many MHW events would be detected per year and how it
% changes with time.
[mean_freq,annual_freq,trend_freq,p_freq]=mean_and_trend_new(MHW,mhw_ts,1982,’Metric’,’Frequency’);
% These four outputs separately represent the total mean, annual mean,
% annual trend and associated p value of frequency.
% This function could detect mean states and trends for six different
% variables (Frequency, mean intensity, max intensity, duration and total
% MHW/MCs days).
metric_used={‘Frequency’,’MeanInt’,’MaxInt’,’CumInt’,’Duration’,’Days’};
for i=1:6;
eval([‘[mean_’ metric_used{i} ‘,annual_’ metric_used{i} ‘,trend_’ metric_used{i} ‘,p_’ metric_used{i} ‘]=mean_and_trend_new(MHW,mhw_ts,1982,’ ”” ‘Metric’ ”” ‘,’ ‘metric_used{i}’ ‘);’])
end
% plot mean and trend
% It could be detected that, as a global hotspot, the oceanic region off
% eastern Tasmania exhibits significant positive trends of MHW metrics.
figure(‘pos’,[10 10 1500 1500]);
m_proj(‘mercator’,’lat’,[-45 -37],’lon’,[147 155]);
for i=1:6;
subplot(2,6,i);
eval([‘mean_here=mean_’ metric_used{i} ‘;’]);
eval([‘t_here=trend_’ metric_used{i} ‘;’]);
m_pcolor(lon_used,lat_used,mean_here’);
shading interp
m_coast(‘patch’,[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar(‘location’,’southoutside’);
title(metric_used{i},’fontname’,’consolas’,’fontsize’,12);

subplot(2,6,i+6);
eval([‘mean_here=mean_’ metric_used{i} ‘;’]);
eval([‘t_here=trend_’ metric_used{i} ‘;’]);
m_pcolor(lon_used,lat_used,t_here’);
shading interp
m_coast(‘patch’,[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar(‘location’,’southoutside’);
title([‘Trend-‘ metric_used{i}],’fontname’,’consolas’,’fontsize’,12);
end
%% 5. Applying cluster algoirthm to MHW – A kmeans example.
% We get so many MHWs now…. Could we distinguish them into different
% gropus based on their metrics?
% Change it to matrix;
MHW_m=MHW{:,:};
% Extract mean, max, cumulative intensity and duration.
MHW_m=MHW_m(:,[3 4 5 7]);
[data_for_k,mu,sd]=zscore(MHW_m);
% Determine suitable groups of kmeans cluster.
index_full=[];
cor_full=[];
for i=2:20;
k=kmeans(data_for_k,i,’Distance’,’cityblock’,’maxiter’,200);
k_full=[];
for j=1:i;
k_full=[k_full;nanmean(data_for_k(k==j,:))];
end

k_cor=k_full(k,:);
k_cor=k_cor(:);

[c,p]=corr([data_for_k(:) k_cor]);

index_full=[index_full;2];
cor_full=[cor_full;c(1,2)];

end
figure(‘pos’,[10 10 1500 1500]);
subplot(1,2,1);
plot(2:20,cor_full,’linewidth’,2);
hold on
plot(9*ones(1000,1),linspace(0.6,1,1000),’r–‘);
xlabel(‘Number of Groups’,’fontsize’,16,’fontweight’,’bold’);
ylabel(‘Correlation’,’fontsize’,16,’fontweight’,’bold’);
title(‘Correlation’,’fontsize’,16);
set(gca,’xtick’,[5 9 10 15 20],’fontsize’,16);
subplot(1,2,2);
plot(3:20,diff(cor_full),’linewidth’,2);
hold on
plot(9*ones(1000,1),linspace(-0.02,0.14,1000),’r–‘);
xlabel(‘Number of Groups’,’fontsize’,16,’fontweight’,’bold’);
ylabel(‘First difference of Correlation’,’fontsize’,16,’fontweight’,’bold’);
title(‘First Difference of Correlation’,’fontsize’,16);
set(gca,’fontsize’,16);
% Use 9 groups.
k=kmeans(data_for_k,9,’Distance’,’cityblock’,’maxiter’,200);
k_9=[];
prop_9=[];
for i=1:9;
data_here=data_for_k(k==i,:);
data_here=nanmean(data_here);
data_here=data_here.*sd+mu;
k_9=[k_9;data_here];
prop_9=[prop_9;nansum(k==i)./size(data_for_k,1)];
end
loc_x=[1.5 1.5 1.5 2.5 2.5 2.5 3.5 3.5 3.5];
loc_y=[1.5 2.5 3.5 1.5 2.5 3.5 1.5 2.5 3.5];
text_used={[‘1: ‘ num2str(round(prop_9(1)*100)) ‘%’],[‘2: ‘ num2str(round(prop_9(2)*100)) ‘%’],[‘3: ‘ num2str(round(prop_9(3)*100)) ‘%’],…
[‘4: ‘ num2str(round(prop_9(4)*100)) ‘%’],[‘5: ‘ num2str(round(prop_9(5)*100)) ‘%’],[‘6: ‘ num2str(round(prop_9(6)*100)) ‘%’],…
[‘7: ‘ num2str(round(prop_9(7)*100)) ‘%’],[‘8: ‘ num2str(round(prop_9(8)*100)) ‘%’],[‘9: ‘ num2str(round(prop_9(9)*100)) ‘%’]};
figure(‘pos’,[10 10 1500 1500]);
h=subplot(2,2,1);
data_here=k_9(:,1);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
set(h,’ydir’,’reverse’);
axis off
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar;
title(‘Durations’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,2);
data_here=k_9(:,2);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,’ydir’,’reverse’);
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘MaxInt’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,3);
data_here=k_9(:,3);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,’ydir’,’reverse’);
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘MeanInt’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,4);
data_here=k_9(:,4);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
[x,y]=meshgrid(1:4,1:4);
pcolor(1:4,1:4,data_here);
set(h,’ydir’,’reverse’);
axis off
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘CumInt’,’fontsize’,16,’fontweight’,’bold’);
% Their associated SSTA patterns
% Calculate SSTA
time_used=datevec(datenum(1982,1,1):datenum(2016,12,31));
m_d_unique=unique(time_used(:,2:3),’rows’);
ssta_full=NaN(size(sst_full));
for i=1:size(m_d_unique);
date_here=m_d_unique(i,:);
index_here=find(time_used(:,2)==date_here(1) & time_used(:,3)==date_here(2));
ssta_full(:,:,index_here)=sst_full(:,:,index_here)-nanmean(sst_full(:,:,index_here),3);
end
sst_1993_2016=ssta_full(:,:,(datenum(1993,1,1):datenum(2016,12,31))-datenum(1982,1,1)+1);
time_used=MHW{:,:};
time_used=time_used(:,1:2);
start_full=datenum(num2str(time_used(:,1)),’yyyymmdd’)-datenum(1993,1,1)+1;
end_full=datenum(num2str(time_used(:,2)),’yyyymmdd’)-datenum(1993,1,1)+1;
for i=1:9;
start_here=start_full(k==i);
end_here=end_full(k==i);

index_here=[];

for j=1:length(start_here);
period_here=start_here(j):end_here(j);
index_here=[index_here;period_here(:)];
end

eval([‘sst_’ num2str(i) ‘=nanmean(sst_1993_2016(:,:,index_here),3);’])

end
color_used=hot;
color_used=color_used(end:-1:1,:);
figure(‘pos’,[10 10 1500 1500]);
plot_index=[1 4 7 2 5 8 3 6 9];
for i=1:9;
eval([‘plot_here=sst_’ num2str(i) ‘;’]);
subplot(3,3,plot_index(i));
eval([‘data_here=sst_’ num2str(i) ‘;’])

m_contourf(lon_used,lat_used,data_here’,-3:0.1:3,’linestyle’,’none’);

if i~=3;
m_grid(‘xtick’,[],’ytick’,[]);
else;
m_grid(‘linestyle’,’none’);
end
m_gshhs_h(‘patch’,[0 0 0]);
colormap(color_used);
caxis([0 2]);

title([‘Group (‘ num2str(i) ‘):’ num2str(round(prop_9(i)*100)) ‘%’],’fontsize’,16);

end
hp4=get(subplot(3,3,9),’Position’);
s=colorbar(‘Position’, [hp4(1)+hp4(3)+0.025 hp4(2) 0.025 0.85],’fontsize’,14);
s.Label.String=’^{o}C’;
Index exceeds the number of array elements. Index must not exceed 12784.
Error in event_line (line 84)
y1=m90_plot(x1-datenum(data_start,1,1)+1);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in an_example2 (line 53)
event_line(sst_full,MHW,mclim,m90,[1 2],1982,[2015 9 1],[2016 5 1]);% Please note that this section requires the toolbox m_map
% Now we would like to know the mean states and annual trends of MHW
% frequency, i.e. how many MHW events would be detected per year and how it
% changes with time.
[mean_freq,annual_freq,trend_freq,p_freq]=mean_and_trend_new(MHW,mhw_ts,1982,’Metric’,’Frequency’);
% These four outputs separately represent the total mean, annual mean,
% annual trend and associated p value of frequency.
% This function could detect mean states and trends for six different
% variables (Frequency, mean intensity, max intensity, duration and total
% MHW/MCs days).
metric_used={‘Frequency’,’MeanInt’,’MaxInt’,’CumInt’,’Duration’,’Days’};
for i=1:6;
eval([‘[mean_’ metric_used{i} ‘,annual_’ metric_used{i} ‘,trend_’ metric_used{i} ‘,p_’ metric_used{i} ‘]=mean_and_trend_new(MHW,mhw_ts,1982,’ ”” ‘Metric’ ”” ‘,’ ‘metric_used{i}’ ‘);’])
end
% plot mean and trend
% It could be detected that, as a global hotspot, the oceanic region off
% eastern Tasmania exhibits significant positive trends of MHW metrics.
figure(‘pos’,[10 10 1500 1500]);
m_proj(‘mercator’,’lat’,[-45 -37],’lon’,[147 155]);
for i=1:6;
subplot(2,6,i);
eval([‘mean_here=mean_’ metric_used{i} ‘;’]);
eval([‘t_here=trend_’ metric_used{i} ‘;’]);
m_pcolor(lon_used,lat_used,mean_here’);
shading interp
m_coast(‘patch’,[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar(‘location’,’southoutside’);
title(metric_used{i},’fontname’,’consolas’,’fontsize’,12);

subplot(2,6,i+6);
eval([‘mean_here=mean_’ metric_used{i} ‘;’]);
eval([‘t_here=trend_’ metric_used{i} ‘;’]);
m_pcolor(lon_used,lat_used,t_here’);
shading interp
m_coast(‘patch’,[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar(‘location’,’southoutside’);
title([‘Trend-‘ metric_used{i}],’fontname’,’consolas’,’fontsize’,12);
end
%% 5. Applying cluster algoirthm to MHW – A kmeans example.
% We get so many MHWs now…. Could we distinguish them into different
% gropus based on their metrics?
% Change it to matrix;
MHW_m=MHW{:,:};
% Extract mean, max, cumulative intensity and duration.
MHW_m=MHW_m(:,[3 4 5 7]);
[data_for_k,mu,sd]=zscore(MHW_m);
% Determine suitable groups of kmeans cluster.
index_full=[];
cor_full=[];
for i=2:20;
k=kmeans(data_for_k,i,’Distance’,’cityblock’,’maxiter’,200);
k_full=[];
for j=1:i;
k_full=[k_full;nanmean(data_for_k(k==j,:))];
end

k_cor=k_full(k,:);
k_cor=k_cor(:);

[c,p]=corr([data_for_k(:) k_cor]);

index_full=[index_full;2];
cor_full=[cor_full;c(1,2)];

end
figure(‘pos’,[10 10 1500 1500]);
subplot(1,2,1);
plot(2:20,cor_full,’linewidth’,2);
hold on
plot(9*ones(1000,1),linspace(0.6,1,1000),’r–‘);
xlabel(‘Number of Groups’,’fontsize’,16,’fontweight’,’bold’);
ylabel(‘Correlation’,’fontsize’,16,’fontweight’,’bold’);
title(‘Correlation’,’fontsize’,16);
set(gca,’xtick’,[5 9 10 15 20],’fontsize’,16);
subplot(1,2,2);
plot(3:20,diff(cor_full),’linewidth’,2);
hold on
plot(9*ones(1000,1),linspace(-0.02,0.14,1000),’r–‘);
xlabel(‘Number of Groups’,’fontsize’,16,’fontweight’,’bold’);
ylabel(‘First difference of Correlation’,’fontsize’,16,’fontweight’,’bold’);
title(‘First Difference of Correlation’,’fontsize’,16);
set(gca,’fontsize’,16);
% Use 9 groups.
k=kmeans(data_for_k,9,’Distance’,’cityblock’,’maxiter’,200);
k_9=[];
prop_9=[];
for i=1:9;
data_here=data_for_k(k==i,:);
data_here=nanmean(data_here);
data_here=data_here.*sd+mu;
k_9=[k_9;data_here];
prop_9=[prop_9;nansum(k==i)./size(data_for_k,1)];
end
loc_x=[1.5 1.5 1.5 2.5 2.5 2.5 3.5 3.5 3.5];
loc_y=[1.5 2.5 3.5 1.5 2.5 3.5 1.5 2.5 3.5];
text_used={[‘1: ‘ num2str(round(prop_9(1)*100)) ‘%’],[‘2: ‘ num2str(round(prop_9(2)*100)) ‘%’],[‘3: ‘ num2str(round(prop_9(3)*100)) ‘%’],…
[‘4: ‘ num2str(round(prop_9(4)*100)) ‘%’],[‘5: ‘ num2str(round(prop_9(5)*100)) ‘%’],[‘6: ‘ num2str(round(prop_9(6)*100)) ‘%’],…
[‘7: ‘ num2str(round(prop_9(7)*100)) ‘%’],[‘8: ‘ num2str(round(prop_9(8)*100)) ‘%’],[‘9: ‘ num2str(round(prop_9(9)*100)) ‘%’]};
figure(‘pos’,[10 10 1500 1500]);
h=subplot(2,2,1);
data_here=k_9(:,1);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
set(h,’ydir’,’reverse’);
axis off
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar;
title(‘Durations’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,2);
data_here=k_9(:,2);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,’ydir’,’reverse’);
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘MaxInt’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,3);
data_here=k_9(:,3);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,’ydir’,’reverse’);
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘MeanInt’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,4);
data_here=k_9(:,4);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
[x,y]=meshgrid(1:4,1:4);
pcolor(1:4,1:4,data_here);
set(h,’ydir’,’reverse’);
axis off
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘CumInt’,’fontsize’,16,’fontweight’,’bold’);
% Their associated SSTA patterns
% Calculate SSTA
time_used=datevec(datenum(1982,1,1):datenum(2016,12,31));
m_d_unique=unique(time_used(:,2:3),’rows’);
ssta_full=NaN(size(sst_full));
for i=1:size(m_d_unique);
date_here=m_d_unique(i,:);
index_here=find(time_used(:,2)==date_here(1) & time_used(:,3)==date_here(2));
ssta_full(:,:,index_here)=sst_full(:,:,index_here)-nanmean(sst_full(:,:,index_here),3);
end
sst_1993_2016=ssta_full(:,:,(datenum(1993,1,1):datenum(2016,12,31))-datenum(1982,1,1)+1);
time_used=MHW{:,:};
time_used=time_used(:,1:2);
start_full=datenum(num2str(time_used(:,1)),’yyyymmdd’)-datenum(1993,1,1)+1;
end_full=datenum(num2str(time_used(:,2)),’yyyymmdd’)-datenum(1993,1,1)+1;
for i=1:9;
start_here=start_full(k==i);
end_here=end_full(k==i);

index_here=[];

for j=1:length(start_here);
period_here=start_here(j):end_here(j);
index_here=[index_here;period_here(:)];
end

eval([‘sst_’ num2str(i) ‘=nanmean(sst_1993_2016(:,:,index_here),3);’])

end
color_used=hot;
color_used=color_used(end:-1:1,:);
figure(‘pos’,[10 10 1500 1500]);
plot_index=[1 4 7 2 5 8 3 6 9];
for i=1:9;
eval([‘plot_here=sst_’ num2str(i) ‘;’]);
subplot(3,3,plot_index(i));
eval([‘data_here=sst_’ num2str(i) ‘;’])

m_contourf(lon_used,lat_used,data_here’,-3:0.1:3,’linestyle’,’none’);

if i~=3;
m_grid(‘xtick’,[],’ytick’,[]);
else;
m_grid(‘linestyle’,’none’);
end
m_gshhs_h(‘patch’,[0 0 0]);
colormap(color_used);
caxis([0 2]);

title([‘Group (‘ num2str(i) ‘):’ num2str(round(prop_9(i)*100)) ‘%’],’fontsize’,16);

end
hp4=get(subplot(3,3,9),’Position’);
s=colorbar(‘Position’, [hp4(1)+hp4(3)+0.025 hp4(2) 0.025 0.85],’fontsize’,14);
s.Label.String=’^{o}C’;
Index exceeds the number of array elements. Index must not exceed 12784.
Error in event_line (line 84)
y1=m90_plot(x1-datenum(data_start,1,1)+1);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in an_example2 (line 53)
event_line(sst_full,MHW,mclim,m90,[1 2],1982,[2015 9 1],[2016 5 1]); % Please note that this section requires the toolbox m_map
% Now we would like to know the mean states and annual trends of MHW
% frequency, i.e. how many MHW events would be detected per year and how it
% changes with time.
[mean_freq,annual_freq,trend_freq,p_freq]=mean_and_trend_new(MHW,mhw_ts,1982,’Metric’,’Frequency’);
% These four outputs separately represent the total mean, annual mean,
% annual trend and associated p value of frequency.
% This function could detect mean states and trends for six different
% variables (Frequency, mean intensity, max intensity, duration and total
% MHW/MCs days).
metric_used={‘Frequency’,’MeanInt’,’MaxInt’,’CumInt’,’Duration’,’Days’};
for i=1:6;
eval([‘[mean_’ metric_used{i} ‘,annual_’ metric_used{i} ‘,trend_’ metric_used{i} ‘,p_’ metric_used{i} ‘]=mean_and_trend_new(MHW,mhw_ts,1982,’ ”” ‘Metric’ ”” ‘,’ ‘metric_used{i}’ ‘);’])
end
% plot mean and trend
% It could be detected that, as a global hotspot, the oceanic region off
% eastern Tasmania exhibits significant positive trends of MHW metrics.
figure(‘pos’,[10 10 1500 1500]);
m_proj(‘mercator’,’lat’,[-45 -37],’lon’,[147 155]);
for i=1:6;
subplot(2,6,i);
eval([‘mean_here=mean_’ metric_used{i} ‘;’]);
eval([‘t_here=trend_’ metric_used{i} ‘;’]);
m_pcolor(lon_used,lat_used,mean_here’);
shading interp
m_coast(‘patch’,[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar(‘location’,’southoutside’);
title(metric_used{i},’fontname’,’consolas’,’fontsize’,12);

subplot(2,6,i+6);
eval([‘mean_here=mean_’ metric_used{i} ‘;’]);
eval([‘t_here=trend_’ metric_used{i} ‘;’]);
m_pcolor(lon_used,lat_used,t_here’);
shading interp
m_coast(‘patch’,[0.7 0.7 0.7]);
m_grid;
colormap(jet);
s=colorbar(‘location’,’southoutside’);
title([‘Trend-‘ metric_used{i}],’fontname’,’consolas’,’fontsize’,12);
end
%% 5. Applying cluster algoirthm to MHW – A kmeans example.
% We get so many MHWs now…. Could we distinguish them into different
% gropus based on their metrics?
% Change it to matrix;
MHW_m=MHW{:,:};
% Extract mean, max, cumulative intensity and duration.
MHW_m=MHW_m(:,[3 4 5 7]);
[data_for_k,mu,sd]=zscore(MHW_m);
% Determine suitable groups of kmeans cluster.
index_full=[];
cor_full=[];
for i=2:20;
k=kmeans(data_for_k,i,’Distance’,’cityblock’,’maxiter’,200);
k_full=[];
for j=1:i;
k_full=[k_full;nanmean(data_for_k(k==j,:))];
end

k_cor=k_full(k,:);
k_cor=k_cor(:);

[c,p]=corr([data_for_k(:) k_cor]);

index_full=[index_full;2];
cor_full=[cor_full;c(1,2)];

end
figure(‘pos’,[10 10 1500 1500]);
subplot(1,2,1);
plot(2:20,cor_full,’linewidth’,2);
hold on
plot(9*ones(1000,1),linspace(0.6,1,1000),’r–‘);
xlabel(‘Number of Groups’,’fontsize’,16,’fontweight’,’bold’);
ylabel(‘Correlation’,’fontsize’,16,’fontweight’,’bold’);
title(‘Correlation’,’fontsize’,16);
set(gca,’xtick’,[5 9 10 15 20],’fontsize’,16);
subplot(1,2,2);
plot(3:20,diff(cor_full),’linewidth’,2);
hold on
plot(9*ones(1000,1),linspace(-0.02,0.14,1000),’r–‘);
xlabel(‘Number of Groups’,’fontsize’,16,’fontweight’,’bold’);
ylabel(‘First difference of Correlation’,’fontsize’,16,’fontweight’,’bold’);
title(‘First Difference of Correlation’,’fontsize’,16);
set(gca,’fontsize’,16);
% Use 9 groups.
k=kmeans(data_for_k,9,’Distance’,’cityblock’,’maxiter’,200);
k_9=[];
prop_9=[];
for i=1:9;
data_here=data_for_k(k==i,:);
data_here=nanmean(data_here);
data_here=data_here.*sd+mu;
k_9=[k_9;data_here];
prop_9=[prop_9;nansum(k==i)./size(data_for_k,1)];
end
loc_x=[1.5 1.5 1.5 2.5 2.5 2.5 3.5 3.5 3.5];
loc_y=[1.5 2.5 3.5 1.5 2.5 3.5 1.5 2.5 3.5];
text_used={[‘1: ‘ num2str(round(prop_9(1)*100)) ‘%’],[‘2: ‘ num2str(round(prop_9(2)*100)) ‘%’],[‘3: ‘ num2str(round(prop_9(3)*100)) ‘%’],…
[‘4: ‘ num2str(round(prop_9(4)*100)) ‘%’],[‘5: ‘ num2str(round(prop_9(5)*100)) ‘%’],[‘6: ‘ num2str(round(prop_9(6)*100)) ‘%’],…
[‘7: ‘ num2str(round(prop_9(7)*100)) ‘%’],[‘8: ‘ num2str(round(prop_9(8)*100)) ‘%’],[‘9: ‘ num2str(round(prop_9(9)*100)) ‘%’]};
figure(‘pos’,[10 10 1500 1500]);
h=subplot(2,2,1);
data_here=k_9(:,1);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
set(h,’ydir’,’reverse’);
axis off
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar;
title(‘Durations’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,2);
data_here=k_9(:,2);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,’ydir’,’reverse’);
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘MaxInt’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,3);
data_here=k_9(:,3);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
pcolor(1:4,1:4,data_here);
axis off
set(h,’ydir’,’reverse’);
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘MeanInt’,’fontsize’,16,’fontweight’,’bold’);
h=subplot(2,2,4);
data_here=k_9(:,4);
data_here=reshape(data_here,3,3);
data_here(:,end+1)=data_here(:,end);
data_here(end+1,:)=data_here(end,:);
[x,y]=meshgrid(1:4,1:4);
pcolor(1:4,1:4,data_here);
set(h,’ydir’,’reverse’);
axis off
colormap(jet);
text(loc_x,loc_y,text_used,’fontsize’,16,’horiz’,’center’,’fontweight’,’bold’);
colorbar
title(‘CumInt’,’fontsize’,16,’fontweight’,’bold’);
% Their associated SSTA patterns
% Calculate SSTA
time_used=datevec(datenum(1982,1,1):datenum(2016,12,31));
m_d_unique=unique(time_used(:,2:3),’rows’);
ssta_full=NaN(size(sst_full));
for i=1:size(m_d_unique);
date_here=m_d_unique(i,:);
index_here=find(time_used(:,2)==date_here(1) & time_used(:,3)==date_here(2));
ssta_full(:,:,index_here)=sst_full(:,:,index_here)-nanmean(sst_full(:,:,index_here),3);
end
sst_1993_2016=ssta_full(:,:,(datenum(1993,1,1):datenum(2016,12,31))-datenum(1982,1,1)+1);
time_used=MHW{:,:};
time_used=time_used(:,1:2);
start_full=datenum(num2str(time_used(:,1)),’yyyymmdd’)-datenum(1993,1,1)+1;
end_full=datenum(num2str(time_used(:,2)),’yyyymmdd’)-datenum(1993,1,1)+1;
for i=1:9;
start_here=start_full(k==i);
end_here=end_full(k==i);

index_here=[];

for j=1:length(start_here);
period_here=start_here(j):end_here(j);
index_here=[index_here;period_here(:)];
end

eval([‘sst_’ num2str(i) ‘=nanmean(sst_1993_2016(:,:,index_here),3);’])

end
color_used=hot;
color_used=color_used(end:-1:1,:);
figure(‘pos’,[10 10 1500 1500]);
plot_index=[1 4 7 2 5 8 3 6 9];
for i=1:9;
eval([‘plot_here=sst_’ num2str(i) ‘;’]);
subplot(3,3,plot_index(i));
eval([‘data_here=sst_’ num2str(i) ‘;’])

m_contourf(lon_used,lat_used,data_here’,-3:0.1:3,’linestyle’,’none’);

if i~=3;
m_grid(‘xtick’,[],’ytick’,[]);
else;
m_grid(‘linestyle’,’none’);
end
m_gshhs_h(‘patch’,[0 0 0]);
colormap(color_used);
caxis([0 2]);

title([‘Group (‘ num2str(i) ‘):’ num2str(round(prop_9(i)*100)) ‘%’],’fontsize’,16);

end
hp4=get(subplot(3,3,9),’Position’);
s=colorbar(‘Position’, [hp4(1)+hp4(3)+0.025 hp4(2) 0.025 0.85],’fontsize’,14);
s.Label.String=’^{o}C’;
Index exceeds the number of array elements. Index must not exceed 12784.
Error in event_line (line 84)
y1=m90_plot(x1-datenum(data_start,1,1)+1);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in an_example2 (line 53)
event_line(sst_full,MHW,mclim,m90,[1 2],1982,[2015 9 1],[2016 5 1]); error, heat MATLAB Answers — New Questions

​

hi, I want to choose as value in a table field ‘Fil’ or ‘Stat’

cat=categorical({‘Fil’;’Stat’});
name={‘A’,’B’};
nrows=numel(name);
VNAMES={‘Draw’;’Filt_Stat’};
VTYPES=[{‘logical’},{cat}];
T=table(‘Size’,[nrows,numel(VNAMES)],’VariableTypes’,VTYPES,’VariableNames’,VNAMES); % empty t
T.Draw=repmat({true},nrows,1);
ct=cell(1, numel(name));
ct(:) = {cat{2}};
T.Filt_Stat=ct’;
app.Portfolio_UITable.Data=T;

i get this error:
Error using table
Specify variable types as a string array or a cell array of character
vectors, such as ["string", "datetime", "double"].hi, I want to choose as value in a table field ‘Fil’ or ‘Stat’

cat=categorical({‘Fil’;’Stat’});
name={‘A’,’B’};
nrows=numel(name);
VNAMES={‘Draw’;’Filt_Stat’};
VTYPES=[{‘logical’},{cat}];
T=table(‘Size’,[nrows,numel(VNAMES)],’VariableTypes’,VTYPES,’VariableNames’,VNAMES); % empty t
T.Draw=repmat({true},nrows,1);
ct=cell(1, numel(name));
ct(:) = {cat{2}};
T.Filt_Stat=ct’;
app.Portfolio_UITable.Data=T;

i get this error:
Error using table
Specify variable types as a string array or a cell array of character
vectors, such as ["string", "datetime", "double"]. hi, I want to choose as value in a table field ‘Fil’ or ‘Stat’

cat=categorical({‘Fil’;’Stat’});
name={‘A’,’B’};
nrows=numel(name);
VNAMES={‘Draw’;’Filt_Stat’};
VTYPES=[{‘logical’},{cat}];
T=table(‘Size’,[nrows,numel(VNAMES)],’VariableTypes’,VTYPES,’VariableNames’,VNAMES); % empty t
T.Draw=repmat({true},nrows,1);
ct=cell(1, numel(name));
ct(:) = {cat{2}};
T.Filt_Stat=ct’;
app.Portfolio_UITable.Data=T;

i get this error:
Error using table
Specify variable types as a string array or a cell array of character
vectors, such as ["string", "datetime", "double"]. how to use categorical in uitable MATLAB Answers — New Questions

​

Hi,
I am having a weird issue with the Arduino via Matlab code.
this is the simple code I have
a = arduino
configurePin(a,’D12′,’PWM’)
writePWMVoltage(a,’D12′,3)
writePWMDutyCycle(a,’D12′,0.5)
When I run this code, it is just outputting 1.65V which is 0.5 of 3.3V (Due) instead of doing the cycle of the voltage.
Anyone made this cycle to work?Hi,
I am having a weird issue with the Arduino via Matlab code.
this is the simple code I have
a = arduino
configurePin(a,’D12′,’PWM’)
writePWMVoltage(a,’D12′,3)
writePWMDutyCycle(a,’D12′,0.5)
When I run this code, it is just outputting 1.65V which is 0.5 of 3.3V (Due) instead of doing the cycle of the voltage.
Anyone made this cycle to work? Hi,
I am having a weird issue with the Arduino via Matlab code.
this is the simple code I have
a = arduino
configurePin(a,’D12′,’PWM’)
writePWMVoltage(a,’D12′,3)
writePWMDutyCycle(a,’D12′,0.5)
When I run this code, it is just outputting 1.65V which is 0.5 of 3.3V (Due) instead of doing the cycle of the voltage.
Anyone made this cycle to work? arduino, matlab, pwm MATLAB Answers — New Questions

​

How do I name a variable in simulink so that i can later use it to compute a value based on it, such sin(psi) or cos(phi) or tan(theta)?How do I name a variable in simulink so that i can later use it to compute a value based on it, such sin(psi) or cos(phi) or tan(theta)? How do I name a variable in simulink so that i can later use it to compute a value based on it, such sin(psi) or cos(phi) or tan(theta)? simulink names variables MATLAB Answers — New Questions

​