-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathcalculateForces.m
executable file
·173 lines (159 loc) · 8.5 KB
/
calculateForces.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
function forceArray = calculateForces(distanceMatrixXY,distanceMatrix,rc,...
headings,reversals,segmentLength,v_target,k_l,k_theta,theta_0,phaseOffset,...
sigma_LJ,r_LJcutoff,eps_LJ,LJnodes,LJmode,angleNoise,ri,r_overlap,f_hapt,haptotaxisMode,roamingLogInd,f_align)
% updates object directions according to update rules
% issues/to-do's:
% - mixed periodic boundary conditions can be quite slow
% - calculate forces without loops?
% - refactor individual forces into their own functions
% - if using noise for head angle, better use von mises random number?
% - if changing head orientation based on forces instantaneously, move
% angleNoise generation to updateOscillators function
% short-hand for indexing coordinates
x = 1;
y = 2;
N = size(distanceMatrixXY,1);
M = size(distanceMatrixXY,2);
distanceMatrixXY = permute(distanceMatrixXY,[3 4 5 1 2]); % this will make indexing later on faster without need for squeeze()
% format of distanceMatrixXY is now N by M by [x y] by N by M
forceArray = zeros(N,M,2); % preallocate forces
for objCtr = 1:N
% check if worm is currently reversing
if ~reversals(objCtr)
headInd = 1;
bodyInd = 2:M;
else
headInd = M;
bodyInd = (M-1):-1:1;
end
movState = 1 - 2*reversals(objCtr); % =-1 if worm is reversing, 1 if not
% calculate force contributions
% motile force
Fm = NaN(M,2);
ds = NaN(M,2); %change in node positon along worm
% estimate tangent vectors
if M>1
% head tangent: direction from next node's position
ds(1,:) = distanceMatrixXY(objCtr,2,[x y],objCtr,1);
% calculate direction towards previous node's position
for nodeCtr = 2:M
ds(nodeCtr,:) = distanceMatrixXY(objCtr,nodeCtr,[x y],objCtr,nodeCtr-1);
end
% before estimating tangent vectors, calculate polar unit vectors and
% bending angles
phi = atan2(ds(2:M,y),ds(2:M,x));
e_phi = [-sin(phi) cos(phi)]; % unit vector in direction of phi, size M-1 by 2
l = sqrt(sum(ds(2:M,:).^2,2)); % length between node and prev node, length M-1
bendingAngles = unwrap(atan2(ds(:,y),ds(:,x)));
% body tangents: estimate tangent direction by average of directions towards previous node and from next node
ds(2:M-1,:) = (ds(2:M-1,:) + ds(3:M,:))/2;
% (tail tangent: towards previous node)
end
% head motile force
headAngle = headings(objCtr,headInd) + angleNoise*randn();
Fm(headInd,:) = [cos(headAngle), sin(headAngle)];
% haptotaxis - move to the direction of other worms (should this be added before angular noise?)
if f_hapt~=0 % haptotaxis could be attractive or repulsive
if strcmp(haptotaxisMode,'weighted_additive')%&&~roamingLogInd(objCtr) % only update those worms which are not roaming
Fm(headInd,:) = Fm(headInd,:) ...
+ calculateHaptotaxisWeightedAdditive(distanceMatrixXY(:,:,:,objCtr,headInd),...
distanceMatrix(:,:,objCtr,headInd),objCtr,ri,r_overlap,f_hapt);
elseif strcmp(haptotaxisMode,'weighted')%&&~roamingLogInd(objCtr) % only update those worms which are not roaming
Fm(headInd,:) = Fm(headInd,:) ...
+ calculateHaptotaxisWeighted(distanceMatrixXY(:,:,:,objCtr,headInd),...
distanceMatrix(:,:,objCtr,headInd),objCtr,ri,r_overlap,f_hapt);
elseif strcmp(haptotaxisMode,'constant')
Fm(headInd,:) = Fm(headInd,:) ...
+ calculateHaptotaxis(distanceMatrixXY(:,:,:,objCtr,headInd),...
distanceMatrix(:,:,objCtr,headInd),objCtr,ri,r_overlap,f_hapt);
end
end
% Vicsek-type alignment force - only really used for demonstration
if f_align~=0
Fm(headInd,:) = Fm(headInd,:) ...
+ calculateAlignmentForce(headings,distanceMatrix(:,:,objCtr,headInd),objCtr,headInd,ri,f_align,true);
end
% body motile force
Fm(bodyInd,:) = movState*ds(bodyInd,:);
% fix magnitue of motile force to give target velocity
% % Fm = v_target(objCtr).*Fm./sqrt(sum(Fm.^2,2));
Fm = v_target(objCtr).*bsxfun(@rdivide,Fm,sqrt(sum(Fm.^2,2)));
% length constraint
if k_l>0&&M>1
dl = NaN(M-1,2);
for nodeCtr = 1:M-1
dl(nodeCtr,:) = distanceMatrixXY(objCtr,nodeCtr,[x y],objCtr,nodeCtr+1); % direction to next node
end
% % dl = dl./l.*(l - segmentLength);% normalised for segment length and deviation from segmentLength
dl = bsxfun(@rdivide,dl,l./(l - segmentLength));
dl_max = segmentLength;
dl_nl = min(abs(l - segmentLength),0.99*dl_max);% non-linear contribution saturates to avoid force decrease for
% l - segL > segLmax, which could happen for finite timestep
nl = 1./abs(1 - (dl_nl/dl_max).^2); % non-linear part of spring
% % Fl = k_l.*(cat(1,dl.*nl,zeros(1,2)) - cat(1,zeros(1,2),dl.*nl)); % add forces to next and previous nodes shifted
Fl = k_l.*(cat(1,bsxfun(@times,dl,nl),zeros(1,2)) ...
- cat(1,zeros(1,2),bsxfun(@times,dl,nl))); % add forces to next and previous nodes shifted
else
Fl = zeros(size(Fm));
end
% bending constraints - rotational springs with changing rest angle due to active undulations
if k_theta(objCtr)>0&&M>2
dK_max = pi;
if theta_0>0
dK = (diff(bendingAngles) - diff(sin(phaseOffset(objCtr,:))'));
else
dK = diff(bendingAngles);
end
dK_nl = min(abs(dK),0.99*dK_max);% non-linear contribution saturates to avoid force decrease for
% dK > dK_max, which could happen for finite timestep
torques = k_theta(objCtr).*dK./(1 - dK_nl.^2/dK_max^2);
torques = torques(2:end); % no rotational springs at head and tail node
% % uncomment for debugging...
% plot(gradient(bodyAngles)), hold on, plot(gradient(sin(phaseOffset(objCtr,:))'))
% % momentsfwd = torques.*l(1:end-1).*e_phi(1:end-1,:);
momentsfwd = bsxfun(@times,torques.*l(1:end-1),e_phi(1:end-1,:));
% % momentsbwd = torques.*l(2:end).*e_phi(2:end,:);
momentsbwd = bsxfun(@times,torques.*l(2:end),e_phi(2:end,:));
% F_theta = NaN(M,2); % pre-allocate to index nodes in order depending on movement state
F_theta = [momentsfwd; 0 0; 0 0] ... % rotational force from node n+1 onto n
+ [0 0; 0 0; momentsbwd] ...% rotational force from node n-1 onto n
+ [0 0; -(momentsfwd + momentsbwd); 0 0];% reactive force on node n (balancing forces exerted onto nodes n+1 and n -1
else
F_theta = zeros(size(Fm));
end
% sum force contributions
forceArray(objCtr,:,:) = Fm + Fl + F_theta;
end
% resolve contact forces
if rc>0||r_LJcutoff>0
if N==40&&M==49 % check if we can use compiled mex function
Fc = resolveContactsLoop_mex(forceArray,distanceMatrixXY,...
distanceMatrix,2*rc,sigma_LJ,r_LJcutoff,eps_LJ, LJnodes, LJmode);
elseif N==40&&M==36 % check if we can use compiled mex function
% Fc = resolveContactsLoop_N40_M36_mex(forceArray,distanceMatrixXY,...
% distanceMatrix,2*rc,sigma_LJ,r_LJcutoff,eps_LJ, LJnodes, LJmode);
Fc = resolveContactsLoop(forceArray,distanceMatrixXY,...
distanceMatrix,2*rc,sigma_LJ,r_LJcutoff,eps_LJ, LJnodes, LJmode);
elseif N==40&&M==18 % check if we can use compiled mex function
% % Fc = resolveContactsLoop_N40_M18_mex(forceArray,distanceMatrixXY,...
% % distanceMatrix,2*rc,sigma_LJ,r_LJcutoff,eps_LJ, LJnodes, LJmode);
Fc = resolveContactsLoop(forceArray,distanceMatrixXY,...
distanceMatrix,2*rc,sigma_LJ,r_LJcutoff,eps_LJ, LJnodes, LJmode);
elseif N==200&&M==18 % check if we can use compiled mex function
Fc = resolveContactsLoop_N200_M18_mex(forceArray,distanceMatrixXY,...
distanceMatrix,2*rc,sigma_LJ,r_LJcutoff,eps_LJ, LJnodes, LJmode);
elseif N==200&&M==36 % check if we can use compiled mex function
Fc = resolveContactsLoop_N200_M36_mex(forceArray,distanceMatrixXY,...
distanceMatrix,2*rc,sigma_LJ,r_LJcutoff,eps_LJ, LJnodes, LJmode);
else
Fc = NaN(N,M,2);
for objCtr = 1:N
for nodeCtr = 1:M
Fc(objCtr,nodeCtr,:) = resolveContacts(forceArray,distanceMatrixXY(:,:,:,objCtr,nodeCtr),...
distanceMatrix(:,:,objCtr,nodeCtr),objCtr,nodeCtr,2*rc,...
sigma_LJ,r_LJcutoff,eps_LJ,LJnodes,LJmode); % factor of two so that rc is node radius
end
end
end
forceArray = forceArray + Fc;
end