Contour maps matlab for a number of poles

Hi with this code:

damp_theta_vec = [0.1:0.1:0.6];

wn_theta_aux = [-3:0.1:-0.1];

for k=1:length(wn_theta_aux)

for j = 1:length(damp_theta_vec)

wn_theta = wn_theta_aux(k)/(-damp_theta_vec(j));

damp_theta = damp_theta_vec(j);

couple_aux(k,j) = wn_theta;

sys_theta = tf(wn_theta^2, [1 2*damp_theta*wn_theta wn_theta^2 0]);

%% Actuator parameters

wn_act = 150; %[rad/s]

damp_act = 0.2; %[]

max_def = 30*pi/180;

min_def = -30*pi/180;

rate_limit_act = 100*pi/180; %[rad/s]

%% Criterio de short-term response (Figura 5)

%obtenemos datos del bode de la función de transferencia

[mag_theta,phase_theta,wout_theta]=bode(sys_theta,{0.01,100});

%encontramos la frequencia y la magnitud donde la fase es -180

for i = 1:size(wout_theta,1)

if phase_theta(i)<-180

w_180 = wout_theta(i);

mag_180 = mag_theta(i);

break

end

end

w_180_out(k,j) = w_180;

mag_180_out(k,j) = mag_180;

%encontramos la frecuencia donde la fase es -180 + 45 = -135

for i = 1:size(wout_theta,1)

if phase_theta(i)<-135

w_bd_phase = wout_theta(i);

break

end

end

w_bd_phase_out(k,j) = w_bd_phase;

%encontramos la fase para la frecuencia 2 x w_180 y calculamos el Phase Delay

for i = 1:size(wout_theta,1)

if wout_theta(i)>2*w_180

delta_phase_aux = -180 - phase_theta(i);

phase_delay_theta = delta_phase_aux/57.3/(2*w_180);

break

end

end

phase_delay_theta_out(k,j)=phase_delay_theta;

%encontramos la frecuencia que tiene 6dB más que la magnitud a w_180

for i = 1:size(wout_theta,1)

if (mag_180+6)>mag_theta(i)

w_bd_mag = wout_theta(i);

break

end

end

w_bd_mag_out(k,j) = w_bd_mag;

%criterio de aceptación w_bd_mag < w_bd_phase (Figura 6) ¿Cumple?

if (w_bd_mag < w_bd_phase)

disp('Criterio de aceptación en respuesta en el corto plazo se cumple');

crit_1(k,j) = 1;

else

disp('Criterio de aceptación respuesta en el corto plazo no se cumple');

crit_1(k,j)=0;

end

%criterio de aceptación w_bd (w_bd = w_bd_phase) y phase_delay_theta de acuerdo con Figura 5 ¿Cumple?

a_aux = (0.4-0.17)/(4.5-2);

b_aux = 0.17-a_aux*2;

if (w_bd_phase >= 2.0) && (phase_delay_theta <= 0.17)

disp('Bandwidth pertenece al Level 1');

crit_2(k,j)=1;

else if (w_bd_phase >= 2.0) && (w_bd_phase <= 4.0) && (phase_delay_theta <= a_aux*w_bd_phase+b_aux)

disp('Bandwidth pertenece al Level 1');

crit_2(k,j)=1;

else

disp('Bandwidth no pertenece al Level 1');

crit_2(k,j)=0;

end

end

%% Criterio de sobreoscilación (Requirements for moderate-amplitude) (Figura 8)

%obtenemos q_pk (derivamos la función de transferencia multiplicando por s)

q_pk = max(impulse(sys_theta*tf([1 0],1)));

q_pk_out(k,j) = q_pk;

%obtenemos Theta_pk

[Theta_pk,i_aux] = max(impulse(sys_theta));

Theta_pk_out(k,j) = Theta_pk;

[val_aux,t_aux] = impulse(sys_theta);

Theta_min = min(val_aux(i_aux:end));

peak = q_pk/Theta_pk;

peak_out(k,j) = peak;

%caso de función de transferencia y modelo lineal.

%Criterio de aceptación sobreoscilación Figura 8 ¿Cumple?

if (peak >=1.55)

disp('Criterio de sobreoscilación. Pertenece al Level 1');

crit_3(k,j)=1;

else

disp('Criterio de sobreoscilación. No pertenece al Level 1');

crit_3(k,j) =0;

end

end

end

How can I create a grid in the complex plane. With the pole that is in a specific node of said grid, calculate the dynamic response of the helicopter, and find magnitudes that characterize said response (amplitudes, angular velocities, ...). Do it for each point of the grid. Create several contour maps where, for each position of the pole in the complex plane, draw these magnitudes. Obtain the information that justifies which is the most suitable location for the pole, from the point of view that it meets the criteria.