#% Function to Construct MONOMIAL-INTEGRATION NODES and WEGHTS
# This function returns the integration nodes using different deterministic
# integration rules:

# Type = Integration rule M1/M2
# N_in = number of shocks  ->  Need to make the routine more general.
#
# Updated: January 20th, 2015
# Julia version: March 21st, 2019
# ==========================================================================

function get_Monomial_nodes(type,par,N_in::Int64)

sig_r = par["sig_r"];
sig_z = par["sig_z"];
sig_g = par["sig_g"];
sig_d = par["sig_d"];

if type=="M1"

    N_nodes = 2*(N_in);

    ER_nodes = sqrt(N_in)*sig_r*([1;-1; 0; 0; 0; 0; 0; 0]);
    EZ_nodes = sqrt(N_in)*sig_z*([0; 0; 1;-1; 0; 0; 0; 0]);
    EG_nodes = sqrt(N_in)*sig_g*([0; 0; 0; 0; 1;-1; 0; 0]);
    ED_nodes = sqrt(N_in)*sig_d*([0; 0; 0; 0; 0; 0; 1;-1]);


    WEIGHT  = (1/N_nodes).*ones(N_nodes);

elseif type=="M2";

    error("Not Ready for this Yet!!!");

    N_nodes = 2*(N_in^2) + 1;

    R = sqrt(2+N_in);
    D = sqrt((2+N_in)/2);

    ER_nodes = sig_r*([0 R -R 0 0 0 0 D D -D -D D D -D -D 0 0 0 0 ]');
    EZ_nodes = sig_z*([0 0 0 R -R 0 0 D -D D -D 0 0 0 0 D D -D -D ]');
    EG_nodes = sig_g*([0 0 0 0 0 R -R 0 0 0 0 D -D D -D D -D D -D ]');

    # Weights
    w1 = 2 / (2+N_in);
    w2 = (4-N_in) / (2*(2+N_in)^2);
    w3 = 1 / (N_in + 2)^2;

    WEIGHT = [w1 w2 w2 w2 w2 w2 w2 w3 w3 w3 w3 w3 w3 w3 w3 w3 w3 w3 w3]';
end

return ER_nodes, EZ_nodes, EG_nodes, ED_nodes, WEIGHT, N_nodes

end
