[7837] | 1 | """boundary.py - Classes for implementing boundary conditions |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | class Boundary: |
---|
| 6 | """Base class for boundary conditions. |
---|
| 7 | Specific boundary conditions must provide values for |
---|
| 8 | the conserved_quantities |
---|
| 9 | |
---|
| 10 | A boundary object has one neighbour; the one it |
---|
| 11 | serves. |
---|
| 12 | """ |
---|
| 13 | |
---|
| 14 | def __init__(self): |
---|
| 15 | pass |
---|
| 16 | |
---|
| 17 | def evaluate(self, vol_id=None, edge_id=None): |
---|
| 18 | msg = 'Generic class Boundary must be subclassed' |
---|
| 19 | raise msg |
---|
| 20 | |
---|
| 21 | class Reflective_boundary(Boundary): |
---|
| 22 | """Reflective boundary returns same conserved quantities as |
---|
| 23 | those present in its neighbour volume but reflected. |
---|
| 24 | |
---|
| 25 | This class is specific to the shallow water equation as it |
---|
| 26 | works with the momentum quantities assumed to be the second |
---|
| 27 | and third conserved quantities. |
---|
| 28 | """ |
---|
| 29 | |
---|
| 30 | def __init__(self, domain = None): |
---|
| 31 | Boundary.__init__(self) |
---|
| 32 | |
---|
| 33 | if domain is None: |
---|
| 34 | msg = 'Domain must be specified for reflective boundary' |
---|
| 35 | raise msg |
---|
| 36 | |
---|
| 37 | #Handy shorthands |
---|
| 38 | self.normals = domain.normals |
---|
| 39 | self.stage = domain.quantities['stage'].vertex_values |
---|
| 40 | self.xmom = domain.quantities['xmomentum'].vertex_values |
---|
| 41 | self.bed = domain.quantities['elevation'].vertex_values |
---|
| 42 | self.height = domain.quantities['height'].vertex_values |
---|
| 43 | self.velocity = domain.quantities['velocity'].vertex_values |
---|
| 44 | |
---|
| 45 | from Numeric import zeros, Float |
---|
| 46 | self.evolved_quantities = zeros(5, Float) |
---|
| 47 | |
---|
| 48 | def __repr__(self): |
---|
| 49 | return 'Reflective_boundary' |
---|
| 50 | |
---|
| 51 | def evaluate(self, vol_id, edge_id): |
---|
| 52 | """Reflective boundaries reverses the outward momentum |
---|
| 53 | of the volume they serve. |
---|
| 54 | """ |
---|
| 55 | q = self.evolved_quantities |
---|
| 56 | q[0] = self.stage[vol_id, edge_id] |
---|
| 57 | q[1] = -self.xmom[vol_id, edge_id] |
---|
| 58 | q[2] = self.bed[vol_id, edge_id] |
---|
| 59 | q[3] = self.height[vol_id, edge_id] |
---|
| 60 | q[4] = -self.velocity[vol_id, edge_id] |
---|
| 61 | return q |
---|
| 62 | |
---|
| 63 | |
---|
| 64 | class Transmissive_boundary(Boundary): |
---|
| 65 | """Transmissive boundary returns same conserved quantities as |
---|
| 66 | those present in its neighbour volume. |
---|
| 67 | |
---|
| 68 | Underlying domain must be specified when boundary is instantiated |
---|
| 69 | """ |
---|
| 70 | |
---|
| 71 | def __init__(self, domain = None): |
---|
| 72 | Boundary.__init__(self) |
---|
| 73 | |
---|
| 74 | if domain is None: |
---|
| 75 | msg = 'Domain must be specified for transmissive boundary' |
---|
| 76 | raise msg |
---|
| 77 | |
---|
| 78 | self.domain = domain |
---|
| 79 | |
---|
| 80 | def __repr__(self): |
---|
| 81 | return 'Transmissive_boundary(%s)' %self.domain |
---|
| 82 | |
---|
| 83 | def evaluate(self, vol_id, vertex_id): |
---|
| 84 | """Transmissive boundaries return the edge values |
---|
| 85 | of the volume they serve. |
---|
| 86 | """ |
---|
| 87 | |
---|
| 88 | q = self.domain.get_conserved_quantities(vol_id, vertex = vertex_id)#(vol_id, edge = edge_id) |
---|
| 89 | return q |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | class Dirichlet_boundary(Boundary): |
---|
| 93 | """Dirichlet boundary returns constant values for the |
---|
| 94 | conserved quantities |
---|
| 95 | """ |
---|
| 96 | |
---|
| 97 | |
---|
| 98 | def __init__(self, conserved_quantities=None): |
---|
| 99 | Boundary.__init__(self) |
---|
| 100 | |
---|
| 101 | if conserved_quantities is None: |
---|
| 102 | msg = 'Must specify one value for each conserved quantity' |
---|
| 103 | raise msg |
---|
| 104 | |
---|
| 105 | from Numeric import array, Float |
---|
| 106 | self.conserved_quantities=array(conserved_quantities).astype(Float) |
---|
| 107 | |
---|
| 108 | def __repr__(self): |
---|
| 109 | return 'Dirichlet boundary (%s)' %self.conserved_quantities |
---|
| 110 | |
---|
| 111 | def evaluate(self, vol_id=None, edge_id=None): |
---|
| 112 | return self.conserved_quantities |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | |
---|
| 116 | class Time_boundary(Boundary): |
---|
| 117 | """Time dependent boundary returns values for the |
---|
| 118 | conserved quantities as a function of time. |
---|
| 119 | Must specify domain to get access to model time and a function f(t) |
---|
| 120 | which must return conserved quantities as a function time |
---|
| 121 | """ |
---|
| 122 | |
---|
| 123 | |
---|
| 124 | def __init__(self, domain = None, f = None): |
---|
| 125 | Boundary.__init__(self) |
---|
| 126 | |
---|
| 127 | try: |
---|
| 128 | q = f(0.0) |
---|
| 129 | except Exception, e: |
---|
| 130 | msg = 'Function for time boundary could not be executed:\n%s' %e |
---|
| 131 | raise msg |
---|
| 132 | |
---|
| 133 | |
---|
| 134 | from Numeric import array, Float |
---|
| 135 | try: |
---|
| 136 | q = array(q).astype(Float) |
---|
| 137 | except: |
---|
| 138 | msg = 'Return value from time boundary function could ' |
---|
| 139 | msg += 'not be converted into a Numeric array of floats.\n' |
---|
| 140 | msg += 'Specified function should return either list or array.' |
---|
| 141 | raise msg |
---|
| 142 | |
---|
| 143 | msg = 'ERROR: Time boundary function must return a 1d list or array ' |
---|
| 144 | assert len(q.shape) == 1, msg |
---|
| 145 | |
---|
| 146 | d = len(domain.evolved_quantities) |
---|
| 147 | |
---|
| 148 | #print 'd = ',d |
---|
| 149 | msg = 'Return value for f must be a list or an array of length %d' %d |
---|
| 150 | assert len(q) == d, msg |
---|
| 151 | |
---|
| 152 | self.f = f |
---|
| 153 | self.domain = domain |
---|
| 154 | |
---|
| 155 | def __repr__(self): |
---|
| 156 | return 'Time boundary' |
---|
| 157 | |
---|
| 158 | def evaluate(self, vol_id=None, edge_id=None): |
---|
| 159 | return self.f(self.domain.time) |
---|
| 160 | |
---|
| 161 | |
---|
| 162 | class File_boundary_time(Boundary): |
---|
| 163 | """Boundary values obtained from file and interpolated. |
---|
| 164 | conserved quantities as a function of time. |
---|
| 165 | |
---|
| 166 | Assumes that file contains a time series. |
---|
| 167 | |
---|
| 168 | No spatial info assumed. |
---|
| 169 | """ |
---|
| 170 | |
---|
| 171 | #FIXME: Is this still necessary |
---|
| 172 | |
---|
| 173 | def __init__(self, filename, domain): |
---|
| 174 | import time |
---|
| 175 | from Numeric import array |
---|
| 176 | from config import time_format |
---|
| 177 | from util import File_function |
---|
| 178 | |
---|
| 179 | Boundary.__init__(self) |
---|
| 180 | |
---|
| 181 | self.F = File_function(filename, domain) |
---|
| 182 | self.domain = domain |
---|
| 183 | |
---|
| 184 | #Test |
---|
| 185 | q = self.F(0) |
---|
| 186 | |
---|
| 187 | d = len(domain.conserved_quantities) |
---|
| 188 | msg = 'Values specified in file must be a list or an array of length %d' %d |
---|
| 189 | assert len(q) == d, msg |
---|
| 190 | |
---|
| 191 | |
---|
| 192 | def __repr__(self): |
---|
| 193 | return 'File boundary' |
---|
| 194 | |
---|
| 195 | def evaluate(self, vol_id=None, edge_id=None): |
---|
| 196 | """Return linearly interpolated values based on domain.time |
---|
| 197 | |
---|
| 198 | vol_id and edge_id are ignored |
---|
| 199 | """ |
---|
| 200 | |
---|
| 201 | t = self.domain.time |
---|
| 202 | return self.F(t) |
---|
| 203 | |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | |
---|
| 207 | class File_boundary(Boundary): |
---|
| 208 | """Boundary values obtained from file and interpolated. |
---|
| 209 | conserved quantities as a function of time. |
---|
| 210 | |
---|
| 211 | Assumes that file contains a time series and possibly |
---|
| 212 | also spatial info. |
---|
| 213 | See docstring for File_function in util.py for details about |
---|
| 214 | admissible file formats |
---|
| 215 | |
---|
| 216 | The full solution history is not exactly the same as if |
---|
| 217 | the models were coupled: |
---|
| 218 | File boundary must read and interpolate from *smoothed* version |
---|
| 219 | as stored in sww and cannot work with the discontinuos triangles. |
---|
| 220 | |
---|
| 221 | """ |
---|
| 222 | |
---|
| 223 | def __init__(self, filename, domain, verbose = False): |
---|
| 224 | import time |
---|
| 225 | from Numeric import array, zeros, Float |
---|
| 226 | from config import time_format |
---|
| 227 | from util import file_function |
---|
| 228 | |
---|
| 229 | Boundary.__init__(self) |
---|
| 230 | |
---|
| 231 | #Get x,y vertex coordinates for all triangles |
---|
| 232 | V = domain.vertex_coordinates |
---|
| 233 | |
---|
| 234 | #Compute midpoint coordinates for all boundary elements |
---|
| 235 | #Only a subset may be invoked when boundary is evaluated but |
---|
| 236 | #we don't know which ones at this stage since this object can be attached to |
---|
| 237 | #any tagged boundary later on. |
---|
| 238 | |
---|
| 239 | if verbose: print 'Find midpoint coordinates of entire boundary' |
---|
| 240 | self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float) |
---|
| 241 | boundary_keys = domain.boundary.keys() |
---|
| 242 | |
---|
| 243 | |
---|
| 244 | xllcorner = domain.geo_reference.get_xllcorner() |
---|
| 245 | yllcorner = domain.geo_reference.get_yllcorner() |
---|
| 246 | |
---|
| 247 | |
---|
| 248 | #Make ordering unique #FIXME: should this happen in domain.py? |
---|
| 249 | boundary_keys.sort() |
---|
| 250 | |
---|
| 251 | |
---|
| 252 | #Record ordering #FIXME: should this also happen in domain.py? |
---|
| 253 | self.boundary_indices = {} |
---|
| 254 | for i, (vol_id, edge_id) in enumerate(boundary_keys): |
---|
| 255 | |
---|
| 256 | x0 = V[vol_id, 0]; y0 = V[vol_id, 1] |
---|
| 257 | x1 = V[vol_id, 2]; y1 = V[vol_id, 3] |
---|
| 258 | x2 = V[vol_id, 4]; y2 = V[vol_id, 5] |
---|
| 259 | |
---|
| 260 | #Compute midpoints |
---|
| 261 | if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2]) |
---|
| 262 | if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2]) |
---|
| 263 | if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2]) |
---|
| 264 | |
---|
| 265 | #Convert to absolute UTM coordinates |
---|
| 266 | m[0] += xllcorner |
---|
| 267 | m[1] += yllcorner |
---|
| 268 | |
---|
| 269 | #Register point and index |
---|
| 270 | self.midpoint_coordinates[i,:] = m |
---|
| 271 | self.boundary_indices[(vol_id, edge_id)] = i |
---|
| 272 | |
---|
| 273 | |
---|
| 274 | if verbose: print 'Initialise file_function' |
---|
| 275 | self.F = file_function(filename, domain, |
---|
| 276 | interpolation_points=self.midpoint_coordinates, |
---|
| 277 | verbose = verbose) |
---|
| 278 | self.domain = domain |
---|
| 279 | |
---|
| 280 | #Test |
---|
| 281 | q = self.F(0, point_id=0) |
---|
| 282 | |
---|
| 283 | d = len(domain.conserved_quantities) |
---|
| 284 | msg = 'Values specified in file %s must be ' %filename |
---|
| 285 | msg += ' a list or an array of length %d' %d |
---|
| 286 | assert len(q) == d, msg |
---|
| 287 | |
---|
| 288 | |
---|
| 289 | def __repr__(self): |
---|
| 290 | return 'File boundary' |
---|
| 291 | |
---|
| 292 | |
---|
| 293 | def evaluate(self, vol_id=None, edge_id=None): |
---|
| 294 | """Return linearly interpolated values based on domain.time |
---|
| 295 | |
---|
| 296 | vol_id and edge_id are ignored |
---|
| 297 | """ |
---|
| 298 | |
---|
| 299 | t = self.domain.time |
---|
| 300 | |
---|
| 301 | if vol_id is not None and edge_id is not None: |
---|
| 302 | i = self.boundary_indices[ vol_id, edge_id ] |
---|
| 303 | return self.F(t, point_id = i) |
---|
| 304 | else: |
---|
| 305 | #raise 'Boundary call without point_id not implemented' |
---|
| 306 | #FIXME: What should the semantics be? |
---|
| 307 | return self.F(t) |
---|
| 308 | |
---|
| 309 | |
---|
| 310 | |
---|
| 311 | |
---|
| 312 | |
---|
| 313 | |
---|
| 314 | |
---|
| 315 | |
---|
| 316 | #THIS FAR (10/8/4) |
---|
| 317 | class Connective_boundary(Boundary): |
---|
| 318 | """Connective boundary returns values for the |
---|
| 319 | conserved quantities from a volume as defined by a connection table |
---|
| 320 | mapping between tuples of (volume id, face id) for volumes that |
---|
| 321 | have their boundaries connected. |
---|
| 322 | |
---|
| 323 | FIXME: Perhaps include possibility of mapping between |
---|
| 324 | different domains as well |
---|
| 325 | |
---|
| 326 | FIXME: In case of shallow water we may want to have a |
---|
| 327 | special version that casts this in terms of height rather than stage |
---|
| 328 | """ |
---|
| 329 | |
---|
| 330 | |
---|
| 331 | def __init__(self, table): |
---|
| 332 | from domain import Volume |
---|
| 333 | |
---|
| 334 | Boundary.__init__(self) |
---|
| 335 | |
---|
| 336 | self.connection_table = table |
---|
| 337 | self.Volume = Volume |
---|
| 338 | |
---|
| 339 | def __repr__(self): |
---|
| 340 | return 'Connective boundary' |
---|
| 341 | |
---|
| 342 | #FIXME: IF we ever need to get field_values from connected volume, |
---|
| 343 | #that method could be overridden here (using same idea as in |
---|
| 344 | #get_conserved_quantities |
---|
| 345 | #def get_field_values() |
---|
| 346 | |
---|
| 347 | def get_conserved_quantities(self, volume, face=0): |
---|
| 348 | |
---|
| 349 | id = volume.id |
---|
| 350 | if self.connection_table.has_key((id, face)): |
---|
| 351 | other_id, other_face = self.connection_table[(id, face)] |
---|
| 352 | |
---|
| 353 | other_volume = self.Volume.instances[other_id] |
---|
| 354 | cmd = 'q = other_volume.conserved_quantities_face%d' %face; |
---|
| 355 | exec(cmd) |
---|
| 356 | return q |
---|
| 357 | else: |
---|
| 358 | msg = 'Volume, face tuple ($d, %d) has not been mapped'\ |
---|
| 359 | %(id, face) |
---|
| 360 | raise msg |
---|
| 361 | |
---|
| 362 | |
---|
| 363 | |
---|
| 364 | |
---|
| 365 | |
---|
| 366 | #FIXME: Add a boundary with a general function of x,y, and t |
---|
| 367 | |
---|
| 368 | #FIXME: Add periodic boundaries e.g.: |
---|
| 369 | # Attempt at periodic conditions from advection_spik. Remember this |
---|
| 370 | # |
---|
| 371 | #first = 2*(N-1)*N |
---|
| 372 | #for i in range(1,2*N+1,2): |
---|
| 373 | # k = first + i-1# |
---|
| 374 | # |
---|
| 375 | # print i,k |
---|
| 376 | # |
---|
| 377 | # domain[i].faces[2].neighbour = domain[k].faces[1] |
---|
| 378 | # domain[k].faces[1].neighbour = domain[i].faces[2] |
---|
| 379 | |
---|
| 380 | |
---|
| 381 | |
---|
| 382 | class General_boundary(Boundary): |
---|
| 383 | """General boundary which can compute conserved quantities based on |
---|
| 384 | their previous value, conserved quantities of its neighbour and model time. |
---|
| 385 | |
---|
| 386 | Must specify initial conserved quantities, |
---|
| 387 | neighbour, |
---|
| 388 | domain to get access to model time |
---|
| 389 | a function f(q_old, neighbours_q, t) which must return |
---|
| 390 | new conserved quantities q as a function time |
---|
| 391 | |
---|
| 392 | FIXME: COMPLETE UNTESTED - JUST AN IDEA |
---|
| 393 | """ |
---|
| 394 | |
---|
| 395 | def __init__(self, neighbour=None, conserved_quantities=None, domain=None, f=None): |
---|
| 396 | Boundary.__init__(self, neighbour=neighbour, conserved_quantities=conserved_quantities) |
---|
| 397 | |
---|
| 398 | self.f = f |
---|
| 399 | self.domain = domain |
---|
| 400 | |
---|
| 401 | |
---|
| 402 | def get_conserved_quantities(self, volume=None, face=0): |
---|
| 403 | return self.f(self.conserved_quantities, |
---|
| 404 | neighbour.conserved_quantities, |
---|
| 405 | self.domain.time) |
---|
| 406 | |
---|
| 407 | |
---|
| 408 | |
---|
| 409 | |
---|
| 410 | |
---|
| 411 | class File_boundary_old(Boundary): |
---|
| 412 | """Boundary values obtained from file and interpolated. |
---|
| 413 | conserved quantities as a function of time. |
---|
| 414 | |
---|
| 415 | Assumes that file contains a time series. |
---|
| 416 | |
---|
| 417 | No spatial info assumed. |
---|
| 418 | """ |
---|
| 419 | |
---|
| 420 | |
---|
| 421 | def __init__(self, domain = None, filename = None): |
---|
| 422 | import time |
---|
| 423 | from Numeric import array |
---|
| 424 | from config import time_format |
---|
| 425 | |
---|
| 426 | Boundary.__init__(self) |
---|
| 427 | |
---|
| 428 | |
---|
| 429 | try: |
---|
| 430 | fid = open(filename) |
---|
| 431 | except Exception, e: |
---|
| 432 | msg = 'Boundary file %s could not be opened: %s\n' %(filename, e) |
---|
| 433 | raise msg |
---|
| 434 | |
---|
| 435 | |
---|
| 436 | line = fid.readline() |
---|
| 437 | fid.close() |
---|
| 438 | fields = line.split(',') |
---|
| 439 | msg = 'File %s must have the format date, values' |
---|
| 440 | assert len(fields) == 2, msg |
---|
| 441 | |
---|
| 442 | try: |
---|
| 443 | starttime = time.mktime(time.strptime(fields[0], time_format)) |
---|
| 444 | except ValueError: |
---|
| 445 | msg = 'First field in file %s must be' %filename |
---|
| 446 | msg += ' date-time with format %s.\n' %time_format |
---|
| 447 | msg += 'I got %s instead.' %fields[0] |
---|
| 448 | raise msg |
---|
| 449 | |
---|
| 450 | #Split values |
---|
| 451 | values = [] |
---|
| 452 | for value in fields[1].split(): |
---|
| 453 | values.append(float(value)) |
---|
| 454 | |
---|
| 455 | q = array(values) |
---|
| 456 | |
---|
| 457 | msg = 'ERROR: File boundary function must return a 1d list or array ' |
---|
| 458 | assert len(q.shape) == 1, msg |
---|
| 459 | |
---|
| 460 | d = len(domain.conserved_quantities) |
---|
| 461 | msg = 'Values specified in file must be a list or an array of length %d' %d |
---|
| 462 | assert len(q) == d, msg |
---|
| 463 | |
---|
| 464 | self.filename = filename |
---|
| 465 | self.domain = domain |
---|
| 466 | self.starttime = starttime |
---|
| 467 | |
---|
| 468 | if domain.starttime is None: |
---|
| 469 | domain.starttime = starttime |
---|
| 470 | else: |
---|
| 471 | msg = 'Start time as specified in domain (%s) is earlier ' |
---|
| 472 | 'than the starttime of file %s: %s'\ |
---|
| 473 | %(domain.starttime, self.filename, self.starttime) |
---|
| 474 | if self.starttime > domain.starttime: |
---|
| 475 | raise msg |
---|
| 476 | |
---|
| 477 | self.read_time_boundary() #Now read all times in. |
---|
| 478 | |
---|
| 479 | |
---|
| 480 | def read_time_boundary(self): |
---|
| 481 | from Numeric import zeros, Float, alltrue |
---|
| 482 | from config import time_format |
---|
| 483 | import time |
---|
| 484 | |
---|
| 485 | fid = open(self.filename) |
---|
| 486 | lines = fid.readlines() |
---|
| 487 | fid.close() |
---|
| 488 | |
---|
| 489 | N = len(lines) |
---|
| 490 | d = len(self.domain.conserved_quantities) |
---|
| 491 | |
---|
| 492 | T = zeros(N, Float) |
---|
| 493 | Q = zeros((N, d), Float) |
---|
| 494 | |
---|
| 495 | |
---|
| 496 | for i, line in enumerate(lines): |
---|
| 497 | fields = line.split(',') |
---|
| 498 | real_time = time.mktime(time.strptime(fields[0], time_format)) |
---|
| 499 | |
---|
| 500 | T[i] = real_time - self.starttime |
---|
| 501 | |
---|
| 502 | |
---|
| 503 | for j, value in enumerate(fields[1].split()): |
---|
| 504 | Q[i, j] = float(value) |
---|
| 505 | |
---|
| 506 | msg = 'Time boundary must list time as a monotonuosly ' |
---|
| 507 | msg += 'increasing sequence' |
---|
| 508 | |
---|
| 509 | assert alltrue( T[1:] - T[:-1] > 0 ), msg |
---|
| 510 | |
---|
| 511 | self.T = T #Time |
---|
| 512 | self.Q = Q #Boundary values |
---|
| 513 | self.index = 0 #Initial index |
---|
| 514 | |
---|
| 515 | |
---|
| 516 | def __repr__(self): |
---|
| 517 | return 'File boundary' |
---|
| 518 | |
---|
| 519 | def evaluate(self, vol_id=None, edge_id=None): |
---|
| 520 | """Return linearly interpolated values based on domain.time |
---|
| 521 | """ |
---|
| 522 | |
---|
| 523 | t = self.domain.time |
---|
| 524 | |
---|
| 525 | msg = 'Time given in File boundary does not match model time' |
---|
| 526 | if t < self.T[0]: raise msg |
---|
| 527 | if t > self.T[-1]: raise msg |
---|
| 528 | |
---|
| 529 | oldindex = self.index |
---|
| 530 | |
---|
| 531 | #Find slot |
---|
| 532 | while t > self.T[self.index]: self.index += 1 |
---|
| 533 | while t < self.T[self.index]: self.index -= 1 |
---|
| 534 | |
---|
| 535 | #if oldindex != self.index: |
---|
| 536 | # print 'Changing from %d to %d' %(oldindex, self.index) |
---|
| 537 | |
---|
| 538 | |
---|
| 539 | #t is now between index and index+1 |
---|
| 540 | ratio = (t - self.T[self.index])/\ |
---|
| 541 | (self.T[self.index+1] - self.T[self.index]) |
---|
| 542 | |
---|
| 543 | return self.Q[self.index,:] +\ |
---|
| 544 | ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) |
---|
| 545 | |
---|
| 546 | |
---|
| 547 | |
---|
| 548 | |
---|