public ComplexMatrix getLocalPropagationMatrix(int index, double energy) {

        MatterInterval inl = intervals.get(index);
        MatterInterval next = intervals.get(index + 1);

        Complex k_ratio;
        // Double.POSITIVE_INFINITY is handled here and it assigns to real part of Complex if appears
        k_ratio = Complex.sqrt((energy - next.getY()) / (energy - inl.getY()) * inl.getMass() / next.getMass());

        Complex kcp = k_ratio.add(1); // 1 + km_ratio
        Complex kcm = k_ratio.reverseSign().add(1); // 1 - km_ratio
        Complex kd = inl.getWaveNumber(energy).times(inl.getWidth());

        if (VERBOSE) {
            System.out.println("k_ratio = " + k_ratio);
            System.out.println("kcp = " + kcp);
            System.out.println("kcm = " + kcm);
            System.out.println("kd = " + kd);
        }

        Complex i = new Complex(0, 1);
        Complex ni = i.reverseSign();
        try {
            return new ComplexMatrix(new Complex[][]{
                    {kcp.times(Complex.exp(ni.times(kd))), kcm.times(Complex.exp(ni.times(kd)))},
                    {kcm.times(Complex.exp(i.times(kd))), kcp.times(Complex.exp(i.times(kd)))}}).multiplyEach(0.5);
        } catch (Exception e) {
            e.printStackTrace();
            return null;
        }
    }