%%if you want to use this code please cite to this paper ''Mojtaba Ghodsi, H. Ziaiefar, M. Mohammadzaheri, and P. Soltani, 
%%"Effect of Damping on Performance of Magnetostrictive Vibration Energy Harvester", World Academy of Science, Engineering and Technology, 
%%International Journal of Mechanical and Mechatronics Engineering, vol.14, No:4, pp.147-151 (2020).%% 
clc;clear all;close all;
E=72*(10^9); % eleasticity modulus of Galfenol USEd for Harvested Voltage
I=6.1342/E;
IE=E*I; % This IE has been measured and calculted by  FL3/3EI Used for Resonace frequency
%1i= imaginary number=squrt(-1)
mu=(300*4*pi)*10^(-7); %Variable for different materilas 
d=26E-12; % magnetostriction coefficient
re=0.05;  %external damping 
ri=0.05;  %internal damping (Althoougth it has high effect, but we ignore)
N=80; % number of coil turns
L=90E-3; %length of cantilever beam
Lp=30E-3; %active length of MSM
L0=0; % first point of permendur
h=3E-3; % Centroid of permandur to center of beam
A=8.5E-6; % cross section area
z=90E-3;%position of force
m=21.3E-3; % tip mass consist of motor and solid mass
mq=6.55E-7; %Unbalance mass*radiuse
M=14.8E-3% mass of beam
u=M/L;    % mass per unit length, by considering density of Galfenol=7970 kg/m3
FR=300; %Frequency bandwidth
%% Analytic Deflection + Frequency Response
           q=0;
           for f=1:1:FR
                o=2*pi*f;% Driveing frequency
                F0=mq*o^2; %Applied force
           q=q+1;  
           p=0; 
                for x=0:0.0001:L
                   p=p+1;
                   k=(((-u*o^2)+1i*o*re)/(IE+1i*o*ri))^(1/4);
                   B1=(-(sin(k*L)+sinh(k*L)/(cos(k*L)+cosh(k*L))))*(sin(k*L)-sinh(k*L))-cos(k*L)-cosh(k*L);
                   B21=(sin(k*L)-sinh(k*L));
                   B22=(sin(k*L)+sinh(k*L))/(cos(k*L)+cosh(k*L));
                   B23=(cos(k*L)-cosh(k*L));
                   B2=B21-B22*B23;
                   W1=(F0*L^3/IE)/((B1*(k^3)*(L^3)+(B2*m*o^2*L^3)/(IE)));
                   W2=sin(k*x)-sinh(k*x);
                   W3=((sin(k*L)+sinh(k*L))/(cos(k*L)+cosh(k*L)));
                   W4=(cos(k*x)-cosh(k*x));
                   W(1,p)=((-W1)*(W2-W3*W4));
                   X(1,p)=x;
                end
                            
                   MW(1,q)=max (abs(W));
                   F(1,q)=f;
           end 
             
Fres=0;
MMM=max (MW);
      for i=1:1:FR
        if MMM==MW(1,i)  
                 Fres=i;
        end
      end
           
      
      
%% Analytic Voltage
o=2*pi*Fres; % change to rad/sed
% Trapezoidal integration method 
ddX=0;
cc=0;
 for x=L0:0.0000001:Lp
 cc=cc+1;
X1=(F0*L^3/IE)/((B1*(k^3)*(L^3)+(B2*m*o^2*(L^3)/IE)));
X2=-sin(k*x)-sinh(k*x);
X3=((sin(k*L)+sinh(k*L))/(cos(k*L)+cosh(k*L)));
X4=(cos(k*x)+cosh(k*x));
ddX(1,cc)=abs(((-X1)*k^2*(X2+X3*X4)));
 end
 NN=length (ddX);
 section=(Lp-L0)/(NN-1);
 Int= (section/2)*(2*sum(ddX)-ddX(1,1)-ddX(1,NN));

 ppp=0;
 qqq=0;
 
for R=0.001:.1:10
        qqq=qqq+1;
        ppp=0;
     for t=1:0.001:1+(10/FR)
          ppp=ppp+1;
          G1=N^2*A*mu/(Lp*R);
          G2=(1i*o*N*A*d*E*h/Lp)*(Int);
          voltage=real((G2/(1+(1i*o*G1)))*(exp(1i*o*t)-exp(-t/G1)));
          Volt(ppp,1)=(voltage);
          time(ppp,1)=t;
     end
       MVolt(qqq,1)=max(Volt);
       Power(qqq,1)=((MVolt(qqq,1))^2)/R;
       Resi(qqq,1)=R;
end



 %% Plot
figure
plot(F,MW,'blue')
title('Stiff Displacement vs. frequency')
xlabel('Frequency (Hz)')
ylabel('Displacement (m)')
hold on

figure
plot (Resi,MVolt)
xlabel('External Resistance (\Omega)')
ylabel('Voltage (V)')

figure
plot (Resi,Power)
title('Analitical Power vs. Resistanc')
xlabel('External Resistance (\Omega)')
ylabel('Power (W)')