function [ theta_cut, thetap_pl, thetap_mn ] = phase_d_energy( En, Ueff, theta )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

j=1;
for i=1:length(theta)
    if  (En-Ueff(i))>0
        theta_cut(j)=theta(i);
        thetap_pl(j)=+sqrt(2*( En - Ueff(i) ));
        thetap_mn(j)=-sqrt(2*( En - Ueff(i) ));
        j=j+1;
    end
end

