function propmat=SaderCalcProps(paramatrix)

%Program that calculates the most important parameters 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.
%
%%syntax:    input:paramatrix= [length, width, thickness1, Youngs modulus 1, density1, thickness2, Young's modulus2, density2]
%            output: [k,f_a,f_w,q_a,q_w,force_noise_a,force_noise_w]
%                   k,          spring constant in N/m
%                   f_a,        resonance frequency in air, in Hz
%                   q_a,        quality factor in air.  
%                   force_noise_a, force noise in air, in N/sqrt(Hz)
%                   subscript _w are the analogs for water          
% constants
%
%uses:      SaderAlfaFunc.m, and through that SaderHydroRect.m, should be
%           somewhere in the paths Matlab looks for its functions

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
force_noise_w=sqrt(W_w_sq(2))*k; %force noise is defined as (zero frequency PSD of motion)*(spring constant)
force_noise_a=sqrt(W_a_sq(2))*k;
[peakair,f_a_index]=max(W_a_sq);
f_a=freq(f_a_index); %resonant frequency in air
[peakwater,f_w_index]=max(W_w_sq);
f_w=freq(f_w_index); %resonant frequency in water
[diffa,q_a_index1]=min(abs(W_a_sq(1:f_a_index-1)-(peakair/2)));
[diffw,q_w_index1]=min(abs(W_w_sq(1:f_w_index-1)-(peakwater/2)));
[diffa,q_a_index2]=min(abs(W_a_sq(f_a_index:end)-(peakair/2)));
[diffw,q_w_index2]=min(abs(W_w_sq(f_w_index:end)-(peakwater/2)));
q_a=f_a_index/(q_a_index2-q_a_index1+f_a_index-1); % q is defined as freq/FWHM, where FWHM is calculated for the square of the motion PSD
q_w=f_w_index/(q_w_index2-q_w_index1+f_w_index-1);
propmat=[k,f_a,f_w,q_a,q_w,force_noise_a,force_noise_w];