|
77 | 77 | " interaction_distance = 0.5\n", |
78 | 78 | " velocity = -0.04 # predefined attracting velocity\n", |
79 | 79 | "\n", |
80 | | - " lon_a = particles.lon[np.where(particles.attractor)[0]]\n", |
| 80 | + " # Boolean mask and coordinates of attractors\n", |
| 81 | + " attractor_mask = particles.attractor.astype(bool)\n", |
| 82 | + " lon_a = particles.lon[attractor_mask]\n", |
| 83 | + " lat_a = particles.lat[attractor_mask]\n", |
| 84 | + "\n", |
| 85 | + " # Pairwise differences and distances (n_attractors × n_particles)\n", |
81 | 86 | " dx = particles.lon - lon_a[:, None]\n", |
82 | | - " lat_a = particles.lat[np.where(particles.attractor)[0]]\n", |
83 | 87 | " dy = particles.lat - lat_a[:, None]\n", |
| 88 | + " distances = np.sqrt(dx**2 + dy**2)\n", |
84 | 89 | "\n", |
85 | | - " distances = np.sqrt(dx**2 + dy**2) # TODO use np.linalg.norm?\n", |
| 90 | + " # Mask dx, dy by interaction range\n", |
| 91 | + " within = distances < interaction_distance\n", |
| 92 | + " dx = np.where(within, dx, 0.0)\n", |
| 93 | + " dy = np.where(within, dy, 0.0)\n", |
86 | 94 | "\n", |
87 | | - " dx = np.where(distances < interaction_distance, dx, 0.0)\n", |
88 | | - " dy = np.where(distances < interaction_distance, dy, 0.0)\n", |
| 95 | + " with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n", |
| 96 | + " inv_dist = np.where(distances > 0, 1.0 / distances, 0.0)\n", |
| 97 | + " dx_norm = dx * inv_dist\n", |
| 98 | + " dy_norm = dy * inv_dist\n", |
89 | 99 | "\n", |
90 | | - " dx_norm = np.where(distances > 0, dx / distances, 0)\n", |
91 | 100 | " particles.dlon += np.sum(dx_norm, axis=0) * velocity * particles.dt\n", |
92 | | - "\n", |
93 | | - " dy_norm = np.where(distances > 0, dy / distances, 0)\n", |
94 | 101 | " particles.dlat += np.sum(dy_norm, axis=0) * velocity * particles.dt" |
95 | 102 | ] |
96 | 103 | }, |
|
0 commit comments