
function [nodes,weights]=define_normal_rule(n,m,a,b,psi,phi)

% Rule direction "x".
abn=r_jacobi(n,0,0); xw=gauss(n,abn); % Gauss-Legendre
x=xw(:,1); w=xw(:,2); x=(a+b)/2+(b-a)*x/2; w=(b-a)*w/2; 

% Rule direction "y".
abm=r_jacobi(m,0,0); yv=gauss(m,abm); % Gauss-Legendre

% Rule on the normal domain
y=yv(:,1); v=yv(:,2);

nodes=[]; weights=[];

for i=1:n
    psi_i=feval(psi,x(i)); phi_i=feval(phi,x(i));
    y_i=(psi_i+phi_i)/2+((phi_i-psi_i)/2)*y; % scaled nodes
    v_i=((phi_i-psi_i)/2)*v; % scaled weights
    
    nodes_add=[x(i)*ones(size(y_i)) y_i]; % rule nodes/weights to add
    nodes=[nodes; nodes_add];
    weights_add=w(i)*v_i;
    weights=[weights; weights_add];
end