function [freq,k,W_a,dW_a,W_w,dW_w] =SaderCalcFull(paramatrix)
%program that calculates the full motion as a function of frequency of a
%rectangular cantilever, consisting of a maximum of two layers, for the
%entire length of the cantilever. The first four flexural bending modes are
%considered, and the full viscous model is used.
%syntax:    input:paramatrix= [length, width, thickness1, E1, rho1, thickness2, E2, rho2]
%           output: freq,    contains x (frequency) values
%                   k,       spring constant 
%                   W_a,     2D variable, thermal motion power spectral density for measurement in air as function of position 
%                            along cantilever and frequency
%                   dW_a,    space derivative of W_a, linearly related to what you measure with a laser deflection method
%                   W_w,dW_w analog to W_a, dW_a for water

%uses:      CN10.dat, file with first 4roots of 1+cos(Cn).*cosh(Cn), should be in the current directory
%           SaderAlfaFunc.m, and through that SaderHydroRect.m, should be
%           somewhere in the paths Matlab looks for its functions

%paramatrix=[60e-6,30e-6,0.125e-6,200e9,3200,50e-9,75e9,19200]; %parameters for a biolever, used when testing 

nmodes=4;   %number of modes to calculate

% constants
load CN10.dat; %file with first 4 roots of 1+cos(Cn).*cosh(Cn)
CN10=CN10(1:nmodes);
C1=CN10(1); % smallest positive root of 1+cos(Cn) cosh(Cn) =0
kBT=300*1.38e-23; %thermal energy

%material constants
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;

x=1; %positions to calculate for, 0=at chip, 1=end
npositions=length(x);
delta=0.001; %relative to omega_0, is integration step
omega=omega_0*(0.005:delta:35);
freq=omega/(2*pi);
nfreq=length(omega);
for poscount=1:npositions
    y(2*poscount-1)=0.999*x(poscount);
    y(2*poscount)=x(poscount);
end

%start the real calculations
    %mode shapes & their derivatives, are the same for all media
for poscount2=1:npositions
    phi_help(:,poscount2)=((cos(CN10')+cosh(CN10'))./(sin(CN10')+sinh(CN10')));    %help function to create the right size
    phi_y_help(:,poscount2)=phi_help(:,poscount2);                                              %of matrices and keep the lines short
    phi_y_help(:,poscount2+npositions)=phi_help(:,poscount2);
end
phi= (cos(CN10'*x)-cosh(CN10'*x))+phi_help.*(sinh(CN10'*x)-sin(CN10'*x));      %matrix contains undamped modes of the beam,
yphi=(cos(CN10'*y)-cosh(CN10'*y))+phi_y_help.*(sinh(CN10'*y)-sin(CN10'*y));    %each row = 1 mode, each column = 1 position
for modecount=1:nmodes
    dyphi_0(modecount,:)= [diff(yphi(modecount,:))./diff(y)];
end

for poscount3=1:npositions
    dphi(:,poscount3)=dyphi_0(:,2*poscount3-1);
end

 % thermal motion for air
alfa_a     = SaderAlfafunc(rho_a,eta_a,omega,omega_0,rho_c,W,t,CN10'); %matrix of alfa values, each row = 1 mode, each column= one frequency
intalfa_a = delta.*omega_0.*trapz(abs((alfa_a').^2))';  %produces a column vector with frequency integral of abs square alfa for each mode
FD_a=(3*pi*kBT./k)./(((CN10)'.^4).*(intalfa_a)); % column vector, contains absolute square of the force spectral density (? should k be in there if it really was?)

for modecount2=1:nmodes
    W_a(modecount2,:,:)=sqrt((FD_a(modecount2,:).*phi(modecount2,:).^2)' * abs(alfa_a(modecount2,:)).^2);
   dW_a(modecount2,:,:)=sqrt((FD_a(modecount2,:).*dphi(modecount2,:).^2)' * abs(alfa_a(modecount2,:)).^2);
end

W_a=squeeze(sum(W_a,1));
dW_a=squeeze(sum(dW_a,1));

    %thermal motion for water
alfa_w= SaderAlfafunc(rho_w,eta_w,omega,omega_0,rho_c,W,t,CN10'); %matrix of alfa values, each row = 1 mode, each column= one frequency
intalfa_w=delta.*omega_0.*trapz(abs((alfa_w').^2))'; %produces a column vector with frequency integral of abs square alfa for each mode
FD_w=(3*pi*kBT./k)./(((CN10)'.^4).*(intalfa_w)); % column vector, contains absolute square of the force spectral density (? should k be in there if it really was?)

for modecount2=1:nmodes
    W_w(modecount2,:,:)=sqrt((FD_w(modecount2,:).*phi(modecount2,:).^2)' * abs(alfa_w(modecount2,:)).^2);
   dW_w(modecount2,:,:)=sqrt((FD_w(modecount2,:).*dphi(modecount2,:).^2)' * abs(alfa_w(modecount2,:)).^2);
end

W_w=squeeze(sum(W_w,1));
dW_w=squeeze(sum(dW_w,1));


