function curves=SaderCalcCurves(paramatrix,plotting)

%Program that calculates the thermal motion of a
%rectangular cantilever, consisting of a maximum of two layers. The full
%viscous model is used, for air and water. Only the first bending mode is
%considered.
%If plotting is 1,the motion is plotted in figure 1, for air and water. The values for the curves are output to the
%main program if the function is properly called
%
%%syntax:    input:(paramatrix=[length, width, thickness1, Youngs modulus1, density1, thickness2, Young's modulus2, density2],plotting)
%            output: [freq, W_a_end, W_w_end]
%                   freq,       x-axis vector, contains frequency values 
%                   W_a_end,    Power spectral density of the thermal motion in air, in m/sqrt(Hz)
%                   W_w_end,    Power spectral density of the thermal motion in water, in m/sqrt(Hz)
%
%uses:      SaderAlfaFunc.m, which uses SaderHydroRect.m, they should be
%           somewhere in the paths Matlab looks for its functions


% constants
C1=1.8751041e+000;  % smallest positive root of 1+cos(Cn) cosh(Cn) =0
kBT=300*1.38e-23;   %thermal energy

%material constants for air & fluid
rho_a=1.18;
eta_a=1.86e-5;
rho_w=997;
eta_w=8.59e-4;

%extract properties from input matrix
L=paramatrix(1);W=paramatrix(2);t1=paramatrix(3);t2=paramatrix(6);
t=t1+t2;
I_1=W*(t1^3)/12;I_2=W*(t2^3)/12;
h1=t2/2; h2=t1/2;
E_cI_c=(paramatrix(4))*(I_1+W*t1*h1^2) + paramatrix(7)*(I_2+W*t2*h2^2);
rho_c=(t1*paramatrix(5)+t2*paramatrix(8))/t;

% calculate k and omega_0 in vacuum
omega_0=(C1^2/L^2)*sqrt(E_cI_c/(rho_c*W*t));
k=3*E_cI_c/L^3;

delta=0.0003; %relative to omega_0, is integration step
omega=omega_0*(0.001:delta:2);

%start the real calculations
phi= cos(C1)-cosh(C1)+ ((cos(C1)+cosh(C1))/(sin(C1)+sinh(C1))) *(sinh(C1)-sin(C1)); %mode shape at the end

% air
alfa_a= SaderAlfafunc(rho_a,eta_a,omega,omega_0,rho_c,W,t,C1); %vector of alfa values for each frequency
intalfa_a=delta.*omega_0.*trapz(abs((alfa_a').^2))';  %frequency integral of abs square alfa
FD_a=(3*pi*kBT./k)./(((C1)'.^4).*(intalfa_a)); %eq 27 in the paper. What is this really, F suggests force, but this should not depend on k?
W_a_sq=(FD_a.*(abs(alfa_a).^2))'*phi.^2; %square displacement spectral density 

%water
alfa_w= SaderAlfafunc(rho_w,eta_w,omega,omega_0,rho_c,W,t,C1); %vector of alfa values for each frequency
intalfa_w=delta.*omega_0.*trapz(abs((alfa_w').^2))';  %frequency integral of abs square alfa
FD_w=(3*pi*kBT./k)./(((C1)'.^4).*(intalfa_w)); %eq 27 in the paper. What is this really, F suggests force, but this should not depend on k?
W_w_sq=(FD_w.*(abs(alfa_w).^2))'*phi.^2; %square displacement spectral density 

%finally extract the information you want
freq=(omega/(2*pi))'; %frequency vector
W_a_end=sqrt(W_a_sq);W_w_end=sqrt(W_w_sq);
if plotting==1
figure(1);plotwater=plot(freq,W_w_end);hold on;plotair=plot(freq,W_a_end);hold off
end
curves=[freq, W_a_end, W_w_end];