function [freq,k,omega_0,W_a,dW_a,W_w,dW_w] =SaderCalcCilinder(paramatrix)
%syntax: input:length, diameter, E, rho, nmodes
%        output: osc_prop= f_0_vac, f_0_air, f_0_water, k, Q_air, Q_water
%uses:   cn10.dat, file with first 4roots of 1+cos(Cn).*cosh(Cn)

%paramatrix=[14e-6,130e-9,200e9,3200,3]; %parameters used when testing 

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);
E_cI_c=paramatrix(3).*(pi*W^4/64); % for a cilinder: I= pi/4 R^4
rho_c=paramatrix(4);
nmodes=paramatrix(5);   %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

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

x=1; %positions to calculate for, 0=at chip, 1=end
npositions=length(x);
delta=0.0005; %relative to omega_0, is integration step
omega=omega_0*(delta:delta:20);
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     = SaderAlfafuncCil(rho_a,eta_a,omega,omega_0,rho_c,W,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= SaderAlfafuncCil(rho_w,eta_w,omega,omega_0,rho_c,W,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));

% place code to extract information here
