# Function to Construct GAUSS-HERMITE Integration Nodes
# This function returns the integration nodes using a tensor product rule
#
# Updated (RH): 01/29/2015
# ==========================================================================

function get_GH_nodes(Q_in::Int64,sig_1::Float64,sig_2::Float64,sig_3::Float64,sig_4::Float64)

(quadpoint,weight) = ghpoints(Q_in);    # quadrature nodes and weights

# The following are the points for the case of 4 uncorrelated random process
# with mean zero and variance sig_

er_prime = (sqrt(2)*sig_1*quadpoint);
ez_prime = (sqrt(2)*sig_2*quadpoint);
eg_prime = (sqrt(2)*sig_3*quadpoint);
ed_prime = (sqrt(2)*sig_4*quadpoint);

# Create product rule %

# For two shocks only


ER_nodes = kron(er_prime,ones(Q_in^3));                       # All the er'
EZ_nodes = kron(ones(Q_in),kron(ez_prime,ones(Q_in^2)));    # All the ez'
EG_nodes = kron(ones(Q_in^2),kron(eg_prime,ones(Q_in)));    # All the eg'
ED_nodes = kron(ones(Q_in^3),ed_prime);                       # All the ed'


WEIGHT  = ((1/sqrt(pi))^4)*(kron(weight,kron(weight,kron(weight,weight))));  # weights for the product rule integral

# WEIGHT  = (1/sqrt(pi))*weight;  % weights for the product rule integral

#%%%% FOR 2 SHOCK VERSION %%%%%
#  ER_nodes = kron(er_prime,ones(Q_in,1));                       # All the er'
#  EG_nodes = kron(ones(Q_in,1),eg_prime);    # All the ez'
#
#  WEIGHT  = ((1/sqrt(pi))^2)*(kron(weight,weight));  # weights for the product rule integral
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


#%%%% FOR 3 SHOCK VERSION %%%%%
# ER_nodes = kron(er_prime,ones(Q_in^2,1));                      # All the er'
# EZ_nodes = kron(ones(Q_in,1),kron(ez_prime,ones(Q_in,1)));     # All the ez'
# EG_nodes = kron(ones(Q_in^2,1),eg_prime);                      # All the eg'
#
# WEIGHT  = ((1/sqrt(pi))^3)*(kron(weight,kron(weight,weight)));            # weights for the product rule integral
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

return EG_nodes, ER_nodes, EZ_nodes, ED_nodes, WEIGHT
end
