%spinning_magnet
% Load Excel sheet containing voltage vs time data recorded via a
% datalogger for the EMF induced across the terminals of a 1000 turn
% coil, with a magnet spun inside, mounted on a spindle. As the magnet
% slows down, the rate of change of magnetic flux linked reduces, and hence
% the EMF reduces by Faraday's law. This code analyses the (t,V) data, and
% explores the veracity of Faraday's law.
%
% LAST UPDATED by Andy French Feb 2026.

function spinning_magnet

%Fontsize for graphs
fsize = 18;

%Remove data after stop time (s)
tstop = 5.5;

%Remove data before start time (s), and reset time to start of allowed samples.
tstart = 2.5;

%Time window for rolling DFT (s)
t_DFT = 0.1;

%Max frequency for DFT (Hz)
max_f = 30;

%Time after tstart for estimating frequencies (s)
dt_f_estimate = 2;

%%

%Load (times,voltage) data from Excel. Time in seconds, voltage in volts.
[num,txt,raw] = xlsread('Magnet spinning in a coil.xlsx'); t = num(:,2).'; V = num(:,1).';

%Truncate data
i = (t > tstart) & (t<tstop); t = t(i); t = t-t(1); V = V(i);

%Interpolate and smooth data
tt = linspace(0,t(end),10000); V = interp1( t,V,tt); t = tt; V = smooth( V, 0.5 );

%Determine times of EMF maxima and minima
Vmaxmin = []; tmaxmin = [];
for n=2:length(t)-2
    %Check for local maxima
    if ( V(n)>V(n-1) ) && ( V(n)>V(n+1) )
        tmaxmin = [tmaxmin,t(n)];
        Vmaxmin = [Vmaxmin,V(n)];
        
        %Check for local minima
    elseif ( V(n)<V(n-1) ) && ( V(n)<V(n+1) )
        tmaxmin = [tmaxmin,t(n)];
        Vmaxmin = [Vmaxmin,V(n)];
    end
end

%Remove times beyond dt_f_estimate
i = find(tmaxmin > dt_f_estimate );
tmaxmin(i) = []; Vmaxmin(i) = [];

%Determine times between maxima and minima, assuming one follows from the
%other... and hence estimate frequencies (Hz). Assume between max and min
%this is half a period.
fmm = 0.5./diff(tmaxmin);

%Set times and voltages associated with fmm to be the mean average of the
%t,V values used to calculate fmm
for n=2:length(tmaxmin)
    tmm(n-1) = 0.5*( tmaxmin(n) + tmaxmin(n-1) );
    Vmm(n-1) = 0.5*( abs( Vmaxmin(n) ) + abs( Vmaxmin(n-1) ) );
end

%Perform line of best fit of tmm vs f
[yfit,xfit,r,m,c,dm,dc,yupper,ylower,s] = bestfitc(tmm,fmm);

%Plot frequency vs time and overlay line of best fit
plot(xfit,yfit,'k','linewidth',2); hold on; set(gca,'fontsize',fsize);
plot( tmm, fmm,'r+','markersize',16 );
xlabel('Time (s)'); ylabel('Frequency (Hz)'); grid on; box on;
title(['(f/Hz) = ',num2str(c,3),' - ',num2str(abs(m),3),'( t/s ). PMCC = ',num2str(r)]);
print(gcf,'Magnet in coil f vs t.png','-dpng','-r300'); close(gcf);

%Estimate angle vs time
te = linspace(0, dt_f_estimate, 1000 );
theta = 2*pi*c*te + pi*m*(te.^2) - 2*pi*0.1;
Ve = -sin(theta).*( c + m*te );
Ve = Ve*max(Vmm)/max(Ve);

%Perform line of best fit of f vs Vmm to investigate
%veracity of Faraday's law
[yfit,xfit,r,m,dm,yupper,ylower,s] = bestfit(fmm,Vmm);
xfit = linspace(0,max(fmm),1000);
yfit = m*xfit;

%Plot absolute value of Vmaxmin vs f
plot( xfit,yfit,'k-','linewidth',2 ); hold on;
plot(fmm,Vmm,'r-','linewidth',2); set(gca,'fontsize',fsize);
ylabel('Voltage (V)'); xlabel('Frequency (Hz)'); grid on; box on;
title(['EMF = ',num2str(m,3),'( f/Hz ). PMCC=',num2str(r)]);
print(gcf,'EMF vs f.png','-dpng','-r300'); close(gcf);

%Plot t vs V
plot(t,V,'r','linewidth',2); hold on; set(gca,'fontsize',fsize);
plot(te,Ve,'b-','linewidth',1);
xlabel('Time (s)'); ylabel('Voltage induced (V)'); grid on; box on;
plot(tmaxmin,Vmaxmin,'b*'); legend({'Data','Model'})
title('Magnet spinning in 1000 turn coil of diameter 7.4cm');
print(gcf,'Magnet in coil V vs t.png','-dpng','-r300'); close(gcf);

%Compute spectrograph
[fmax,tt,ff,DFT] = spectrograph(t,V,t_DFT,max_f);

%Plot spectrograph
pcolor( tt,ff,abs(DFT) ); set(gca,'fontsize',fsize);
colormap('jet'); xlabel('Time /s'); ylabel('Frequency /Hz'); colorbar('fontsize',fsize);
hold on; interp_colormap(1000); shading interp; box on;
plot(t,fmax,'w','linewidth',2); plot(tmm,fmm,'y--','linewidth',2);
title('Spectrograph of magnet spinning in a coil');
print(gcf,'Magnet in coil spectrograph.png','-dpng','-r300'); close(gcf);

%%

%Determine spectrogaph of (t,V) data, using a rolling Discrete Fourier
%Transform (DFT). fmax is the frequency corresponding to the largest
%Fourier component of the signal. DFT is the output of the DFT, with tt
%(time in s) and ff (frequency in HZ) similarly sized arrays.
%It is assumed (t,V) are produced with at constant time intervals, i.e. a
%fixed sample rate. t_DFT is the time (in s) of the samples used to form
%the DFT.
function [fmax,tt,ff,DFT] = spectrograph(t,V,t_DFT,max_f)

%Number of frequencies, from 0 to max_f (in Hz)
Nf = 1000;

%

%Sample rate for data (Hz)
fs = 1/( t(2) - t(1));

%Number of samples into rolling DFT
N = floor( t_DFT*fs );

%Determine time array (s) which corresponds to DFT
M = floor( length(t)/N ); tt = linspace(0,t(end),M); tt = repmat( tt,[Nf,1] );

%Determine frequency array (Hz) which corresponds to DFT
f = linspace(0,max_f,Nf); ff = repmat(f.',[1,M]);

%Truncate V so length is N*M
V = V(1:N*M);

%Transform V into a matrix of N rows and M columns
V = reshape(V,[N,M]);

%Determine DFT samples weighting function
w = dolphcheb(N,40); w = repmat( w, [Nf,1]) ;
w = ones(Nf,N);

%Determine rolling DFT
fff = repmat( f.',[1,N] );
DFT_index = repmat( 0:(N-1),[Nf,1] );
for m=1:M
    x = repmat( V(:,m).', [Nf,1] );
    y = sum( w.*x.*exp( -2*pi*1i*( DFT_index.*fff/fs )),2 )/N;
    DFT(:,m) = y(:);
    
    %Determine frequency of largest Fourier component
    imax = find( abs(y) == max(abs(y)) ); imax = imax(1);
    fmax(m) = f(imax);
end

%Interpolate fmax a t values
fmax = interp1( linspace(0,t(end),M), fmax, t );

%%

%DOLPHCHEB
% Function that computes Dolph-Chebyshev windowing function. The main role
% of the window (as with all windows) is to damp out the effects of the
% Gibbs phenomenon that results from truncation of an infinite series (e.g.
% a Discrete Fourier Transform). Note this window is complex and can
% therefore only really be visualised when the amplitude of the window DFT
% is plotted (in dB).
% N      - Length of window sl_level. This ifference between peak and
%          sidelobes /dB in the DFT of the window.
% w      - Window 'weights'. These are real, positive numbers.
function w=dolphcheb(N,sl_level)
gamma=sl_level/20;
if N>1
    M=N-1;
    alpha=cosh(acosh(10^gamma)/M);
    for n=1:M
        A(n)=((-1)^(n-1))*cos(M*acos(alpha*cos(pi*(n-1)/M)));
    end
    B=fft(A); B(1)=0.5*B(1); B(N)=B(1);
    w=B/max(B);
else
    w=1;
end

%%

%smooth
% Smoothing function.
% x is the input vector of values to smooth
% K is the smoothing constant 0 < K < 1. K=1 implies total smoothing, K=0
% implies no smoothing at all.
function y = smooth( x, K )
y=zeros(size(x));
Y=x(1);
for n=1:length(x)
    y(n) = (1-K) * x(n) + K * Y;
    Y = y(n);
end

%%

%dB
% Function which converts a quantity to decibels i.e. a logarithmic scale.
% The maximum value is taken to be 0 dB.
function y = dB(x)
y = NaN(size(x));
x = abs(x);
maxx = max(max(x));
if maxx > 0
    y(x>0) = 10*log10( x(x>0) ) - 10*log10( maxx );
end

%%

%interp_colormap
% Function which interpolates current colourmap to yield better graduated
% shading. N is number of possible colours.
function interp_colormap(N)

%Get current colourmap, and Initialise new colormap
map = colormap;
new_map = ones(N,3);

%Get size of current colormap and initalise red,green,blue vectors
dim = size(map);
R = ones(1,dim(1)); G = ones(1,dim(1)); B = ones(1,dim(1));
RR = ones(1,N); GG = ones(1,N); BB = ones(1,N);

%Populate these with current colormap
R(:) = map(:,1); G(:) = map(:,2); B(:) = map(:,3);

%Interpolate to yield new colour map
x = linspace( 1, dim(1), N );
RR = interp1( 1:dim(1), R, x ); GG = interp1( 1:dim(1), G, x );
BB = interp1( 1:dim(1), B, x );
new_map(:,1) = RR(:); new_map(:,2) = GG(:); new_map(:,3) = BB(:);

%Set colormap to be new map
colormap( new_map );

%%

%Line of best fit function yfit = m*x +c, with product moment correlation
%coefficient r
function [yfit,xfit,r,m,c,dm,dc,yupper,ylower,s] = bestfitc(x,y)

%Find any x or y values that are NaN or Inf
ignore = isnan(abs(x)) | isnan(abs(y)) | isinf(abs(x)) | isinf(abs(y));
x(ignore) = []; y(ignore) = [];
n = length(x);

%Compute line of best fit
xbar = mean(x); ybar = mean(y); xybar = mean(x.*y);
xxbar = mean(x.^2 ); yybar = mean(y.^2 );
Vx = xxbar - xbar^2; Vy = yybar - ybar^2;
COVxy = xybar - xbar*ybar;
m = COVxy/Vx; c = ybar - m*xbar; r = COVxy/sqrt( Vx*Vy );
[x,i] = sort(x); y = y(i);
yfit = m*x + c; xfit = x;

%Compute errors in gradient m and intercept c
s = sqrt( (1/(n-2))*sum( (y - yfit).^2 ) );
dm = s/sqrt(n*Vx);
dc = ( s/sqrt(n) )*sqrt( 1 + (xbar^2)/Vx );

%Determine envelope lines
yupper = (m+dm)*x + c - dc;
ylower = (m-dm)*x + c + dc;

%%

%Line of best fit function yfit = m*x, with product moment correlation
%coefficient r
function [yfit,xfit,r,m,dm,yupper,ylower,s] = bestfit(x,y)

%Find any x or y values that are NaN or Inf
ignore = isnan(abs(x)) | isnan(abs(y)) | isinf(abs(x)) | isinf(abs(y));
x(ignore) = []; y(ignore) = [];
n = length(x);

%Compute line of best fit
xbar = mean(x); ybar = mean(y); xybar = mean(x.*y);
xxbar = mean(x.^2 ); yybar = mean(y.^2 );
Vx = xxbar - xbar^2; Vy = yybar - ybar^2;
COVxy = xybar - xbar*ybar;
m = xybar/xxbar; r = COVxy/sqrt( Vx*Vy );
[x,i] = sort(x); y = y(i);
yfit = m*x; xfit = x;

%Compute error in gradient m
s = sqrt( (1/(n-1))*sum( (y - yfit).^2 ) );
dm = s/sqrt(n*Vx);

%Determine envelope lines
yupper = (m+dm)*x;
ylower = (m-dm)*x;

%End of code

