Homec4science
Diffusion yearII_reports (master) Browse

rSTICAZZI/labo/A4/phase.m~
851e56311b98master

Last Change
Raw File

phase.m~
View Options

fs = [2.5 2.75 2.9 3 3.2 3.4 3.5 3.6 3.8 4 4.05 4.08 4.09 4.1 4.11 4.12 4.15 4.2 4.25 4.3 4.4 4.5 4.6 4.8 5 5.5 6 6.5 7 7.5];
n = length(fs);
omega = 4.02; % m = 0
lambd = 0.093921; % f = 90, m = 0
files = cell(1, n);
phiset = zeros(n, 1);
Oset = zeros(n, 1);
Amax = zeros(n, 1);
for i=1:n
files{i} = sprintf('rawdata/force%g_f90.dat', fs(i));
[data, N] = load_stable(files{i}, 5);
t = data(:,1);
f = data(:,3);
f = adjust_height(f, 20);
[tp, p] = real_peaks(t, f, 3);
Oset(i) = Omega_factor(tp, 20, 0);
phiset(i) = phi(lambd, Oset(i), omega);
tp = tp(20:end);
p = p(20:end);
amax = max(p);
Amax(i) = amax;
end
%% UTILS %%
Ofit = polyfit(fs', Oset, 1);
Olinfit = Ofit(1) * fs + Ofit(2);
figure
plot(fs, Oset, 'k+', fs, Olinfit)
xlabel('fs')
ylabel(sprintf('Omega'))
grid on
phifit = zeros(n, 1);
for i=1:n
phifit(i) = phi(lambd, Olinfit(i), omega);
end
figure
plot(fs, phiset, 'r*', fs, phifit)
xlabel('fs')
ylabel(sprintf('phi'))
grid on
[Afit, p, o, l] = amax_fit(fs', Amax, omega, lambd);
fpr
figure
plot(fs, Amax, 'b*', fs, Afit(fs))
xlabel('fs')
ylabel(sprintf('Amax'))
grid on
%% UTILS %%
function [Afit, p_, o_, l_] = amax_fit(fs, Amax, omega, lambd)
x0 = [0.5 omega lambd];
fitfun = fittype( @(p, o, l, x) p * ((o^2 - x.^2).^2 + 4 * l^2 * x.^2).^(-0.5) );
[curve, ~] = fit(fs,Amax,fitfun,'StartPoint',x0);
coeff = coeffvalues(curve);
p_ = coeff(1);
o_ = coeff(2);
l_ = coeff(3);
Afit = curve;
end
function p = phi(lambd, Omega, omega)
p = atan(2 * lambd * Omega / (omega^2 - Omega^2));
if p < 0
p = p + pi;
end
end
function o = Omega_factor(tp, first, last)
n = length(tp) - last - (first-1);
diff = zeros(n-1, 1);
T = 0;
for i=1:(n-1)
diff(i) = tp((first-1) + i+1) - tp((first-1) + i);
T = T + diff(i);
end
T = T / (n-1);
o = 2 * pi / T;
end
function [data, N] = load_stable(file, filter)
data = load(file);
N = length(data(:,1));
% filter out zeros
zero = N;
while abs(data(zero,1)) < 1e-14
zero = zero - 1;
end
% filter out initial domain
filt = 1;
while data(filt,2) < filter
filt = filt + 1;
end
data = data(filt:zero,2:end);
N = length(data);
end
function f = adjust_height(F, until)
middle = 0;
a = F(end-until:end);
N = length(a);
for i=1:N
middle = middle + a(i);
end
middle = middle / N;
f = F - middle;
end
function [T, M] = real_peaks(t, f, a)
smoothed = smooth(f);
[~, inds] = findpeaks(smoothed);
s = length(inds);
T = zeros(s, 1);
M = zeros(s, 1);
i = 1;
while i <= s
try
fit = polyfit(t(inds(i)-a : inds(i)+a), f(inds(i)-a : inds(i)+a), 2);
T(i) = - fit(2) / (2*fit(1));
M(i) = fit(3) - fit(2)^2 / (4*fit(1));
i = i + 1;
catch
T(i) = [];
M(i) = [];
inds(i) = [];
s = s - 1;
end
end
end
c4science · Help