Geometry differs; all logic below is identical
hap: H = R − AG − AP · normal radial outward | jet: H = Z − AP · normal vertical (Z axis)
Initialise particle (J)
Compute HFRIC & HMIN
brentq root-find; set DELTAMAX
while ATTACHK == 0
Compute forces
VDW · EDL · Born · AB · Steric · Drag · Diff · Gravity · Lift
Heterodomain EDL (HFLAG > 1)
Compute AFRACT; blend FEDL if SCOV > 0
Which zone? H vs HNS / HFRIC
H > HNS
HFRIC < H ≤ HNS
H ≤ HFRIC
Bulk (HFLAG=1)
dT = MULTB × dTMRT
Near-surface (HFLAG=2)
dT = MULTNS × dTMRT; IREF1++
Contact (HFLAG=3)
dT = MULTC × dTMRT; IREF2++
Exit domain?
hap: R > RB | jet: R > REXIT
NS slow-motion? (IREF1 = 1000)
3D Δpos < DFACTNS × diffusion scale
ATTMODE?
0 = perfect sink; 1 = torque-balance
yes
yes
ATTACHK = 1
Exit domain
ATTACHK = 5
Retained — near-surface slow-motion
0
1
ATTACHK=2
Perfect sink
Torque balance
Update UT via
rolling eqn.
UT < 0?
Clamp; set ARRESTFLAG
FADH ≈ FREP?
±0.5%; ARRESTFLAG=1
reset if not yet at equil.
yes
ATTACHK = 2
Arrest
no / continue
Contact slow-motion?
IREF2=1000; tangential
Δpos < DFACTC × diffusion
yes
ATTACHK = 4
Contact slow-motion
H < H0?
ATTACHK = 6 crash
no attachment yet
PTIMEF > TTIME?
Simulation clock expired
yes
ATTACHK = 3
Remaining — time expired
Integrate motion
Update UN, UT, OMEGA, X, Y, Z, R, H
next step
no
no → continue
Per-step speed: traj_jet.py vs traj_hap.py
traj_jet.py is faster per timestep due to simpler geometry — not a code quality difference
Flow field
jet: polynomial, 14 coeffs pre-computed once
hap: analytic Stokes evaluated every step
Gravity
jet: scalar added to UZ
hap: gravvect() — full normal/tangential decomposition per step
Velocity
jet: 3 scalar updates (UX, UY, UZ)
hap: decompose → normal/tangential update → recompose
Heterodomains
jet: 2D polar on flat surface (no rotation)
hap: quaternion 3D rotation to sphere tangent plane — evaluated every NS step
Outcomes
1 Exit
2 Attached — arrest/sink
3 Remaining — timeout
4 Attached — contact slow-motion
5 Retained — near-surface slow-motion
6 Crashed — H < H0