Data structure and algorithms for 1D velocity model using layers?

Question

This is for a geophysical analysis program I am creating. I already have code to do all this, but I am looking for inspirations and ideas (good datastructures and algorithms).

What I want to model:

• Velocity as a function of depth (z)
• The model is built up from multiple layers (<10)
  • Every layer is accessible by an index going from 0 for the top most layer to n for the bottom most layer
• Every layer has velocity as a linear function of depth (gradient a_k and axis intercept b_k of the kth layer)
• Every layer has a top and bottom depth (z_k-1 and z_k)
• The model is complete, there is no space between layers. The point directly between two layers belongs to the lower layer

Requirements:

• Get velocity at an arbitrary depth within the model. This will be done on the order of 1k to 10k times, so it should be well optimized.
• Access to the top and bottom depths, gradients and intercepts of a layer by the layer index

What I have so far:
I have working Python code where every layer is saved as a numpy array with the values of z_k (bottom depth), z_k-1 (top depth), a_k (velocity gradient) and b_k (axis intercept). To evaluate the model at a certain depth, I get the layer index (, use that to get the parameters of the layer and pass them to a function that evaluates the linear velocity gradient.

Answer 1

So you have piecewise linear dependence, where z-coordinates of pieces ends go irrregular, and want to get function value at given z.

Note that there is no sense to use binary search for 10 pieces (3-4 rounds of BS might be slower than 9 simple comparisons).

But what precision have your depth queries? Note that you can store a table both for 1-meter resolution and for 1 millimeter too - only 10^7 entries provide O(1) access to any precalculated velocity value

For limited number of pieces it is possible to make long formula (involving integer division) but results perhaps should be slower.

Example for arbitrary three-pieces polyline with border points 2 and 4.5:

f = f0 + 0.2*int(z/2.0)*(z-2.0) + 0.03*int(z/4.5)*(z-4.5)