% fitDischarge.m
% Fits data for the current during a condenser discharge in a LCR circuit
% J.M.Soler, Feb.2023

% Set previous data
clear all
R = 0.350;     % resistance (Ohm)
DR = 0.005;    % resistance uncertainty (Ohm)

% Get discharge data
dat = load('discharge.dat');
t = dat(:,1);   % times, in s
I = dat(:,2);   % positions, in m
dt = t(2)-t(1);
tmax = t(end);

% Set initial parameter values for nlinfit
T = 0.1e-3;        % approximate oscillation period, in s
w = 2*pi/T;        % approximate angular frequency, in rad/s
I0 = 0.15;         % approximate initial current amplitude, in A
beta0(1) = I0/w;   % approximate initial charge, in C
beta0(2) = w;
beta0(3) = 0.5e-3; % approximate decay time, in s

% Fita data
[beta,~,~,covPar,fitErr2] = nlinfit( t, I, @dischargeIntensity_JMS, beta0 );
Q0  = beta(1);
w   = beta(2);
tau = beta(3);
DQ0  = sqrt(covPar(1,1));
Dw   = sqrt(covPar(2,2));
Dtau = sqrt(covPar(3,3));
fitError = sqrt(fitErr2);

% Plot data and fit
figure(11)
Ifit = dischargeIntensity_JMS(beta,t);
plot(t,I,'o', t,Ifit,'LineWidth',1.5)
xlabel('Time (s)')
ylabel('Current (A)')
set(gca,'FontSize',15)
grid on

% Find L and C
L = 0.5*R*tau;
DL = L*(DR/R+Dtau/tau);
C = 1/L/w^2;
DC = C*(DL/L+2*Dw/w);

% Print results
fprintf('Fit error = %6.1e A\n',fitError)
fprintf(' Q0 = (%6.3f +/- %5.3f)*10^-6 C\n',    Q0/1e-6,  DQ0/1e-6)
fprintf('  w = (%6.3f +/- %5.3f)*10^3 rad/s\n',  w/1e3,    Dw/1e3) 
fprintf('tau = (%6.3f +/- %5.3f)*10^-3 s\n',   tau/1e-3, Dtau/1e-3)
fprintf('  L = (%6.3f +/- %5.3f)*10^-6 H\n',     L/1e-6,   DL/1e-6)
fprintf('  C = (%6.3f +/- %5.3f)*10^-6 F\n',     C/1e-6,   DC/1e-6)
fprintf('  R = (%6.3f +/- %5.3f) Ohm\n',         R/1e0,    DR/1e0)
