close all
clear

%% ---Exercise 1---%%
y = linspace(-4,4,100);
J = 0.5:1:2.5;

figure
title('Exercise 1')
hold on %graphic parameter
grid on %graphic parameter
for j = J 
    B = Brillouin(j,y); %get the value of the Brillouin function
    plot(y,B,'DisplayName',['J = ',num2str(j)]) % plotting the Brillouin function
    myB = myBrillouin(j,y); %get the value of the simplified Brillouin function
    plot(y,myB,'DisplayName',['J = ',num2str(j)]) % plotting the simplified Brillouin function
end
set(gca,'FontSize',14);
legend % add the legend to the plot
xlabel('x'); ylabel('Bj');
ylim([-1,1])

%% ---Exercise 2---%%
T = linspace(2,300,100);
B = [0.0001,20]; 

figure
title('Exercise 2')
hold on %graphic parameter
grid on %graphic parameter
for b = B 
    X = generateData2(b,T); %get the value of the Brillouin function
    plot(T,1./X,'DisplayName',['B = ',num2str(b)]) % plotting the Brillouin function
end
set(gca,'FontSize',14);
legend % add the legend to the plot
xlabel('T [°K]'); ylabel('1/X');
xlim([min(T),max(T)]);

%% ---Exercise 3---%%
T = linspace(2,300,100);

figure
subplot(1,2,1)
title('Exercise 3 - Fitting')
hold on %graphic parameter
grid on %graphic parameter
X = generateData3(T); %get the value of the Brillouin function
ft = fittype('x/(a+b*x)'); 
myFit = fit(T',1./X',ft,'StartPoint',[1e-6,1e-8]);


S = 0.5; %free electron estimation
L = 0;
J = 0.5;
g = 3/2 + (S*(S+1)-L*(L+1))/2/J/(J+1); %g-landé factor
u0 = 4*pi*1e-7; % void permeability
kB = 1.380649e-23; % Boltzmann constant
uB = 9.274009994e-24; % Bohr magneton
n = 6.022140857e23;

disp(['C = ',num2str(myFit.a)])
disp(['C linear = ',num2str(n*u0*g^2*J*(J+1)*uB^2/3/kB)])
disp(['X0 = ',num2str(myFit.b)])
disp(['X0 goal = ',num2str(2e-8)])

hold on
grid on
plot(myFit,T,1./X)
set(gca,'FontSize',14);
xlabel('T [°K]'); ylabel('1/X');
xlim([min(T),max(T)]);

subplot(1,2,2)
plot(T,1./(X-myFit.b),'.b')
title('Exercise 3 - Corrected')
grid on
set(gca,'FontSize',14);
xlabel('T [°K]'); ylabel('1/(X-X0)');
xlim([min(T),max(T)]);

%% ---Your function--- %%
function B = myBrillouin(J,x)
    B = (J+1)/3/J*x;
end
%% ---Given function--- %%
function B = Brillouin(J,x)
    % Brillouin function:
    % input :
    %     J : angular momentum
    %     x : g*uB*B*J/kB/T
    % output : 
    %     B : the value used to determine magnetic susceptibility
    B = (2*J+1)/2/J*coth((2*J+1)/2/J*x)-1/2/J*coth(1/2/J*x);
end

function X=generateData2(B,T)
    % simulation of a susceptibility measurement in a magnetic field B and a
    % temperature T. No noise.
    % input :
    %     B : the magnetic field applied during measurement in Tesla
    %     T : the temperature applied during measurement in Kelvin
    % output : 
    %     X : The resulting magnetic susceptibility per mole
    S = 0.5; %free electron estimation
    L = 0;
    J = 0.5;
    g = 3/2 + (S*(S+1)-L*(L+1))/2/J/(J+1); %g-landé factor
    uB = 9.274009994e-24; % Bohr magneton
    u0 = 4*pi*1e-7; % void permeability
    kB = 1.380649e-23; % Boltzmann constant
    n = 6.022140857e23;
    x = g*uB*J*B/kB./T;
    M = n*g*uB*J*Brillouin(J,x); %we suppose one mole per m3
    X = M./(B/u0-M); %corresponding susceptibility from B = u0*(H+M) and M = X*H
end

function X=generateData3(T)
    % simulation of a susceptibility measurement in a magnetic field B=0.1mT and a
    % temperature T. experimental noise added and default susceptibility.
    % input :
    %     T : the temperature applied during measurement in Kelvin
    % output : 
    %     X : The resulting magnetic susceptibility per mole
    
    B = 1e-1; %Applied field in Tesla
    X0 = 2e-8; %default shift in the susceptibility
    S = 0.5; %free electron estimation
    L = 0;
    J = 0.5;
    g = 3/2 + (S*(S+1)-L*(L+1))/2/J/(J+1); %g-landé factor
    uB = 9.274009994e-24; % Bohr magneton
    u0 = 4*pi*1e-7; % void permeability
    kB = 1.380649e-23; % Boltzmann constant
    n = 6.022140857e23;
    x = g*uB*J*B/kB./T;
    M = n*g*uB*J*Brillouin(J,x); %we suppose one mole per m3
    M = M + 0.05*M.*randn(size(T)); % 5% error
    X = M./(B/u0-M)+X0; %corresponding susceptibility from B = u0*(H+M) and M = X*H
end