A visualization of a curve that tries to effectively trap light (slow it down quickly)

light_horn.py

import pyqtgraph as pg
import numpy as np

import sys
from PyQt6.QtCore import Qt
from PyQt6.QtWidgets import QApplication, QHBoxLayout, QLabel, QSizePolicy, QSlider, QSpacerItem, \
    QVBoxLayout, QWidget



class Slider(QWidget):
    def __init__(self, minimum, maximum, parent=None):
        super(Slider, self).__init__(parent=parent)
        self.verticalLayout = QVBoxLayout(self)
        self.label = QLabel(self)
        self.verticalLayout.addWidget(self.label)
        self.horizontalLayout = QHBoxLayout()
        spacerItem = QSpacerItem(0, 20, QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Minimum)
        self.horizontalLayout.addItem(spacerItem)
        self.slider = QSlider(self)
        self.slider.setOrientation(Qt.Orientation.Vertical)
        self.horizontalLayout.addWidget(self.slider)
        spacerItem1 = QSpacerItem(0, 20, QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Minimum)
        self.horizontalLayout.addItem(spacerItem1)
        self.verticalLayout.addLayout(self.horizontalLayout)
        self.resize(self.sizeHint())

        self.minimum = minimum
        self.maximum = maximum
        self.slider.valueChanged.connect(self._setLabelValue)
        self.x = None
        self._setLabelValue(self.slider.value())

    def setValue(self, value : float):
        self.slider.setValue(100 * (value - self.minimum) / (self.maximum - self.minimum))
        # self.slider.setValue(value)
        self._setLabelValue(self.slider.value())


    def _setLabelValue(self, value):
        self.x = self.minimum + (float(value) / (self.slider.maximum() - self.slider.minimum())) * (
        self.maximum - self.minimum)
        self.label.setText("{0:.4g}".format(self.x))


class Widget(QWidget):
    def __init__(self, parent=None):
        super(Widget, self).__init__(parent=parent)
        self.horizontalLayout = QHBoxLayout(self)
        self.w1 = Slider(0, 1)
        self.w1.setValue(.4141414141414141414141414141)
        # self.w1.setValue(50)
        self.horizontalLayout.addWidget(self.w1)

        self.w2 = Slider(-1, 1)
        self.horizontalLayout.addWidget(self.w2)

        self.w3 = Slider(-10, 10)
        self.horizontalLayout.addWidget(self.w3)

        self.w4 = Slider(-10, 10)
        self.horizontalLayout.addWidget(self.w4)

        self.win = pg.GraphicsWindow(title="Basic plotting examples")
        self.horizontalLayout.addWidget(self.win)
        self.p6 = self.win.addPlot(title="My Plot")
        self.horn_upper = self.p6.plot(pen='r')
        self.horn_lower = self.p6.plot(pen='r')
        self.bounce_path = self.p6.plot(pen='w')
        self.update_plot()

        self.w1.slider.valueChanged.connect(self.update_plot)
        self.w2.slider.valueChanged.connect(self.update_plot)
        self.w3.slider.valueChanged.connect(self.update_plot)
        self.w4.slider.valueChanged.connect(self.update_plot)

    def update_plot(self):

        k = self.w1.x
        def horn(x):
            return np.exp(-k*x)
        def horn_bot(x):
            return -horn(x)

        #create numpy arrays
        #make the numbers large to show that the range shows data from 10000 to all the way 0
        xs = np.linspace(-2.,10, 10000)
        horn_upper = horn(xs)
        horn_lower = -horn_upper

        vel = np.array([1,0])
        bounce_points = [np.array([xs[0], 1])]
        def ray(t: float):
            return bounce_points[-1] + t*vel

        epsilon = 1e-10

        # simulate bounces
        for bounce_i in range(10000):
            bot_t = 0
            top_t = 1
            t = bot_t
            horn_func = horn if vel[1] >= 0 else horn_bot

            p = ray(t)
            initial_side = p[1] < horn_func(p[0])
            t = top_t
            # expand search forwards
            for i in range(64):
                p = ray(top_t)
                side = p[1] < horn_func(p[0])
                if side == initial_side:
                    (bot_t, top_t) = (top_t, top_t+2*(top_t - bot_t))
                else:
                    break
            if side == initial_side:
                print(f'never found a crossover point at bounce_i={bounce_i}')
                # bounce_points.append(p)
                break
            # bisect
            for i in range(64):
                t = (bot_t + top_t) / 2
                p = ray(t)
                side = p[1] < horn_func(p[0])
                if side == initial_side:
                    bot_t = t
                else:
                    top_t = t
                if bot_t == top_t:
                    break
            dx = 2*epsilon
            dy = horn_func(p[0] + epsilon) - horn_func(p[0] - epsilon)
            n = np.array([-dy, dx])
            n /= np.sqrt(n.dot(n))  # normal
            # reflect
            vel = vel - 2*(vel.dot(n))*n
            bounce_points.append(p)
            # (bot_t, top_t)

        bounce_xs = [p[0] for p in bounce_points]
        bounce_ys = [p[1] for p in bounce_points]

        self.horn_upper.setData(x=xs, y=horn_upper)
        self.horn_lower.setData(x=xs, y=horn_lower)
        self.bounce_path.setData(x=bounce_xs, y=bounce_ys)


if __name__ == '__main__':
    app = QApplication(sys.argv)
    w = Widget()
    w.show()
    sys.exit(app.exec())

# if __name__ == '__main__':
#     pg.exec()