Reduction of repeated lines in code

4 views (last 30 days)
sriram amarnath
sriram amarnath on 27 Jun 2019
Commented: sriram amarnath on 27 Jun 2019
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

Sign in to answer this question.

Accepted Answer

KSSV
KSSV on 27 Jun 2019
Edited: KSSV on 27 Jun 2019
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
  1 Comment
sriram amarnath
sriram amarnath on 27 Jun 2019
Thank you for the answer can you also help me in refining the code under %% 4000 Hz it takes too much time to run it. Can you please optimize it
1a velocity 0-1kHz
clc; clear all;
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
fv=10;
func=zeros(fv,length(m));
for i=1:fv
func(i,:)=(i*pi./Lx).^2 +(m*pi/Lz).^2; %
end
% func1=(1*pi./Lx).^2 +(m.*pi./Lz).^2;% summation function n=1 & m= 1:1000
%%% I have used matrix approach to minimize your code. Your func1, func2,
%%% func3, func4 etc., are row vectors of func matrix in line 21.
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),fv,length(m));
for k=1:length(w)
for i=1:fv
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
%vfunc2= 1j.*w(1,1)'.*((sin(2*pi*x0/Lx).*sin(m.*pi*z0./Lz)).^2./(-w(1,1)'.^2+wnm2.^2.*(1+1j*eta))).*exp(1j.*w(1,1).*t');% velocity function
v=(4/(Lx*Lz))*(F/M)*vfunc; %diffrentiated Eq 98 Notes
%v2=(4/(Lx*Lz))*(F/M).*vfunc2;% diffrentiated Eq 98 Notes
%v1000= sum(v,1);% total mode contributions at 1kHz
vSum=squeeze(sum(v,2)); % squeeze-removes the signleton dimension
%sv1=sum(real(v1000)); % real part of all modes
RlvSum=sum(real(vSum),2);
% sv=[sv1 sv2 sv3 sv4 sv5 sv6 sv7 sv8 sv9 sv10];% vector of all modes
figure('NumberTitle','off','Name','Velocity 0-1kHz')
plot(f,RlvSum);
xlabel('Frequency Hz')
ylabel('Velocity m/s')
% I didn't go further, you can use a similar approach.
%% 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
%% Shape function
combo=[2 2; 2 3; 5 3]; % (n=2,m=2)using eq 90 & 91 notes % (n=2,m=3)using eq 90 & 91 notes % (n=5,m=3)using eq 90 & 91 notes
psi=zeros(3,length(xr),length(zr));
freqs=[112 445 850];
%Vsh=zeros(length(m),3);
for i=1:length(combo(:,1))
psi(i,:,:)=sin((combo(i,1)*pi.*X)/Lx).*sin((combo(i,2)*pi.*Z)./Lz);
figure
surf(X,Z,squeeze(psi(i,:,:)),'EdgeColor','interp');
xlabel('Length')
zlabel('Displacement')
ylabel('Width')
colorbar
title(['Shape Function',num2str(freqs(i)), 'Hz'])
Vsh(:,i)=getframe(gcf);% save figure to a video frame
end
close all
%% 1c Between Modes
movie(Vsh,2) % video of 3 resonance modes
Msh=VideoWriter('Resonance modes','Uncompressed AVI');
Msh.FrameRate=1; % frame rate
open(Msh);
writeVideo(Msh,Vsh);% writing video
close(Msh);
%% 4000 Hz
wp4k=2*pi*4000; % driving frequency of 4000Hz
dfunc2=zeros(length(xr),length(zr));
tic
for n=1:fv % index k is equiv to n in act. eqn.
for k=1:length(m) % k equiv to m in eqn. note--> m(k)=k
dfunc2=sin(n*pi*x0/Lx).*sin(k*pi*z0/Lz)...
*sin((n*pi*X)/Lx).*sin(k*pi*Z/Lz)...
.*exp(1j*wp*t)./(-wp4k^2+wnm(n,k).^2+2*1j*eta)+dfunc2; % note--> size of wnm is 10x1000 => its indices should be (k,n)
end
end
toc
ydt4k=(4/(Lx*Lz))*(F/M).*real(dfunc2);% displacement of all modes
figure('NumberTitle','off','Name','Driving Frequency 4000 Hz')
surf(X,Z,ydt4k,'EdgeColor','interp');
xlabel('Length')
zlabel('Displacement')
ylabel('Width')
colorbar

Sign in to comment.

More Answers (0)

Categories

Find more on Read, Write, and Modify Image in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!