磁矩矩阵matlab-m文件的数据分析处理

利用matlab自带quiver3函数可视化磁矩矩阵m—Arrow001()

function [D] = Arrow001(m)
% m(:,:,1) point right
% m(:,:,2) point down
% m(:,:,3) point up perpedicular to screen
% Crop the data for 3D plot

% position in imageJ coordinate
x1 = 0;
y1 = 0;
% size
dx = 1;
dy = 1;

% SPLEEMcolorplot
% contrast = 0.65;
theta = 0;
HSI = spin2hsi2(m,theta);
imcolorplot = hsi2rgb(HSI);
h4 = figure(4);
hold on
set(h4,'outerposition',[300 300 500 500])
imshow(imcolorplot);
rectangle('Position',[x1+0.5 y1+0.5 dx dy],...
    'EdgeColor',"yellow",'LineWidth',1)
set(gca,'position',[0 0 1 1])
hold off

% crop ROI
dx = dx -1;
dy = dy -1;
x1 = x1 +1;
y1 = y1 +1;
x2 = x1 + dx;
y2 = y1 + dy;
alpha = m(:,:,1);
beta  = m(:,:,2);
gamma = m(:,:,3);
U = alpha(y1:y2,x1:x2);
V =  beta(y1:y2,x1:x2);
W = gamma(y1:y2,x1:x2);
% Generate mesh coordinate
N = size(m(:,:,1));
x = 1:1:N(1);
y = 1:1:N(2);
[meshX,meshY]= meshgrid(y,x);
meshZ = zeros(N(1),N(2));
X = meshX(y1:y2,x1:x2);
Y = meshY(y1:y2,x1:x2);
Z = meshZ(y1:y2,x1:x2);

% shift plot center to O point
X_m = ceil((max(max(X)) + min(min(X)))/2);
Y_m = ceil((max(max(Y)) + min(min(Y)))/2);
Z_m = ceil((max(max(Z)) + min(min(Z)))/2);
D.X = X - X_m;
D.Y = Y - Y_m;
D.Z = Z - Z_m;
% D.X = X;
% D.Y = Y;
% D.Z = Z;
D.U = U;
D.V = V;
D.W = W;

% Plot
scale = .8;
figure(1)
quiver3(X,Y,Z,U,V,W,scale,'LineWidth',2)
set(gca,'ydir','reverse')
%axis equal off
%daspect([1 1 3])

end

function [hsi] = spin2hsi2(spin,theta)
% spin space to hsi space
alpha = spin(:,:,1);
beta  = spin(:,:,2);
gamma = spin(:,:,3);
H = atan2(-alpha,-beta) + pi/2;
H(H<0) = H(H<0) + 2*pi;
H = H + theta/180*pi;
H(H<0) = H(H<0) + 2*pi;
H(H>2*pi) = H(H>2*pi) - 2*pi;
% S = ones(size(H));
% I = ones(size(H))/2;
S = sqrt(alpha.^2 + beta.^2);
I = (gamma+1)/2;
m = alpha.^2 + beta.^2 + gamma.^2;
H(m==0) = 0;
S(m==0) = 0;
I(m==0) = 1;
hsi = cat(3,H,S,I);
end

function [rgb] = hsi2rgb(hsi)
% hsi space to rgb space
H = hsi(:,:,1);
S = hsi(:,:,2);
I = hsi(:,:,3);
R = zeros(size(hsi,1),size(hsi,2));
G = zeros(size(hsi,1),size(hsi,2));
B = zeros(size(hsi,1),size(hsi,2));
idx = find((0 <= H)&(H < 2*pi/3));
B(idx) = I(idx).*(1-S(idx));
R(idx) = I(idx).*(1+S(idx).*cos(H(idx))./cos(pi/3-H(idx)));
G(idx) = 3*I(idx)-(B(idx)+R(idx));
idx = find((2*pi/3 <= H)&(H < 4*pi/3));
R(idx) = I(idx).*(1-S(idx));
G(idx) = I(idx).*(1+S(idx).*cos(H(idx)-2*pi/3)...
    ./cos(pi-H(idx)));
B(idx) = 3*I(idx)-(G(idx)+R(idx));
idx = find((4*pi/3 <= H)&(H <= 2*pi));
G(idx) = I(idx).*(1-S(idx));
B(idx) = I(idx).*(1+S(idx).*cos(H(idx)-4*pi/3)...
    ./cos(5*pi/3-H(idx)));
R(idx) = 3*I(idx)-(G(idx)+B(idx));
rgb = cat(3,R,G,B);
rgb = max(min(rgb,1),0);
end

function [colorwheelRGBImage] = Colorwheel(r,azimuth)
x = -r:r;
x = x/r;
y = x;
[X,Y] = meshgrid(x,y);
phi = atan2(Y,X);
phi(phi<0) = phi(phi<0)+2*pi;
phi = phi + azimuth/180*pi;
phi(phi<0) = phi(phi<0) + 2*pi;
phi(phi>2*pi) = phi(phi>2*pi) - 2*pi;
hue = phi/2/pi;
saturation = sqrt(X.^2+Y.^2);
saturation(saturation>1) = 0;
value = ones(size(X));
imHSV = cat(3,hue,saturation,value);
colorwheelRGBImage = hsv2rgb(imHSV);
colorwheelRGBImage = flipud(colorwheelRGBImage);
end

function [] = Colorwheel2(theta)
x0 = 0;
y0 = 0;
t = (0:.02:2)*pi;
r1 = 0:.02:.3;   % interior white
r2 = .3:.02:.8;  % colorwheel
r3 = .8:.02:1;   % exterior black
xp1 = x0+cos(t)'*r1;
yp1 = y0+sin(t)'*r1;
xp2 = x0+cos(t)'*r2;
yp2 = y0+sin(t)'*r2;
xp3 = x0+cos(t)'*r3;
yp3 = y0+sin(t)'*r3;
surface(xp2,yp2,t'*(r2 == r2),'EdgeColor','none');
surface(xp3,yp3,t'*(r3 == r3),'FaceColor',[0 0 0]);
surface(xp1,yp1,t'*(r1 == r1),'FaceColor',[1 1 1],'EdgeColor','none');
colormap(hsv(256))
axis equal off
view(theta,90)
end

function [] = Colorwheel3(azimuth)
r = linspace(0,1,100);
theta = linspace(0,2*pi,1000);
[rg,thg] = meshgrid(r,theta);
[x,y] = pol2cart(thg,rg);
pcolor(x,y,thg);
colormap(hsv)
view(azimuth,90)
shading flat
axis equal off
end

利用自建quiver3_Rendering函数可视化磁矩矩阵m—Arrow002()

function [] = Arrow002(D)
% 3D colorplot for SPLEEM/MC simulation/Mumax simulation results

close all
h1 = figure(1);
% Rendering selector
% 1 ~ HSL
% 2 ~ RWB
% 3 ~ single color
theta = 0;
quiver3_Rendering(D.X,D.Y,D.Z,D.U,D.V,D.W,theta,1);
set(gca,'ydir','reverse')
axis equal off

% camera settings
camlight('left')
camproj('perspective')
% campos([ -5.7932   97.6476  122.0655])
% camtarget([-0.5 -0.5 0.1])
material shiny

%exportgraphics(h1,'SPLEEM_vectorplot.tiff','Resolution',600);

end

function [] = quiver3_Rendering(X,Y,Z,U,V,W,theta,color)

[N1,N2] = size(X);

for i=1:N1
    for j=1:N2

        a0 = X(i,j);
        b0 = Y(i,j);
        c0 = Z(i,j);
        p = [a0,b0,c0];
        alpha0 = U(i,j);
        beta0  = V(i,j);
        gamma0 = W(i,j);
        d = [alpha0,beta0,gamma0];
        % 0.2,0.65,0.08,0.35
        quiver_Refine(1,p,d,theta,color,0.2,0.65,0.08,0.35,1);
        hold on

    end
end

end

function [s1,s2,s3,s4] = quiver_Refine(s,x,V,theta,color,r1,h1,r2,h2,al)
% quiver_Refine(rotation_center,position,direction,color,Cone_r,Cone_h,...
% Cylinder_r,Cylinder_h,alpha)
% s = 1 rotate along the arrow center
% s = 0 rotate along the cylinder center

% color control the colormap
% 1 ~ HSL
% 2 ~ RWB
% 3 ~ Customized

% read data
a = x(1,1);
b = x(1,2);
c = x(1,3);
alpha = V(1,1);
beta  = V(1,2);
gamma = V(1,3);

r1 = 0:0.01:r1;
% Cone
% The cone top position
a1 = a;    
b1 = b;
% Generate cone data
[u,v,w] = cylinder(r1,50);
u = u+a1;
v = v+b1;
w = - w*h1+c+(h1+h2)/2+h1/2-h1/2*s;
% Cone botton
t1 = (0:0.04:2)*pi;
r = max(r1(:));
xc = a1;
yc = b1;
zc = -h1+c+(h1+h2)/2+h1/2-h1/2*s;
% Generate cone botton data
x1 = xc + cos(t1)*r;
y1 = yc + sin(t1)*r;
[m,n] = size(x1);
z1 = repmat(zc,m,n);

% Cylinder
% The cylinder top position
a2 = a;
b2 = b;
c2 = c-h1;
% Generate cylinder data
[x,y,z] = cylinder(r2,50);
x = x+a2;
y = y+b2;
z = -z*h2+c2+(h1+h2)/2+h1/2-h1/2*s;
% Cylinder botton
t2 = (0:0.04:2)*pi;
r = r2;
xc = a2;
yc = b2;
zc = c2-h2+(h1+h2)/2+h1/2-h1/2*s;
% Generate cylinder botton data
x2 = xc + cos(t2)*r;
y2 = yc + sin(t2)*r;
[m,n] = size(x2);
z2 = repmat(zc,m,n);

% Rotate arrow
hold on
origin = [a b c];
theta2 = acos(gamma./sqrt(alpha.^2 + beta.^2 + gamma.^2));
if beta^2 + alpha^2 == 0
    direct = [0 1 0];
else
    direct = [-beta,alpha,0];
end

switch color
    case 1
        % control vector color using HSL
        % Hue:in-plane magnetization
        % Lightness: out-of-plane magnetization
        H = atan2(-beta,alpha); % -beta due to the y direction reverse
        % Get the atan2 space in 0 ~ 2*pi
        if H < 0
            H = 2*pi+H;
        end
        % rotate clockwisely
        H = H + theta/180*pi;
        if H < 0
            H = H + 2*pi;
        elseif H >2*pi
            H = H - 2*pi;
        end
        S = sqrt(alpha.^2 + beta.^2);
        I = (gamma+1)/2;
        % contrast = 1;
        % S = min(1,contrast*(1-gamma));
        % I = min(1,contrast*(1+gamma));
        % change HSL to RGB
        idx = find((0 <= H)&(H < 2*pi/3));
        B(idx) = I(idx).*(1-S(idx));
        R(idx) = I(idx).*(1+S(idx).*cos(H(idx))./cos(pi/3-H(idx)));
        G(idx) = 3*I(idx)-(B(idx)+R(idx));
        idx = find((2*pi/3 <= H)&(H < 4*pi/3));
        R(idx) = I(idx).*(1-S(idx));
        G(idx) = I(idx).*(1+S(idx).*cos(H(idx)-2*pi/3)./cos(pi-H(idx)));
        B(idx) = 3*I(idx)-(G(idx)+R(idx));
        idx = find((4*pi/3 <= H)&(H <= 2*pi));
        G(idx) = I(idx).*(1-S(idx));
        B(idx) = I(idx).*(1+S(idx).*cos(H(idx)-4*pi/3)./cos(5*pi/3-H(idx)));
        R(idx) = 3*I(idx)-(G(idx)+B(idx));
    case 2
        % control vector color using RWB
        % White:in-plane magnetization
        % Red or blue: out-of-plane magnetization
        if gamma < 0
            R = 1;
        else
            R = 1 - gamma;
        end

        if gamma > 0
            B = 1;
        else
            B = 1 + gamma;
        end
        G = 1- abs(gamma);
    case 3
        % control vector color using shades of green color
        R = 0;
        B = 0;
        G = 1;
    otherwise
        R = 0;
        B = 0;
        G = 0;
end

c = [R G B];
c(isnan(c)) = 0;
% Get c in the right section
for i=1:3
    c(i)=max(min(c(i),1),0);
end

%plot
s1 = surf(u,v,w,'Facecolor',c,'Edgecolor','none','FaceAlpha',al);
s2 = surf(x,y,z,'Facecolor',c,'Edgecolor','none','FaceAlpha',al);
s3 = fill3(x1,y1,z1,c,'Edgecolor','none','FaceAlpha',al);
s4 = fill3(x2,y2,z2,c,'Edgecolor','none','FaceAlpha',al);
if theta2 ~= 0
rotate(s1,direct,rad2deg(theta2),origin);
rotate(s2,direct,rad2deg(theta2),origin);
rotate(s3,direct,rad2deg(theta2),origin);
rotate(s4,direct,rad2deg(theta2),origin);
end

% must be enabled to get right axis
axis equal

end

计算磁矩矩阵m中倾角θ的统计图—Theta001()

function [R] = Theta001(m)

% calculate the m_theta distribution
mx = m(:,:,1);
my = m(:,:,2);
mz = m(:,:,3);
% the theta angle between m and normal +z
m_theta = atan2(sqrt(mx.^2 + my.^2),mz)/pi*180;
h1 = figure(1);
hold on
set(h1,'outerposition',[300 300 600 400],'color','w')
xbins = 32;
histogram(m_theta,xbins,'Normalization','probability')
cl=[0.5 0.5 0.5];
plot([45 45],[0 1],'linewidth',1,'color',cl,'linestyle','--')
plot([90 90],[0 1],'linewidth',1,'color',cl,'linestyle','--')
plot([135 135],[0 1],'linewidth',1,'color',cl,'linestyle','--')
xlim([0 180])
ylim([0 0.6])
xticks([0 45 90 135 180])
xlabel('Canting anlge \theta')
ylabel('Probability')
box on
ax = gca;
set(ax,'Fontname','Arial','FontSize',16);
% fold theta
m_theta(m_theta>90) = 180 - m_theta(m_theta>90);
% deleta DW magnetization
% m_theta = filloutliers(m_theta,"previous","mean");
R.theta = mean(m_theta,"all");
R.Q = (1-cos(R.theta/180*pi))/2;
%exportgraphics(h1,'Theta.png','Resolution',600)

end

绘制磁矩矩阵m的Bloch/Neel分布图—Skyrmion001()

function [] = Skyrmion001(m)
% Output Neel/Bloch component distribution
% m(:,:,1) point right
% m(:,:,2) point down
% m(:,:,3) point up
% default theta = 0

% three component plot
figure(1)
imagesc([m(:,:,1) m(:,:,2) m(:,:,3)])
axis equal off
colormap(gray(256))
set(gca,'position',[0 0 1 1])

% plot the Neel/Bloch component
dxTheta = (circshift(m(:,:,3),[-1 0])-circshift(m(:,:,3),[1 0]))/2;
dyTheta = (circshift(m(:,:,3),[0 -1])-circshift(m(:,:,3),[0 1]))/2;
n = atan2(-dxTheta,dyTheta);
t = atan2(dyTheta,dxTheta);
mask = n==0;
phi = atan2(-m(:,:,2),m(:,:,1));
theta = pi/2 - atan2(abs(m(:,:,3)),sqrt(m(:,:,1).^2 + m(:,:,2).^2));
Neel = phi-n;
Bloch = phi-t;
Neel(Neel>pi) = Neel(Neel>pi) - 2*pi;
Neel(Neel<-pi) = Neel(Neel<-pi) + 2*pi;
Bloch(Bloch>pi) = Bloch(Bloch>pi) - 2*pi;
Bloch(Bloch<-pi) = Bloch(Bloch<-pi) + 2*pi;
figure(2)
imagesc([n phi Neel Bloch theta])
axis equal off
colormap(gray(256))
set(gca,'position',[0 0 1 1])

% Neel component colorplot
h3 = figure(3);
set(h3,'outerposition',[300 300 400 400])
Neel = pi - abs(Neel);
Neel = (Neel - pi/2)/pi*2; % Left-handed ~ 1 / Right-handed ~ -1
Neel = Neel.*theta/pi*2; % multiply the in-plane proportion
%Neel(mask) = 0; % set domain to white
imagesc(Neel)
axis equal off
% map = flipud(othercolor('RdBu11'));
map = RWB2;
colormap(map)
clim([-1 1])
%colorbar('Ticks',[],'TickLabels',{})
%set(gca,'position',[0 0 1 1])
Neel = double(Neel);
c = [-1 1];
cmax = max(c);
cmin = min(c);
Neel(Neel>cmax) = cmax;
Neel(Neel<cmin) = cmin;
Neel = double(Neel);
Neel_RGB = 255*(Neel - cmin)/(cmax - cmin);
imwrite(Neel_RGB,map,"NeelComponent.tiff");
%exportgraphics(h3,'colorbar.tiff','Resolution',600)

% Bloch component colorplot
h4 = figure(4);
set(h4,'outerposition',[300 300 400 400])
Bloch = pi - abs(Bloch); % -pi to pi
Bloch = (Bloch-pi/2)/pi*2; % Left-handed ~ 1 / Right-handed ~ -1
Bloch = Bloch.*theta/pi*2; % multiply the in-plane proportion
%Neel(mask) = 0; % set domain to white
imagesc(Bloch)
axis equal off
% map2 = flipud(othercolor('PuOr11'));
map2 = RWB2;
colormap(map2)
clim([-1 1])
%colorbar('Ticks',[],'TickLabels',{})
%set(gca,'position',[0 0 1 1])
Bloch = double(Bloch);
c = [-1 1];
cmax = max(c);
cmin = min(c);
Bloch(Bloch>cmax) = cmax;
Bloch(Bloch<cmin) = cmin;
Bloch = double(Bloch);
Bloch_RGB = 255*(Bloch - cmin)/(cmax - cmin);
imwrite(Bloch_RGB,map2,"BlochComponent.tiff");
%exportgraphics(h4,'colorbar.tiff','Resolution',600)

% % colorplot
% contrast = 0.65;
% theta = 0;
% HSI = spin2hsi(m,theta,contrast);
% %HSI = spin2hsi2(m,theta);
% h4 = figure(4);
% set(h4,'outerposition',[300 300 400 400],'color','w')
% colorplot.RGB = hsi2rgb(HSI);
% imshow(colorplot.RGB);
% set(gca,'position',[0 0 1 1])

% % HSL colorwheel plot
% h3 = figure(3);
% set(h3,'outerposition',[300 300 300 300],'color','w')
% x1 = 0;
% y1 = 0;
% t = (0:.02:2)*pi;
% r1 = 0:.02:0.3; 
% r2 = .3:.02:0.8;
% r3 = 0.8:.02:1;
% r1 = r1*1;
% xp1 = x1+cos(t)'*r1;
% yp1 = y1+sin(t)'*r1;
% xp2 = x1+cos(t)'*r2;
% yp2 = y1+sin(t)'*r2;
% xp3 = x1+cos(t)'*r3;
% yp3 = y1+sin(t)'*r3;
% hold on
% %pcolor(x,y,t'*(r1==r1));
% surface(xp1,yp1,t'*(r1 == r1),'FaceAlpha',0);
% surface(xp2,yp2,t'*(r2 == r2));
% colormap(hsv(256))
% shading interp
% surface(xp3,yp3,t'*(r3 == r3),'FaceColor',[0 0 0]);
% axis equal off
% view(theta,90)
% set(gca,'position',[0 0 1 1])

end

function [hsi] = spin2hsi(spin,theta,contrast)
% spin space to hsi space
alpha = spin(:,:,1);
beta  = spin(:,:,2);
gamma = spin(:,:,3);
H = atan2(beta,alpha);
H(H<0) = H(H<0) + 2*pi;
H = H + theta/180*pi;
H(H<0) = H(H<0) + 2*pi;
H(H>2*pi) = H(H>2*pi) - 2*pi;
S = min(1,contrast*(1-gamma));
I = min(1,contrast*(1+gamma));
m = alpha.^2 + beta.^2 + gamma.^2;
H(m==0) = 0;
S(m==0) = 0;
I(m==0) = 1;
hsi = cat(3,H,S,I);
end

function [hsi] = spin2hsi2(spin,theta)
% spin space to hsi space
alpha = spin(:,:,1);
beta  = spin(:,:,2);
gamma = spin(:,:,3);
H = atan2(beta,alpha);
H(H<0) = H(H<0) + 2*pi;
H = H + theta/180*pi;
H(H<0) = H(H<0) + 2*pi;
H(H>2*pi) = H(H>2*pi) - 2*pi;
S = ones(size(H));
I = ones(size(H))*2/3;
m = alpha.^2 + beta.^2 + gamma.^2;
H(m==0) = 0;
S(m==0) = 0;
I(m==0) = 1;
hsi = cat(3,H,S,I);
end

function rgb = hsi2rgb(hsi)
% hsi space to rgb space
H = hsi(:,:,1);
S = hsi(:,:,2);
I = hsi(:,:,3);

R = zeros(size(hsi,1),size(hsi,2));
G = zeros(size(hsi,1),size(hsi,2));
B = zeros(size(hsi,1),size(hsi,2));

idx = find((0 <= H)&(H < 2*pi/3));
B(idx) = I(idx).*(1-S(idx));
R(idx) = I(idx).*(1+S(idx).*cos(H(idx))./cos(pi/3-H(idx)));
G(idx) = 3*I(idx)-(B(idx)+R(idx));

idx = find((2*pi/3 <= H)&(H < 4*pi/3));
R(idx) = I(idx).*(1-S(idx));
G(idx) = I(idx).*(1+S(idx).*cos(H(idx)-2*pi/3)./cos(pi-H(idx)));
B(idx) = 3*I(idx)-(G(idx)+R(idx));

idx = find((4*pi/3 <= H)&(H <= 2*pi));
G(idx) = I(idx).*(1-S(idx));
B(idx) = I(idx).*(1+S(idx).*cos(H(idx)-4*pi/3)./cos(5*pi/3-H(idx)));
R(idx) = 3*I(idx)-(G(idx)+B(idx));

rgb = cat(3,R,G,B);
rgb = max(min(rgb,1),0);
end

计算磁矩矩阵m的拓扑荷密度以及不同类型DMI下的能量密度矩阵，并计算起沿方位角方向的lineprofile—Skyrmion002()

function [R] = Skyrmion002(m)
% Output different DMI energy densities
% Analysis the DMI energy density versus azimuth angle

mx = m(:,:,1);
my = m(:,:,2);
mz = m(:,:,3);

% calculate the topological charge density
% circshift(m,[Y(row) X(col)])
dmdy = (circshift(m,[1 0])-circshift(m,[-1 0]))/2;
dmdx = (circshift(m,[0 1])-circshift(m,[0 -1]))/2;
chargeDensity = dot(m,cross(dmdx,dmdy),3)/(4*pi);
R.charge = sum(sum(chargeDensity));

% calculate the liftshitz invariants
dmxdx = (circshift(mx,[0 1])-circshift(mx,[0 -1]))/2;
dmzdx = (circshift(mz,[0 1])-circshift(mz,[0 -1]))/2;
dmydy = (circshift(my,[1 0])-circshift(my,[-1 0]))/2;
dmzdy = (circshift(mz,[1 0])-circshift(mz,[-1 0]))/2;
dmydx = (circshift(my,[0 1])-circshift(my,[0 -1]))/2;
dmxdy = (circshift(mx,[1 0])-circshift(mx,[-1 0]))/2;
L_x_xz = mx.*dmzdx-mz.*dmxdx;
L_y_yz = my.*dmzdy-mz.*dmydy;
L_x_zy = mz.*dmydx-my.*dmzdx;
L_y_xz = mx.*dmzdy-mz.*dmxdy;

% Left Neel-type DMI energy density
Neel_L = L_x_xz + L_y_yz;
R.Neel_L = sum(sum(Neel_L));

% Right Neel-type DMI energy density
Neel_R = -L_x_xz - L_y_yz;
R.Neel_R = sum(sum(Neel_R));

% Anti_xRN_yLN DMI energy density
% x-direciotn Right-handed Neel
% y-direction Left-handed Neel
Anti_xRN_yLN = -L_x_xz + L_y_yz;
R.Anti_xRN_yLN = sum(sum(Anti_xRN_yLN));

% Anti_xRN_yLN DMI energy density
% x-direciotn Left-handed Neel
% y-direction Right-handed Neel
Anti_xLN_yRN = +L_x_xz - L_y_yz;
R.Anti_xLN_yRN = sum(sum(Anti_xLN_yRN));

% Left Bloch-type DMI energy density
Bloch_L = L_x_zy + L_y_xz;
R.Bloch_L = sum(sum(Bloch_L));

% Right Bloch-type DMI energy density
Bloch_R = -L_x_zy - L_y_xz;
R.Bloch_R = sum(sum(Bloch_R));

% Anti_xRB_yLB DMI energy density
% x-direciotn Right-handed Bloch
% y-direction Left-handed Bloch
Anti_xRB_yLB = -L_x_zy + L_y_xz;
R.Anti_xRB_yLB = sum(sum(Anti_xRB_yLB));

% Anti_xLB_yRB DMI energy density
% x-direciotn Right-handed Bloch
% y-direction Left-handed Bloch
Anti_xLB_yRB = +L_x_zy - L_y_xz;
R.Anti_xLB_yRB = sum(sum(Anti_xLB_yRB));

% imagesc(Bloch_R)
% colorbar

% Plot
h1 = figure(1);
imagesc([Neel_L Neel_R Anti_xRN_yLN Anti_xLN_yRN; ...
        Bloch_L Bloch_R Anti_xRB_yLB Anti_xLB_yRB])
axis equal off
colormap(RWB)
cmin = -0.2;
cmax = +0.2;
clim([cmin cmax])
% colorbar
% exportgraphics(h1,'DMIEdens.tiff','Resolution',600)
% set(gca,'position',[0 0 1 1])
% output
Neel_L(Neel_L>cmax) = cmax;
Neel_L(Neel_L<cmin) = cmin;
Neel_L = double(Neel_L);
Neel_L = 255*(Neel_L - cmin)/(cmax - cmin);

Neel_R_8bit = Neel_R;
Neel_R_8bit(Neel_R_8bit>cmax) = cmax;
Neel_R_8bit(Neel_R_8bit<cmin) = cmin;
Neel_R_8bit = double(Neel_R_8bit);
Neel_R_8bit = 255*(Neel_R_8bit - cmin)/(cmax - cmin);

Anti_xRN_yLN(Anti_xRN_yLN>cmax) = cmax;
Anti_xRN_yLN(Anti_xRN_yLN<cmin) = cmin;
Anti_xRN_yLN = double(Anti_xRN_yLN);
Anti_xRN_yLN = 255*(Anti_xRN_yLN - cmin)/(cmax - cmin);

Anti_xLN_yRN(Anti_xLN_yRN>cmax) = cmax;
Anti_xLN_yRN(Anti_xLN_yRN<cmin) = cmin;
Anti_xLN_yRN = double(Anti_xLN_yRN);
Anti_xLN_yRN = 255*(Anti_xLN_yRN - cmin)/(cmax - cmin);

Bloch_L(Bloch_L>cmax) = cmax;
Bloch_L(Bloch_L<cmin) = cmin;
Bloch_L = double(Bloch_L);
Bloch_L = 255*(Bloch_L - cmin)/(cmax - cmin);

Bloch_R_8bit = Bloch_R;
Bloch_R_8bit(Bloch_R_8bit>cmax) = cmax;
Bloch_R_8bit(Bloch_R_8bit<cmin) = cmin;
Bloch_R_8bit = double(Bloch_R_8bit);
Bloch_R_8bit = 255*(Bloch_R_8bit - cmin)/(cmax - cmin);

Anti_xRB_yLB_8bit = Anti_xRB_yLB;
Anti_xRB_yLB_8bit(Anti_xRB_yLB_8bit>cmax) = cmax;
Anti_xRB_yLB_8bit(Anti_xRB_yLB_8bit<cmin) = cmin;
Anti_xRB_yLB_8bit = double(Anti_xRB_yLB_8bit);
Anti_xRB_yLB_8bit = 255*(Anti_xRB_yLB_8bit - cmin)/(cmax - cmin);

Anti_xLB_yRB(Anti_xLB_yRB>cmax) = cmax;
Anti_xLB_yRB(Anti_xLB_yRB<cmin) = cmin;
Anti_xLB_yRB = double(Anti_xLB_yRB);
Anti_xLB_yRB = 255*(Anti_xLB_yRB - cmin)/(cmax - cmin);

imwrite(Neel_L,RWB,"Neel_L.tiff");
imwrite(Neel_R_8bit,RWB,"Neel_R.tiff"); %
imwrite(Anti_xRN_yLN,RWB,"Anti_xRN_yLN.tiff");
imwrite(Anti_xLN_yRN,RWB,"Anti_xLN_yRN.tiff");
imwrite(Bloch_L,RWB,"Bloch_L.tiff");
imwrite(Bloch_R_8bit,RWB,"Bloch_R.tiff");%
imwrite(Anti_xRB_yLB_8bit,RWB,"Anti_xRB_yLB.tiff");%
imwrite(Anti_xLB_yRB,RWB,"Anti_xLB_yRB.tiff");

% circle profile of DMI energy density
radius = 9;
Xc = 32.5;
Yc = 32.5;
data_Bloch_R = circleprofile(Bloch_R,Xc,Yc,radius);
data_Neel_R = circleprofile(Neel_R,Xc,Yc,radius);
data_Anti_xRB_yLB = circleprofile(Anti_xRB_yLB,Xc,Yc,radius);

f2 = figure(2);
hold on
box on
L = 3;
c_Bloch = [0.4940 0.1840 0.5560];
c_Neel = [0 0.4470 0.7410];
c_Anti = [0.8500 0.3250 0.0980];
plot(data_Bloch_R(:,1),data_Bloch_R(:,2),"LineWidth",L,"Color",c_Bloch)
plot(data_Neel_R(:,1),data_Neel_R(:,2),"LineWidth",L,"Color",c_Neel)
plot(data_Anti_xRB_yLB(:,1),data_Anti_xRB_yLB(:,2),"LineWidth",L,"Color",c_Anti)
xlabel("Azimuth angle \theta (°)")
ylabel("DMI energy density")
xticks([0 120 240 360])
yticks([])
xlim([0 360])
% ylim([-0.2 0.2]) %Bloch,Neel,Anti
ylim([-0.3 0.3])
legend("Bloch","Neel","Anti")
set(gca,'fontname','Arial','fontsize',21)
set(f2,'outerposition',[200 200 400 420])
%exportgraphics(f2,'legend.pdf')

end

function [data] = circleprofile(image,Xc,Yc,R)

image = double(image);

%let R be the radius. 
%let Xc and Yc be the origin. These do not need to be integer
dR = 1; % select linewidth
N = 128;   %number of locations to extract a
theta = linspace(0, 2*pi, N+1);
%theta(end) = [];  %2pi itself is right back where we started
vals = [];
for R_i = R-dR:1:R+dR
    [x, y] = pol2cart(theta, R_i);
    x = x + Xc;
    y = y + Yc;
    vals_i = interp2(image, x, y);
    vals = [vals;vals_i];
end

vals_mean = mean(vals,1);
data = [theta/pi*180;vals_mean]';

end