#!/usr/bin/env octave

pkg load optim
close all

%% Get data
fname = "nomag-50-500.txt";
file = fopen(fname);
% All measurements txt files have the same format
data = textscan(file, '%d %f %f', 'HeaderLines', 1);
fclose(file);
% The 1st column is not interesting, it's just an index
counts = data{3}; times = data{2};

plot(times, counts, "o;Data;");
hold on;

% Compute the initial values with which nonlin_curvefit will try to fit the
% parameters to. The point where the speed was minimal - When the right most
% column (count column) was at it's lowest value, it's where we should have a
% zero speed
[min_c, min_i] = min(counts);
[max_c, max_i] = max(counts);
guessed_values_struct = struct(...
	'max_amplitude', max_c/2, ...
	'tau', 50.0, ...
	'frequency', 1/(2 * abs(times(max_i) - times(min_i))), ...
	'phase', max_c/2 ...
);

%% This is an attempt that tries to use a struct parameter
%%harmonic_param_func = @(p, t) p.max_amplitude * exp(-t/(2 * p.tau)) * cos(p.frequency*t + p.phase)
%%[p, fy, cvg, outp] = nonlin_curvefit( ...
%%	harmonic_param_func, guessed_values_struct, ...
%%	counts, times ...
%%);

%% This is an attempt that tries to use a normal array / table / matrix ?? I'm
%% not sure how these are called The function I want to fit to
harmonic_param_func = @(t, p) (exp(-t * p(2)) .* (p(1)*sin(2*pi*p(3)*t) + p(4)*cos(2*pi*p(3)*t)) + 0.759);
%% Using `;` instead of `,` doesn't work either
guessed_values_arr = [ ...
	guessed_values_struct.max_amplitude; ...
	guessed_values_struct.tau; ...
	guessed_values_struct.frequency; ...
	guessed_values_struct.phase ...
]
[y, fy, cvg, outp] = leasqr( ...
	times, ...
	counts, ...
	guessed_values_arr, ...
	harmonic_param_func, ...
	1e-9 ...
)
plot(times, y, "-;Fit;")
xlabel("Time")
ylabel("Amplitude")
print("data_fit.svg")
figure(2);
plot(times, (counts-y), '-; Error;')
xlabel("Time")
ylabel("Data - Fit")
print("error.svg")