MATLAB: Reduction of repeated lines in code

loop

Need help making the code more readable by reducing the repeated lines(dfunc lines) in section 1b:
%% 1a velocity 0-1kHz
Lx=0.286; % length of plate(m)
Lz=0.198; % breadth of plate(m)
h=3.25e-3;% thickness of plate(m)(also d)
E=1.4e9 ;%Youngs modulus
M=265*h;%Mass per unit area( rowS surface density)
mu=0.3;%Poisson?s ratio
eta=4e-2;%Damping Factor
x0=0.16;%exc=(0.16,0.083), point of excitation(m)

z0=0.083;%exc=(0.16,0.083), point of excitation(m)
F=1;% Point load assumed to be 1 newton @ (x0,z0)
f=linspace(1000,1,10); % frequency array
w=2*pi.*f;% angular frequency array
t=1; %  time 
m=1:1000;% mode numbers width
n=10; % mode numbers length
func=zeros(n,length(m));% initial values as place holders
for i=1:n 
func(i,:)=(i*pi./Lx).^2 +(m*pi/Lz).^2; % modes for m=1:1000 & n=1:10
end
D=E*h^3/(12*(1-mu^2));% flexural density D eq 86 Notes
wnm=sqrt(D/M)*func; % for loop for wnm is not required as it ia a scalar multiple of func
 
% for all frequencies common vfunc is defined. vfunc is multi-dimensional
% array with dimensions as 10,10,1000. 
vfunc=zeros(length(w),n,length(m));% 3D array initial values as place holders
for k=1:length(w)
for i=1:n
vfunc(k,i,:)= 1j*w(k)*((sin(i*pi*x0/Lx).*sin(m.*pi*z0./Lz)).^2./(-w(k)'.^2 ...
                        +wnm(i,:).^2.*(1+1j*eta))).*exp(1j.*w(k).*t);          % velocity function  
end
end
v=(4/(Lx*Lz))*(F/M)*vfunc; %diffrentiated Eq 98 Notes
vSum=squeeze(sum(v,2)); % squeeze-removes the singleton dimension
RlvSum=sum(real(vSum),2);% vector of all modes at chosen freq array
figure('NumberTitle','off','Name','Velocity 0-1kHz')
plot(f,RlvSum);
xlabel('Frequency Hz')
ylabel('Velocity m/s')
 
 
 
%% 1b Displacement plot 
xr=linspace(0,Lx,1000);% 1000 x co-ordinate points on surface
zr=linspace(0,Lz,1000); % 1000  z co-ordinate points on surface
wp=2*pi*500; % pressure plot @ 500 Hz
[X,Z]=meshgrid(xr,zr);% plane to be plotted on 
dfunc1=(sin(pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(1,:).*(1+1j*eta));% mode n=1 & m=1:1000
dfunc2=(sin(2.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((2.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(2,:).*(1+1j*eta));% mode n=2 & m=1:1000
dfunc3=(sin(3.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((3.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(3,:).*(1+1j*eta));% mode n=3 & m=1:1000
dfunc4=(sin(4.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((4.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(4,:).*(1+1j*eta));% mode n=4 & m=1:1000
dfunc5=(sin(5.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((5.*pi.*X))/Lx).*sin(m.*pi.*Z/Lz).*exp(1j*wp*t)./(-wp^2+wnm(5,:).*(1+1j*eta));% mode n=5 & m=1:1000
dfunc6=(sin(6.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((6.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(6,:).*(1+1j*eta));% mode n=6 & m=1:1000
dfunc7=(sin(7.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((7.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(7,:).*(1+1j*eta));% mode n=7 & m=1:1000
dfunc8=(sin(8.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((8.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(8,:).*(1+1j*eta));% mode n=8 & m=1:1000
dfunc9=(sin(9.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((9.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(9,:).*(1+1j*eta));% mode n=9 & m=1:1000
dfunc10=(sin(10.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((10.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(10,:).*(1+1j*eta));% mode n=10 & m=1:1000
dt1=dfunc1+dfunc2+dfunc3+dfunc4+dfunc5+dfunc6+dfunc7+dfunc8+dfunc9+dfunc10;% total all modes
ydt=(4/(Lx*Lz))*(F/M).*real(dt1);% displacement due to all modes
figure('NumberTitle','off','Name','Displacment @ 500 Hz')
surf(X,Z,ydt,'EdgeColor','interp');
view(-90,90)
xlabel('Length')
zlabel('Displacement')
ylabel('Width')
colorbar

Best Answer

This should do:
dfunc = zeros(1000,1000,10) ; 
for n = 1:10
    dfunc(:,:,i) = (sin(2.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((n.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(n,:).*(1+1j*eta));
end
The whole code:
clc; clear all ;
%% 1a velocity 0-1kHz
Lx=0.286; % length of plate(m)
Lz=0.198; % breadth of plate(m)
h=3.25e-3;% thickness of plate(m)(also d)
E=1.4e9 ;%Youngs modulus
M=265*h;%Mass per unit area( rowS surface density)
mu=0.3;%Poisson?s ratio
eta=4e-2;%Damping Factor
x0=0.16;%exc=(0.16,0.083), point of excitation(m)

z0=0.083;%exc=(0.16,0.083), point of excitation(m)
F=1;% Point load assumed to be 1 newton @ (x0,z0)
f=linspace(1000,1,10); % frequency array
w=2*pi.*f;% angular frequency array
t=1; %  time 
m=1:1000;% mode numbers width
n=10; % mode numbers length
func=zeros(n,length(m));% initial values as place holders
for i=1:n 
func(i,:)=(i*pi./Lx).^2 +(m*pi/Lz).^2; % modes for m=1:1000 & n=1:10
end
D=E*h^3/(12*(1-mu^2));% flexural density D eq 86 Notes
wnm=sqrt(D/M)*func; % for loop for wnm is not required as it ia a scalar multiple of func
 
% for all frequencies common vfunc is defined. vfunc is multi-dimensional
% array with dimensions as 10,10,1000. 
vfunc=zeros(length(w),n,length(m));% 3D array initial values as place holders
for k=1:length(w)
for i=1:n
vfunc(k,i,:)= 1j*w(k)*((sin(i*pi*x0/Lx).*sin(m.*pi*z0./Lz)).^2./(-w(k)'.^2 ...
                        +wnm(i,:).^2.*(1+1j*eta))).*exp(1j.*w(k).*t);          % velocity function  
end
end
v=(4/(Lx*Lz))*(F/M)*vfunc; %diffrentiated Eq 98 Notes
vSum=squeeze(sum(v,2)); % squeeze-removes the singleton dimension
RlvSum=sum(real(vSum),2);% vector of all modes at chosen freq array
figure('NumberTitle','off','Name','Velocity 0-1kHz')
plot(f,RlvSum);
xlabel('Frequency Hz')
ylabel('Velocity m/s')
 
 
 
%% 1b Displacement plot 
xr=linspace(0,Lx,1000);% 1000 x co-ordinate points on surface
zr=linspace(0,Lz,1000); % 1000  z co-ordinate points on surface
wp=2*pi*500; % pressure plot @ 500 Hz
[X,Z]=meshgrid(xr,zr);% plane to be plotted on 
dfunc = zeros(1000,1000,10) ; 
for n = 1:10
    dfunc(:,:,i) = (sin(2.*pi*x0/Lx).*sin(m.*pi.*z0/Lz).*sin((n.*pi.*X)/Lx).*sin(m.*pi.*Z/Lz)).*exp(1j*wp*t)./(-wp^2+wnm(n,:).*(1+1j*eta));
end
dt1 = sum(dfunc,3) ; 
ydt=(4/(Lx*Lz))*(F/M).*real(dt1);% displacement due to all modes
figure('NumberTitle','off','Name','Displacment @ 500 Hz')
surf(X,Z,ydt,'EdgeColor','interp');
view(-90,90)
xlabel('Length')
zlabel('Displacement')
ylabel('Width')
colorbar
Related Question