import matlab models

0_matlab_example.R

library("BacArena")
library("R.matlab")

mod <- readMATmod("Actinobacillus_pleuropneumoniae_L20.mat")

arena <- Arena(n=10,m=10)
arena <- addOrg(arena, Bac(mod), 10)
arena <- addDefaultMed(arena, Bac(mod))
sim <- simEnv(arena, time=5)

Actinobacillus_pleuropneumoniae_L20.mat

%% This code gives the current through the lateral S-N-S junction with spin split superconductors
% The code solves the Usadel equation using that no
% current should flow out of the junction on the sides, only into the
% superconductor.
%% Input: Thouless energy normalized to Delta, length of normal part, Dynes parameter, phase difference between superconductors
%% Output: 2D figure of DOS as function of position and energy.%Define parameters

%Energy in units of Delta. Flip energy array because it is better for convergence to start
%at E = 2\Delta.
Earr = flip(linspace(0,2,400)+1e-21i);
DOS = zeros(length(Earr),1000);
for ee = 1:length(Earr)
ee
%Phase difference between electrodes
Psi = 0;
%Thouless energy E_Th = D/L^2, where 2L is the size of SNS part of junction.
%Normalized to \Delta.
Eth = 0.25;
% Size of normal region, normalized to L
L1 = 1;
%Initialization
Delta =1;
%Boundary parameter S/N interface
gammaS = 15;
%For code it is easier to change to normalization by E_Th, so normalize
%energies by this quantity instead.
Delta = Delta/Eth;
E = Earr(ee)/Eth;
omega = -1iE;
gammaS = gammaS/Eth;
%Define Green's function in the electrode
gs0 = E/sqrt(E^2-Delta^2);
fs0 = 1iDelta/sqrt(E^2-Delta^2);
%Spatial grid. The SNS junction is always between -1 and 1 via
%normalization. By increasing the window you can add extra N parts.
x = linspace(-2,2,1000);

% Superconductivity is only induced for specific regions (below the electrodes)
gamma = @(x) (gammaSmyHeaviside(x-L1/2)+gammaSmyHeaviside(-x-L1/2))*myHeaviside(1-abs(x));

% Differential equations to be solved in the three parts, Usadel equation
% in Ricatti parametrization.
dydx = @(x,y) [y(2);1/((1+y(1)y(3)))(2y(3)y(2)^2)+y(1)(2omega+gamma(x)gs0+0.5gamma(x)fs0y(1))-0.5gamma(x)fs0;y(4);1/((1+y(1)y(3)))(2y(1)y(4)^2)+y(3)(2omega+gamma(x)gs0+0.5gamma(x)fs0y(3))-0.5*gamma(x)*fs0];
%Set tolerance of convergence ()
options = bvpset('AbsTol',1e-8);

%Boundary conditions: no current out of system at -1 and 1 (isolated system).
bc = @(yL,yR) [yL(2); yL(4);yR(2); yR(4)];
% For periodic system change to bc = @(yL,yR) [yL(1)-yR(1);yL(2)-yR(2);yL(3)-yR(3);yL(4)-yR(4)];

%Initial guess. For E = 2\Delta the Ricatti parameter is quite small and the
%zero solution works well. Then use the solution of the previously used
%energy as starting point.
if ee==1
solinit = bvpinit(x,[0 0 0 0]);
else
solinit = sol;
end
%Solve the problems in the three regions
sol = bvp4c(dydx,bc,solinit,options);
%Extract solution in array form
quantities = deval(sol,x);
%Calculate DOS and Current from results
DOS(ee,:) = real((1-quantities(1,:).*quantities(3,:))./(1+quantities(1,:).*quantities(3,:)));
end
%Plot DOS as a function of position and energy with the good font for axes
%and labels
imagesc(x,Earr,DOS)
set(gca,'YDir','normal')
xlim([-1 1])
c = colorbar;
c.Label.String = 'DOS';
c.Label.Interpreter = 'latex';
ax = gca;
ax.FontSize = 14;
box on
xlabel('$x/L$','Interpreter','latex')
ylabel('$E/\Delta$','Interpreter','latex')