0
# Dust Physics Functions
1

2
Standard functions for dust physics including particle dynamics, coagulation kernels, velocities, probabilities, fluxes, and source terms. These functions implement the core physics of dust evolution in protoplanetary disks with performance-critical calculations accelerated by Fortran extensions.
3

4
## Capabilities
5

6
### Core Dust Properties
7

8
Functions that calculate fundamental dust properties from the simulation state.
9

10
```python { .api }
11
from dustpy.std import dust
12

13
def dust.a(sim):
14
    """
15
    Calculate particle size from density and filling factor.
16

17
    Parameters:
18
    - sim: Simulation object
19

20
    Returns:
21
    numpy.ndarray: Particle sizes [cm]
22
    """
23

24
def dust.rho_midplane(sim):
25
    """
26
    Calculate midplane mass density of dust.
27

28
    Parameters:
29
    - sim: Simulation object
30

31
    Returns:
32
    numpy.ndarray: Midplane mass density [g/cm³]
33
    """
34

35
def dust.eps(sim):
36
    """
37
    Calculate vertically integrated dust-to-gas ratio.
38

39
    Parameters:
40
    - sim: Simulation object
41

42
    Returns:
43
    numpy.ndarray: Dust-to-gas ratio
44
    """
45

46
def dust.H(sim):
47
    """
48
    Calculate dust scale height using Dubrulle et al. (1995) prescription.
49

50
    Parameters:
51
    - sim: Simulation object
52

53
    Returns:
54
    numpy.ndarray: Dust scale heights [cm]
55
    """
56

57
def dust.St_Epstein_StokesI(sim):
58
    """
59
    Calculate Stokes number in Epstein and Stokes I regimes.
60

61
    Parameters:
62
    - sim: Simulation object
63

64
    Returns:
65
    numpy.ndarray: Stokes numbers
66
    """
67
```
68

69
### Transport Properties
70

71
Functions for calculating dust diffusivity and transport coefficients.
72

73
```python { .api }
74
def dust.D(sim):
75
    """
76
    Calculate dust diffusivity from turbulent mixing.
77

78
    Parameters:
79
    - sim: Simulation object
80

81
    Returns:
82
    numpy.ndarray: Dust diffusivity [cm²/s]
83
    """
84
```
85

86
### Collision Kernel
87

88
Functions for calculating particle collision rates and coagulation kernels.
89

90
```python { .api }
91
def dust.kernel(sim):
92
    """
93
    Calculate vertically integrated collision kernel.
94

95
    The collision kernel determines the rate at which particles
96
    of different masses collide and potentially stick or fragment.
97

98
    Parameters:
99
    - sim: Simulation object
100

101
    Returns:
102
    numpy.ndarray: Collision kernel [cm²/s]
103
    """
104
```
105

106
### Velocity Components
107

108
Functions for calculating various components of dust particle velocities.
109

110
```python { .api }
111
def dust.vrad(sim):
112
    """
113
    Calculate radial dust velocity including drift and diffusion.
114

115
    Parameters:
116
    - sim: Simulation object
117

118
    Returns:
119
    numpy.ndarray: Radial velocity [cm/s]
120
    """
121

122
def dust.vdriftmax(sim):
123
    """
124
    Calculate maximum drift velocity including back reaction effects.
125

126
    Parameters:
127
    - sim: Simulation object
128

129
    Returns:
130
    numpy.ndarray: Maximum drift velocity [cm/s]
131
    """
132

133
def dust.vrel_tot(sim):
134
    """
135
    Calculate total relative velocity between particles (RMS).
136

137
    Combines contributions from various relative velocity sources.
138

139
    Parameters:
140
    - sim: Simulation object
141

142
    Returns:
143
    numpy.ndarray: Total relative velocity [cm/s]
144
    """
145

146
def dust.vrel_azimuthal_drift(sim):
147
    """
148
    Calculate relative velocity from azimuthal drift differences.
149

150
    Parameters:
151
    - sim: Simulation object
152

153
    Returns:
154
    numpy.ndarray: Relative velocity from azimuthal drift [cm/s]
155
    """
156

157
def dust.vrel_brownian_motion(sim):
158
    """
159
    Calculate relative velocity from Brownian motion.
160

161
    Parameters:
162
    - sim: Simulation object
163

164
    Returns:
165
    numpy.ndarray: Relative velocity from Brownian motion [cm/s]
166
    """
167

168
def dust.vrel_radial_drift(sim):
169
    """
170
    Calculate relative velocity from radial drift differences.
171

172
    Parameters:
173
    - sim: Simulation object
174

175
    Returns:
176
    numpy.ndarray: Relative velocity from radial drift [cm/s]
177
    """
178

179
def dust.vrel_turbulent_motion(sim):
180
    """
181
    Calculate relative velocity from turbulent motion.
182

183
    Parameters:
184
    - sim: Simulation object
185

186
    Returns:
187
    numpy.ndarray: Relative velocity from turbulence [cm/s]
188
    """
189

190
def dust.vrel_vertical_settling(sim):
191
    """
192
    Calculate relative velocity from vertical settling differences.
193

194
    Parameters:
195
    - sim: Simulation object
196

197
    Returns:
198
    numpy.ndarray: Relative velocity from vertical settling [cm/s]
199
    """
200
```
201

202
### Collision Probabilities
203

204
Functions for calculating sticking and fragmentation probabilities.
205

206
```python { .api }
207
def dust.p_stick(sim):
208
    """
209
    Calculate sticking probability for particle collisions.
210

211
    Determines the likelihood that colliding particles will
212
    stick together rather than bounce or fragment.
213

214
    Parameters:
215
    - sim: Simulation object
216

217
    Returns:
218
    numpy.ndarray: Sticking probability (0-1)
219
    """
220

221
def dust.p_frag(sim):
222
    """
223
    Calculate fragmentation probability for particle collisions.
224

225
    Determines the likelihood that colliding particles will
226
    fragment into smaller pieces.
227

228
    Parameters:
229
    - sim: Simulation object
230

231
    Returns:
232
    numpy.ndarray: Fragmentation probability (0-1)
233
    """
234
```
235

236
### Mass Fluxes
237

238
Functions for calculating mass transport through advection and diffusion.
239

240
```python { .api }
241
def dust.F_adv(sim, Sigma=None):
242
    """
243
    Calculate advective mass flux.
244

245
    Parameters:
246
    - sim: Simulation object
247
    - Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
248

249
    Returns:
250
    numpy.ndarray: Advective flux [g/cm/s]
251
    """
252

253
def dust.F_diff(sim, Sigma=None):
254
    """
255
    Calculate diffusive mass flux.
256

257
    Parameters:
258
    - sim: Simulation object
259
    - Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
260

261
    Returns:
262
    numpy.ndarray: Diffusive flux [g/cm/s]
263
    """
264

265
def dust.F_tot(sim, Sigma=None):
266
    """
267
    Calculate total mass flux (advective + diffusive).
268

269
    Parameters:
270
    - sim: Simulation object
271
    - Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
272

273
    Returns:
274
    numpy.ndarray: Total mass flux [g/cm/s]
275
    """
276
```
277

278
### Source Terms
279

280
Functions for calculating mass source and sink terms from various physical processes.
281

282
```python { .api }
283
def dust.S_coag(sim, Sigma=None):
284
    """
285
    Calculate coagulation source terms.
286

287
    Computes mass exchange between size bins due to particle
288
    collisions, sticking, and fragmentation processes.
289

290
    Parameters:
291
    - sim: Simulation object
292
    - Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
293

294
    Returns:
295
    numpy.ndarray: Coagulation source terms [g/cm²/s]
296
    """
297

298
def dust.S_hyd(sim, Sigma=None):
299
    """
300
    Calculate hydrodynamic source terms.
301

302
    Computes source terms from hydrodynamic processes including
303
    radial transport and diffusion.
304

305
    Parameters:
306
    - sim: Simulation object
307
    - Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
308

309
    Returns:
310
    numpy.ndarray: Hydrodynamic source terms [g/cm²/s]
311
    """
312

313
def dust.S_tot(sim, Sigma=None):
314
    """
315
    Calculate total source terms (coagulation + hydrodynamic).
316

317
    Parameters:
318
    - sim: Simulation object
319
    - Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
320

321
    Returns:
322
    numpy.ndarray: Total source terms [g/cm²/s]
323
    """
324
```
325

326
### Initial Conditions
327

328
Functions for setting up initial dust distributions and parameters.
329

330
```python { .api }
331
def dust.MRN_distribution(sim):
332
    """
333
    Set up MRN (Mathis, Rumpl, Nordsieck) initial particle mass distribution.
334

335
    Creates a power-law size distribution commonly used as initial
336
    conditions for dust evolution simulations.
337

338
    Parameters:
339
    - sim: Simulation object
340

341
    Returns:
342
    numpy.ndarray: Initial surface density distribution [g/cm²]
343
    """
344

345
def dust.SigmaFloor(sim):
346
    """
347
    Calculate floor value for dust surface density.
348

349
    Prevents numerical issues by enforcing minimum dust densities.
350

351
    Parameters:
352
    - sim: Simulation object
353

354
    Returns:
355
    numpy.ndarray: Floor values [g/cm²]
356
    """
357
```
358

359
### Integration Functions
360

361
Functions for numerical integration of dust evolution equations.
362

363
```python { .api }
364
def dust.jacobian(sim, x, dx=None):
365
    """
366
    Calculate Jacobian matrix for implicit integration schemes.
367

368
    Parameters:
369
    - sim: Simulation object
370
    - x: Current state vector
371
    - dx: State perturbation (optional)
372

373
    Returns:
374
    numpy.ndarray: Jacobian matrix
375
    """
376

377
def dust.Sigma_deriv(sim, t, Sigma):
378
    """
379
    Calculate time derivative of dust surface density.
380

381
    Core function for explicit integration schemes.
382

383
    Parameters:
384
    - sim: Simulation object
385
    - t: Current time [s]
386
    - Sigma: Current surface density [g/cm²]
387

388
    Returns:
389
    numpy.ndarray: Time derivative [g/cm²/s]
390
    """
391

392
class dust.impl_1_direct:
393
    """
394
    Direct implicit integration scheme for dust evolution.
395

396
    Provides first-order implicit integration with direct
397
    matrix solving for stable evolution of stiff systems.
398
    """
399
    def __init__(self):
400
        """Initialize implicit integration scheme."""
401
        ...
402
```
403

404
### Utility Functions
405

406
Helper functions for dust simulation setup and maintenance.
407

408
```python { .api }
409
def dust.boundary(sim):
410
    """
411
    Set boundary conditions for dust quantities.
412

413
    Applies inner and outer boundary conditions based on
414
    the configured boundary condition types and values.
415

416
    Parameters:
417
    - sim: Simulation object
418
    """
419

420
def dust.enforce_floor_value(sim):
421
    """
422
    Enforce minimum floor values for dust quantities.
423

424
    Prevents numerical issues by ensuring all dust quantities
425
    remain above specified minimum values.
426

427
    Parameters:
428
    - sim: Simulation object
429
    """
430

431
def dust.coagulation_parameters(sim):
432
    """
433
    Calculate and update coagulation-related parameters.
434

435
    Updates collision kernels, probabilities, and other
436
    coagulation parameters based on current simulation state.
437

438
    Parameters:
439
    - sim: Simulation object
440
    """
441

442
def dust.dt(sim):
443
    """
444
    Calculate appropriate time step for dust evolution.
445

446
    Determines stable time step based on current dust
447
    evolution rates and numerical stability criteria.
448

449
    Parameters:
450
    - sim: Simulation object
451

452
    Returns:
453
    float: Time step [s]
454
    """
455

456
def dust.dt_adaptive(sim):
457
    """
458
    Calculate adaptive time step for dust evolution.
459

460
    Uses error estimation to adjust time step dynamically
461
    for optimal balance of accuracy and efficiency.
462

463
    Parameters:
464
    - sim: Simulation object
465

466
    Returns:
467
    float: Adaptive time step [s]
468
    """
469

470
def dust.prepare(sim):
471
    """
472
    Prepare implicit dust integration step.
473

474
    Stores the current value of surface density in a hidden
475
    field for use in implicit integration schemes.
476

477
    Parameters:
478
    - sim: Simulation object
479
    """
480

481
def dust.finalize_explicit(sim):
482
    """
483
    Finalize explicit integration step.
484

485
    Applies boundary conditions and enforces floor values
486
    after an explicit integration step.
487

488
    Parameters:
489
    - sim: Simulation object
490
    """
491

492
def dust.finalize_implicit(sim):
493
    """
494
    Finalize implicit integration step.
495

496
    Applies boundary conditions, enforces floor values, and
497
    updates dependent quantities after implicit integration.
498

499
    Parameters:
500
    - sim: Simulation object
501
    """
502

503
def dust.set_implicit_boundaries(sim):
504
    """
505
    Calculate fluxes at boundaries after implicit integration.
506

507
    Sets boundary source terms based on the change in surface
508
    density during the implicit integration step.
509

510
    Parameters:
511
    - sim: Simulation object
512
    """
513
```
514

515
## Usage Examples
516

517
### Basic Dust Property Calculation
518

519
```python
520
from dustpy import Simulation
521
from dustpy.std import dust
522

523
# Create and initialize simulation
524
sim = Simulation()
525
sim.makegrids()
526
sim.initialize()
527

528
# Calculate dust properties
529
particle_sizes = dust.a(sim)              # [cm]
530
stokes_numbers = dust.St_Epstein_StokesI(sim)
531
dust_scale_heights = dust.H(sim)          # [cm]
532
collision_kernel = dust.kernel(sim)       # [cm²/s]
533
```
534

535
### Velocity Analysis
536

537
```python
538
# Calculate various velocity components
539
radial_velocity = dust.vrad(sim)          # [cm/s]
540
total_relative_v = dust.vrel_tot(sim)     # [cm/s]
541
brownian_v = dust.vrel_brownian_motion(sim)  # [cm/s]
542
drift_v = dust.vrel_radial_drift(sim)     # [cm/s]
543

544
# Maximum drift velocity with back reaction
545
max_drift = dust.vdriftmax(sim)           # [cm/s]
546
```
547

548
### Source Term Analysis
549

550
```python
551
# Calculate source terms
552
coag_sources = dust.S_coag(sim)           # [g/cm²/s]
553
hydro_sources = dust.S_hyd(sim)           # [g/cm²/s]
554
total_sources = dust.S_tot(sim)           # [g/cm²/s]
555

556
# Analyze mass fluxes
557
advective_flux = dust.F_adv(sim)          # [g/cm/s]
558
diffusive_flux = dust.F_diff(sim)         # [g/cm/s]
559
```
560

561
### Custom Time Stepping
562

563
```python
564
# Calculate appropriate time steps
565
dt_dust = dust.dt(sim)                    # [s]
566
dt_adaptive = dust.dt_adaptive(sim)       # [s]
567

568
print(f"Dust time step: {dt_dust/dust.year:.2e} years")
569
print(f"Adaptive time step: {dt_adaptive/dust.year:.2e} years")
570
```