How do I plot these two variables when the values change and graph is empty

I'm very sorry the code is so long but ultimately all I'm trying to do is vary the g_m and plot F vs g_m. To see the effect. When I try and plot the value the graph appears empty. Any ideas would be much appreciated

d=0.413/39.37; %thickness of aluminium
tau=0.099; %pole pitch
m=3; %number of phases
mu0=4*pi*10^-7;
q=3; %number of slots per pole per phase
p=6; %poles
W_s=0.05; %width of primary
N_c=4; %number of turns per slot
k_f=0.7; %fill factor - not used???
% calculated from wire cross section * turns / (slot width * slot height)
w_s=0.0061; %slot width
h_s=0.016; %slot height
k_p=1; %Chording factor. This is =y/tau, where y is the disance between
%windings. Thus, in a perfect machine, y=tau, and the chording factor =1
lambda=0.0111; %slot pitch
wire_width=0.004; %width of flat wire
wire_height=0.002; %height of flat wire
l_ce=0.05; %length of end connection due to the winding structure
A_w=wire_width*wire_height; %cross-sectional area of flat wire
g_m=0.003
while g_m<0.01
 g_m=g_m+0.001
%g_m=0.003; %mechanical airgap
rho_w=1.667*10^-8; %resistivity of copper in COMSOL
rho_r=2.649*10^-8; %resistivity of aluminium in COMSOL
K_w=0.75;%approximation from book was 0.9. This has been adjusted to 0.75
f=262; %frequency
v_r=42; %LIM velocity
I1=400; %primary phase current

%%%%%%%%%%%%%%%%%%%%%%%%% EQUIVALENT CIRCUIT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Primary winding resistance, R1
l_w1=2*(W_s+l_ce);
N_1=N_c*q*p/2;
l_w=N_1*l_w1;
J=I1/A_w;
R1=rho_w*l_w/A_w;
%%%Primary leakage reactance, X1
g=g_m+d/2;%???
lambda_s=(h_s*(1+(3*k_p)))/(12*w_s);
lambda_d=(5*g/w_s)/(5+(4*g/w_s));
lambda_e=0.3*((3*k_p)-1);
bracket=(lambda_s+lambda_d)*(W_s/q);
X1s=(8*mu0*pi*f/p)*(bracket)*(N_1^2);
X1e=(8*mu0*pi*f/p)*(lambda_e*l_ce)*(N_1^2);
X1=X1s+X1e;
X1_check=(8*mu0*pi*f/p)*(bracket+(lambda_e*l_ce))*(N_1^2);
%%%Magnetising reactuance, Xm
alpha=pi/(m*q);
k_d=(sin(q*alpha/2))/(q*(sin(alpha/2)));
constant=w_s/(2*g);
gamma=(4/pi)*((constant*atan(constant))-log(sqrt(1+(constant^2))));
k_c=lambda/(lambda-(gamma*g));
g_e=k_c*g;
W_se=W_s+g_e;
X_m=((24*mu0*f*W_se)/(pi*p*g_e))*((K_w*N_1)^2)*tau;
%%%Secondary resistance, R2
G=(2*mu0*f*(tau^2))/((rho_r/d)*pi*g_e);
R2=X_m/G;
%%%%%%%%%%%%%%%%%%%%% STEADY-STATE CHARACTERISTICS %%%%%%%%%%%%%%%%%%%%%%%%
%%%Total Impedance
v_s=2*f*tau;
s=1-(v_r/v_s);
Z=R1+(X1*j)+((j*X_m*(R2/s))/((j*X_m)+(R2/s)));
Z_real=real(Z);
Z_imag=imag(Z);
Z_mag=sqrt((Z_real^2)+(Z_imag^2));
phi=atan(Z_imag/Z_real);
%%%Primary phase voltage, V1
V1=I1*Z_mag;
%%%Secondary phase current, I2
I2=(X_m/(sqrt(((R2/s)^2)+(X_m^2))))*I1;
%%%Magnetising current, Im
Im=((R2/s)/(sqrt(((R2/s)^2)+(X_m^2))))*I1;
%%%Input active power, P_input
P_input=m*V1*I1*cos(phi);
%%%Output power, P_output
P_output=P_input-(m*(I1^2)*R1)-(m*(I2^2)*R2);
%%%Electromagnetic Force, F
F=P_output/v_r
F_check=(3*(I1^2)*R2)/(2*s*f*tau*(((1/(s*G))^2)+1));
%%%Magnetomotive Force, MMF
MMF=(4*sqrt(2)*m*K_w*N_1*Im)/(pi*p);
%%%Efficiency, eff
eff=P_output/P_input;
end
plot(g_m,F)

Answer 1

[Ameer Hamza]

You are not saving values inside the while-loop. Here is a simple way to solve the issue

d=0.413/39.37; %thickness of aluminium
tau=0.099; %pole pitch
m=3; %number of phases
mu0=4*pi*10^-7;
q=3; %number of slots per pole per phase
p=6; %poles
W_s=0.05; %width of primary
N_c=4; %number of turns per slot
k_f=0.7; %fill factor - not used???
% calculated from wire cross section * turns / (slot width * slot height)
w_s=0.0061; %slot width
h_s=0.016; %slot height
k_p=1; %Chording factor. This is =y/tau, where y is the disance between
%windings. Thus, in a perfect machine, y=tau, and the chording factor =1
lambda=0.0111; %slot pitch
wire_width=0.004; %width of flat wire
wire_height=0.002; %height of flat wire
l_ce=0.05; %length of end connection due to the winding structure
A_w=wire_width*wire_height; %cross-sectional area of flat wire
g_m=0.003;
F = [];
while g_m<0.01
g_m(end+1)=g_m(end)+0.001
%g_m=0.003; %mechanical airgap
rho_w=1.667*10^-8; %resistivity of copper in COMSOL
rho_r=2.649*10^-8; %resistivity of aluminium in COMSOL
K_w=0.75;%approximation from book was 0.9. This has been adjusted to 0.75
f=262; %frequency
v_r=42; %LIM velocity
I1=400; %primary phase current

%%%%%%%%%%%%%%%%%%%%%%%%% EQUIVALENT CIRCUIT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Primary winding resistance, R1
l_w1=2*(W_s+l_ce);
N_1=N_c*q*p/2;
l_w=N_1*l_w1;
J=I1/A_w;
R1=rho_w*l_w/A_w;
%%%Primary leakage reactance, X1
g=g_m(end)+d/2;%???
lambda_s=(h_s*(1+(3*k_p)))/(12*w_s);
lambda_d=(5*g/w_s)/(5+(4*g/w_s));
lambda_e=0.3*((3*k_p)-1);
bracket=(lambda_s+lambda_d)*(W_s/q);
X1s=(8*mu0*pi*f/p)*(bracket)*(N_1^2);
X1e=(8*mu0*pi*f/p)*(lambda_e*l_ce)*(N_1^2);
X1=X1s+X1e;
X1_check=(8*mu0*pi*f/p)*(bracket+(lambda_e*l_ce))*(N_1^2);
%%%Magnetising reactuance, Xm
alpha=pi/(m*q);
k_d=(sin(q*alpha/2))/(q*(sin(alpha/2)));
constant=w_s/(2*g);
gamma=(4/pi)*((constant*atan(constant))-log(sqrt(1+(constant^2))));
k_c=lambda/(lambda-(gamma*g));
g_e=k_c*g;
W_se=W_s+g_e;
X_m=((24*mu0*f*W_se)/(pi*p*g_e))*((K_w*N_1)^2)*tau;
%%%Secondary resistance, R2
G=(2*mu0*f*(tau^2))/((rho_r/d)*pi*g_e);
R2=X_m/G;
%%%%%%%%%%%%%%%%%%%%% STEADY-STATE CHARACTERISTICS %%%%%%%%%%%%%%%%%%%%%%%%
%%%Total Impedance
v_s=2*f*tau;
s=1-(v_r/v_s);
Z=R1+(X1*j)+((j*X_m*(R2/s))/((j*X_m)+(R2/s)));
Z_real=real(Z);
Z_imag=imag(Z);
Z_mag=sqrt((Z_real^2)+(Z_imag^2));
phi=atan(Z_imag/Z_real);
%%%Primary phase voltage, V1
V1=I1*Z_mag;
%%%Secondary phase current, I2
I2=(X_m/(sqrt(((R2/s)^2)+(X_m^2))))*I1;
%%%Magnetising current, Im
Im=((R2/s)/(sqrt(((R2/s)^2)+(X_m^2))))*I1;
%%%Input active power, P_input
P_input=m*V1*I1*cos(phi);
%%%Output power, P_output
P_output=P_input-(m*(I1^2)*R1)-(m*(I2^2)*R2);
%%%Electromagnetic Force, F
F(end+1) = P_output/v_r
F_check=(3*(I1^2)*R2)/(2*s*f*tau*(((1/(s*G))^2)+1));
%%%Magnetomotive Force, MMF
MMF=(4*sqrt(2)*m*K_w*N_1*Im)/(pi*p);
%%%Efficiency, eff
eff=P_output/P_input;
end
plot(g_m(2:end),F)