public ComplexMatrix getLocalPropagationMatrix(int index, double energy) {
        MatterInterval inl = intervals.get(index);
        MatterInterval next = intervals.get(index + 1);
        Complex p;
        // Double.POSITIVE_INFINITY is handled here and it assigns to
        // real part
        // of Complex if appears
        p = Complex.sqrt((energy - next.getY()) / (energy - inl.getY()) * inl.getMass() / next.getMass());

        Complex kcp = p.add(1); // 1 + p
        Complex kcm = p.reverseSign().add(1); // 1 - p
        Complex kd = next.getWaveNumber(energy).times(next.getWidth());

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

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