alpha = 45; beta = 185; gamma_e = 116; G_ei = -4.11; G_ee = 2.07;G_sr = -3.30;G_rs = 0.20;G_es = 0.77;G_re = 0.66;G_se = 7.77;G_sn = 8.10;G_esre = G_es*G_sr*G_re;G_srs = G_sr*G_rs;G_ese = G_es*G_se;G_esn = G_es*G_sn;t_0 = 0.085; r_e = 0.086; a = -60; b = 60; n = 300; f = -40:40 w = 2*pi*f; % angular frequency in radian per second
for k = 1:length(w) L_initial = @(w1) (1-((1i*w1)/alpha))^(-1)* (1-((1i*w1)/beta))^(-1); L_shift = @(w1) (1-((1i*(w(k)-w1))/alpha))^(-1)* (1-((1i*(w(k)-w1))/beta))^(-1); L = @(w1) L_initial(w1)*L_shift(w1); low pass filter for w1*(w-w1)Q_initial = @(w1) (1/(r_e^2))*((1-((1i*w1)/gamma_e))^(2) - (1/(1-G_ei*L_initial(w1)))* (L_initial(w1)*G_ee + (exp(1i*w1*t_0)*(L_initial(w1)^2*G_ese +L_initial(w1).^3*G_esre))/(1-L_initial(w1)^2*G_srs))); Q_shift = @(w1) (1/(r_e^2))*((1-((1i*(w(k)-w1))/gamma_e))^(2) - (1/(1-G_ei*L_shift(w1)))* (L_shift(w1)*G_ee + (exp(1i*(w(k)-w1)*t_0)*(L_shift(w1)^2*G_ese +L_shift(w1)^3*G_esre))/(1-L_shift(w1)^2*G_srs))); Q = @(w1) Q_initial(w1)*Q_shift(w1); p_initial = @(w1) (pi/r_e^4)* (abs((L_initial(w1)^2*G_esn)/((1-L_initial(w1)^2*G_srs)*(1-G_ei*L_initial(w1))))).^2 * abs((atan2((imag(Q_initial(w1))),(real(Q_initial(w1)))))/imag(Q_initial(w1))); p_shift = @(w1) (pi/r_e^4)* (abs((L_shift(w1)^2*G_esn)/((1-L_shift(w1)^2*G_srs)*(1-G_ei*L_shift(w1))))).^2 * abs((atan2((imag(Q_shift(w1))),(real(Q_shift(w1)))))/imag(Q_shift(w1))); p = @(w1) p_initial(w1)*p_shift(w1); I(k) = simprl(p,a,b,n) (simprl is a calling function) endplot(f,I)xlabel('f (frequency in Hz)') ylabel('I (numerical integral values)')
%%%%%%%%%
function [s] = simprl(p,a,b,n) h = (b-a)./n;s1 = 0; s2=0; % loop for odd values in the range
for k = 1:n/2; x = a + h*(2*k-1); s1 = s1+feval(p,x); end% loop for even values in the range
for k = 1:(n/2 - 1); x = a + h*2*k; s2 = s2+feval(p,x); ends = h*(feval(p,a)+feval(p,b)+4*s1+2*s2)/3;
Best Answer