RK4 forward and back trajectories from the initial conditions. Adapted from Bokeh's streamline -uses Runge-Kutta method to fill x and y trajectories then checks length of traj (s in units of axes)
(self, x0, y0)
| 194 | return a0 * (1 - yt) + a1 * yt |
| 195 | |
| 196 | def rk4_integrate(self, x0, y0): |
| 197 | """ |
| 198 | RK4 forward and back trajectories from the initial conditions. |
| 199 | |
| 200 | Adapted from Bokeh's streamline -uses Runge-Kutta method to fill |
| 201 | x and y trajectories then checks length of traj (s in units of axes) |
| 202 | """ |
| 203 | |
| 204 | def f(xi, yi): |
| 205 | dt_ds = 1.0 / self.value_at(self.speed, xi, yi) |
| 206 | ui = self.value_at(self.u, xi, yi) |
| 207 | vi = self.value_at(self.v, xi, yi) |
| 208 | return ui * dt_ds, vi * dt_ds |
| 209 | |
| 210 | def g(xi, yi): |
| 211 | dt_ds = 1.0 / self.value_at(self.speed, xi, yi) |
| 212 | ui = self.value_at(self.u, xi, yi) |
| 213 | vi = self.value_at(self.v, xi, yi) |
| 214 | return -ui * dt_ds, -vi * dt_ds |
| 215 | |
| 216 | def check(xi, yi): |
| 217 | return (0 <= xi < len(self.x) - 1) and (0 <= yi < len(self.y) - 1) |
| 218 | |
| 219 | xb_changes = [] |
| 220 | yb_changes = [] |
| 221 | |
| 222 | def rk4(x0, y0, f): |
| 223 | ds = 0.01 |
| 224 | stotal = 0 |
| 225 | xi = x0 |
| 226 | yi = y0 |
| 227 | xb, yb = self.blank_pos(xi, yi) |
| 228 | xf_traj = [] |
| 229 | yf_traj = [] |
| 230 | while check(xi, yi): |
| 231 | xf_traj.append(xi) |
| 232 | yf_traj.append(yi) |
| 233 | try: |
| 234 | k1x, k1y = f(xi, yi) |
| 235 | k2x, k2y = f(xi + 0.5 * ds * k1x, yi + 0.5 * ds * k1y) |
| 236 | k3x, k3y = f(xi + 0.5 * ds * k2x, yi + 0.5 * ds * k2y) |
| 237 | k4x, k4y = f(xi + ds * k3x, yi + ds * k3y) |
| 238 | except IndexError: |
| 239 | break |
| 240 | xi += ds * (k1x + 2 * k2x + 2 * k3x + k4x) / 6.0 |
| 241 | yi += ds * (k1y + 2 * k2y + 2 * k3y + k4y) / 6.0 |
| 242 | if not check(xi, yi): |
| 243 | break |
| 244 | stotal += ds |
| 245 | new_xb, new_yb = self.blank_pos(xi, yi) |
| 246 | if new_xb != xb or new_yb != yb: |
| 247 | if self.blank[new_yb, new_xb] == 0: |
| 248 | self.blank[new_yb, new_xb] = 1 |
| 249 | xb_changes.append(new_xb) |
| 250 | yb_changes.append(new_yb) |
| 251 | xb = new_xb |
| 252 | yb = new_yb |
| 253 | else: |