
%% Parameters

% a=Sphere Radius, d=distance of sphere edge from NV plane, M Magnetization

a = 1.8; %in µm 

d = 0; %in µm
M = 5e4; %in A/m

% Orientation of magnetization

phi = 23.4 *pi/180;
theta = 150 *pi/180;

% Rotation of the diamond anvil cristallographic axis compared to real axis

beta = -45 *pi/180;

% Parameters of the axis (usually imported from experimental file)

size_pix = 0.35; %in µm
w=51;
h=51;
x_um = (1:h)*size_pix;
y_um = (1:w)*size_pix;

%Localisation of the sphere

x0=x_um(25);
y0=y_um(25);

%Conversion to µm

x = (y_um-y0)*10^-6;
y = (x_um-x0)*10^-6;
a = a*10^-6;
d = d*10^-6;

%% Calculate the magnetic field

Bx=zeros(w,h);By=zeros(w,h);Bz=zeros(w,h);

for i=1:w
    for j=1:h
        [r_pol,theta_pol] = CarttoPol(x(j),y(i));
        [x_cart,y_cart] = PoltoCart(r_pol, theta_pol + beta);        
        B=SphereField(x_cart,y_cart,theta,phi,a,d,M); %ChampSphereLT donne le champ en mT
        Bx(i,j)=B(1);
        By(i,j)=B(2);
        Bz(i,j)=B(3);        
    end
end

%% Convert real axis to camera axis : do x -> -y, y -> x and same on Bx ; By

for i=1:w
    for j=1:h
         Bxint(j,i) = -By(w-i+1,j);
         Byint(j,i) = Bx(w-i+1,j);
         Bzint(j,i) = Bz(w-i+1,j);
    end
end

clear Bx;clear By;clear Bz;

Bx = Bxint;By = Byint;Bz = Bzint;

%% Display results

%fname vient généralement du nom du fichier qu'on cherche à modéliser,
%sinon il suffit d'en choisir un

Bextremtest = 0.5;
Bxmin = -Bextremtest;
Bxmax = Bextremtest;
Bymin = -Bextremtest;
Bymax = Bextremtest;
Bzmin = -Bextremtest;
Bzmax = Bextremtest;

nomSave1 = [fname '-simu'];

fig = figure('Name',nomSave1,'Position',[300,100,900,400]);

sub1 = subplot(1,3,1); 
imagesc(x_um,y_um,real(Bx));
colormap jet;
title({'Bx (mT)';''});
axis('image');
ax = gca;
ax.XAxisLocation = 'bottom';
ax.TickDir = 'out';
xlabel('x (\mum)');
ylabel('y (\mum)');
caxis([Bxmin,Bxmax]);
cbar = colorbar;
cbar.FontSize = 10.5;

sub2 = subplot(1,3,2);  
imagesc(x_um,y_um,real(By));
title({'By (mT)';''});
axis('image');
ax = gca;
ax.XAxisLocation = 'bottom';
ax.TickDir = 'out';
xlabel('x (\mum)');
ylabel('y (\mum)');
caxis([Bymin,Bymax]);
cbar = colorbar;
cbar.FontSize = 10.5;

sub3 = subplot(1,3,3);  
imagesc(x_um,y_um,real(Bz));
title({'Bz (mT)';''});
axis('image');
ax = gca;
ax.XAxisLocation = 'bottom';
ax.TickDir = 'out';
xlabel('x (\mum)');
ylabel('y (\mum)');
caxis([Bzmin,Bzmax]);
cbar = colorbar;
cbar.FontSize = 10.5;
