clear all;
close all;


%Beam properties

rho=5150;                  %beam density
Ebeam=70e9;                %young's modulus of beam
Lbeam=0.003;               %beam length
width=0.001;               %beam and piezo width
Lf=0.200;                  %length to forcing f(t)
thickbeam=0.0002;          %beam thick
Area=thickbeam*width;      %area
I=width*thickbeam^3/12;    %Inertia


%Piezo peoperties
g31=-9.5e-3;               %piezoelectric coefficient
Epiezo=50e9;               %young's modulus of piezo
thickpiezo=0.00005e-3;     %piezo thick
Lpiezo=0.003;              %piezo length
d31=-320e-12;              %m/volt dielectric constant
Rsource = 330000;          %ohm source resistance
c=(Ebeam*I/(rho*Area))^.5;

%Natural Frequencies

beta = [4.2;7.2;12.2;18.2;23.2];
wn=beta.^2*c;
fn=wn/(2*pi);

%ModeShapes

syms x
for i=1:5,
   ModeShape(i) = cosh(beta(i)*x)-cos(beta(i)*x)-(sinh(beta(i)*Lbeam)...
-sin(beta(i)*Lbeam))/(cosh(beta(i)*Lbeam)+cos(beta(i)*Lbeam))...
*(sinh(beta(i)*x)-sin(beta(i)*x));
end

%Harmonic Forcing Function

w=wn;
zeta=0.03;
wd=wn*(1-zeta^2)^.5;
syms t tt
Force=(0.35)*sin(w(1)*tt)*1/(rho*Area)*subs(ModeShape,x,Lf);
for i=1:5
    q(i)=1/wd(i)*exp(-zeta*wn(i)*t)...
        *int(Force(i)*exp(zeta*wn(i)*tt)*sin(wd(i)*(t-tt)),tt,0,t);
end

%Develop voltage equations

for i=1:5
    before_w3(i)=q(i)*ModeShape(i);
end

w3=sum(before_w3);
K_x_t=vpa(diff(w3,x,2),5);
Kavg=vpa((1/Lpiezo)*int(K_x_t,x,0,Lpiezo),5);
tau = [0:0.01:10];
Moment=Ebeam*I.*Kavg;
psi=(Ebeam*thickbeam)/(Epiezo*thickpiezo);
T=thickbeam/thickpiezo;
Veb_sym=(6*g31*psi*(1+T))/(width*thickpiezo*(1+psi^2*T^2+2*psi*(2+3*T+2*T^2))).*Moment;
Veb=subs(Veb_sym,'t',tau);

plot(Veb,'g');
xlabel('time');
ylabel('voltage');