function ODEsimSpirographPlot
% makes a pretty plot of many robots being driven by a constantly rotating
% magnetic field.
global freq vel a
set(0,'defaultaxesfontsize',14);
set(0,'defaulttextfontsize',14);
set(0,'DefaultTextInterpreter', 'latex')
format compact
clc
freq = 10; % as the frequency increases, the lines become straighter, but the differentiation decreases)
vel = 0.7; %mm/s (Estimate by Yan Ou) 700 \mum/s
tF = 10; %final time
%initial conditions for cells (x,t,theta)
[y,x] = meshgrid(1:2,1:4);
x = reshape(x,numel(x),1);
y = reshape(y,numel(y),1);
N = numel(x); %number of robots
a = linspace(4,11,N)'; % a values, should be between 3 and 10
th = 0*ones(N,1); % performance is very dependent on initial orientation for large amplitude
options = odeset('RelTol',1e-6,'AbsTol',1e-6 );
[T,Y] = ode45(@simPyriformis,[0 tF],[x;y;th],options);
%%%%%%%%%%%%%%%%%%%%%%% PLOTTING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%draw path, cells
figure(1)
clf
c = distinguishable_colors(N);
for i = 1:N
s = plot(Y(1,i),Y(1,i+N),'x','color',c(i,:));
uistack(s,'bottom');
hold on
p = plot(Y(:,i),Y(:,i+N),'-','color',c(i,:));
uistack(p,'bottom');
text(Y(1,i),Y(1,i+N)-.2,['$a$ = ',num2str(a(i))],'HorizontalAlignment','center','Interpreter','latex')
end
axis equal
title({[num2str(N),' cells simulated for ',num2str(tF),' s'];},'Interpreter','latex');%['\theta_M = ',num2str(Mag),'sin(',num2str(freq),'t)']});
xlabel('$x$ mm','Interpreter','latex')
ylabel('$y$ mm','Interpreter','latex')
for i = 1:N
drawCell(Y(end,i),Y(end,i+N),Y(end,i+2*N),a(i),0.05);
end
%fix the axis and xticks
set(gca, 'xtick',[-5:10],'ytick',[-5:10])
axis([-.5,4.5,.5,2.8])
%saveas(gcf,'../../pictures/pdf/SpiralsFromRotatingMagneticField.pdf')
format_ticks(gca,{'0','1','2','3','4'},{'1','2'},[0:4],[1:2],0,0,0.01);
xlabel({' ';'$x$ mm'},'Interpreter','latex')
ylabel({'$y$ mm';' '},'Interpreter','latex')
set(gcf,'PaperUnits','inches')
set(gcf,'papersize',[7,4])
set(gcf,'paperposition',[-.5,-.5,8,5])
print -dpdf '../../pictures/pdf/SpiralsFromRotatingMagneticField.pdf'
end
function dq = simPyriformis(t,q)
% simulates a differential equation model for the cells
% q = [x_coordinates; y_coordinates; theta_coordinates,
% magnetic_field_orientation]
global freq vel a
dq = zeros(size(q)); % a column vector for the change in state
N = (numel(q))/3; %number of robots
th = q(2*N+1:3*N); %current orientations of the cells
dq(1:N) = vel*cos(th);
dq(N+1:2*N) = vel*sin(th);
angDif = t*freq-th;
dq(2*N+1:3*N) = a.*sin(angDif);
end