2

Hi I have two numpy arrays (in this case representing depth and percentage depth dose data) as follows:

depth = np.array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ,
                   1.1,  1.2,  1.3,  1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2. ,  2.2,
                   2.4,  2.6,  2.8,  3. ,  3.5,  4. ,  4.5,  5. ,  5.5])
pdd = np.array([  80.40649399,   80.35692155,   81.94323956,   83.78981286,
                  85.58681373,   87.47056637,   89.39149833,   91.33721651,
                  93.35729334,   95.25343909,   97.06283306,   98.53761309,
                  99.56624117,  100.        ,   99.62820672,   98.47564754,
                  96.33163961,   93.12182427,   89.0940637 ,   83.82699219,
                  77.75436857,   63.15528566,   46.62287768,   29.9665386 ,
                  16.11104226,    6.92774817,    0.69401413,    0.58247614,
                   0.55768992,    0.53290371,    0.5205106 ])

which when plotted give the following curve:

Percentage Depth Dose

I need to find the depth at which the pdd falls to a given value (initially 50%). I have tried slicing the arrays at the point where the pdd reaches 100% as I'm only interested in the points after this.

Unfortunately np.interp only appears to work where both x and y values are incresing.

Could anyone suggest where I should go next?

ali_m
  • 71,714
  • 23
  • 223
  • 298
Trigfa
  • 65
  • 1
  • 7

3 Answers3

1

If I understand you correctly, you want to interpolate the function depth = f(pdd) at pdd = 50.0. For the purposes of the interpolation, it might help for you to think of pdd as corresponding to your "x" values, and depth as corresponding to your "y" values.

You can use np.argsort to sort your "x" and "y" by ascending order of "x" (i.e. ascending pdd), then use np.interp as usual:

# `idx` is an an array of integer indices that sorts `pdd` in ascending order
idx = np.argsort(pdd)

depth_itp = np.interp([50.0], pdd[idx], depth[idx])

plt.plot(depth, pdd)
plt.plot(depth_itp, 50, 'xr', ms=20, mew=2)

enter image description here

ali_m
  • 71,714
  • 23
  • 223
  • 298
  • Thanks very much for your help. It works well. Combining this with slicing for values greater than pdd=100% means that it will also work for other values greater than 80% – Trigfa Jan 18 '16 at 15:53
0

This isn't really a programming solution, but it's how you can find the depth. I'm taking the liberty of renaming your variables, so x(i) = depth(i) and y(i) = pdd(i).

In a given interval [x(i),x(i+1)], your linear interpolant is

p_1(X) = y(i) + (X - x(i))*(y(i+1) - y(i))/(x(i+1) - x(i))

You want to find X such that p_1(X) = 50. First find i such that x(i)>50 and x(i+1), then the above equation can be rearranged to give

X = x(i) + (50 - y(i))*((x(i+1) - x(i))/(y(i+1) - y(i)))

For your data (with MATLAB; sorry, no python code) I make it approximately 2.359. This can then be verified with np.interp(X, depth, pdd)

Steve
  • 1,579
  • 10
  • 23
0

There are several methods to carry out interpolation. For your case, you are basically looking for the depth at 50% which is not available in your data. The simplest interpolation is the linear case. I'm using numerical recipes library in C++ for acquiring the interpolated value via several techniques, therefore,

Linear Interpolation: see page 117

interpolated value depth(50%): 2.35915

Polynomial Interpolation: see page 117

interpolated value depth(50%): 2.36017

Cubic Spline Interpolation: see page 120

interpolated value depth(50%): 2.19401

Rational Function Interpolation: see page 124

interpolated value depth(50%): 2.35986

CroCo
  • 5,531
  • 9
  • 56
  • 88