eic-opticks/ana/split_interpolation.py

0001 #!/usr/bin/env python
0002 """
0003
0004 * https://stackoverflow.com/questions/48028766/get-x-values-corresponding-to-y-value-on-non-monotonic-curves
0005
0006 """
0007
0008 import numpy as np
0009 from opticks.ana.mlib import GMaterialLib
0010 mlib = GMaterialLib()
0011
0012 import matplotlib.pyplot as plt
0013 from scipy.interpolate import interp1d
0014
0015 class SplitInterpolation(object):
0016     def __init__(self, x, y ):l
0017
0018
0019
0020
0021 if __name__ == '__main__':
0022
0023     #name = "LS"
0024     name = "Water"
0025     ri = mlib("%s.RINDEX" % name).copy()
0026
0027     # energy(eV) and refractive index in ascending energy
0028
0029     #x,y = mlib.ev[::-1], ri[::-1]
0030     #x,y = mlib.nm, ri
0031     x,y = ri, mlib.nm,
0032
0033
0034     order = np.argsort(x)
0035     xs, ys = x[order], y[order]
0036
0037     fig, axs = plt.subplots(2)
0038     axs[0].plot( x, y )
0039     axs[1].plot(xs, ys )
0040     fig.show()
0041
0042
0043     # compute indices of points where y changes direction
0044     ydir = np.sign(np.diff(ys))
0045     yturn = 1 + np.where(np.diff(ydir) != 0)[0]
0046
0047     # find groups of x and y within which y is monotonic
0048     xgrp = np.split(xs, yturn)
0049     ygrp = np.split(ys, yturn)
0050
0051
0052 if 0:
0053     interps = [interp1d(y, x, bounds_error=False) for y, x in zip(ygrp, xgrp)]
0054
0055     # interpolate all y values
0056     yval = 100
0057     xvals = np.array([interp(yval) for interp in interps])
0058
0059     print(xvals)
0060     # array([          nan,  122.96996037,  207.62395521,           nan])
0061
0062
0063
0064
0065