function [num_gauss,first_gauss_coord,second_gauss_coord,gauss_weights]=GaussQuad(order)
%GAUSSQUAD - get Gauss point for quad elements
%
%   This routine computes gaussian points on a reference qud elements
%   The number of gaussian points returned depends on the order of integration ('order'). The order of 
%   integration can be computed from the polynomial degree p that needs to be integrated. The formula is: 
%   order = (p+1) /2 
%   order=1, num_gauss=1. Can integrate polynomials of degree 0 to 1
%
%   Usage:
%      [num_gauss,first_gauss_coord,second_gauss_coord,gauss_weights]=GaussQuad(order)


if order==1,
	num_gauss=1;
	gauss_weights=4;
	first_gauss_coord= 0;
	second_gauss_coord= 0;

%order=2, num_gauss=2. Can integrate polynomials of degree 0 to 3
elseif order==2,
	num_gauss=4;
	gauss_weights=[ 1, 1, 1, 1];
	first_gauss_coord=[ -0.57735026918962, 0.57735026918962, -0.57735026918962, 0.57735026918962];
	second_gauss_coord=[0.57735026918962, 0.57735026918962, -0.57735026918962, -0.57735026918962];

%order=3, num_gauss=3. Can integrate polynomials of degree 0 to 5
elseif order==3,
	num_gauss=9;
	gauss_weights=[ 0.30864197530864, 0.49382716049383, 0.30864197530864, 0.49382716049383, 0.79012345679012, 0.49382716049383, 0.30864197530864, 0.49382716049383, 0.30864197530864];
	first_gauss_coord=[ -0.77459666924148, 0, 0.77459666924148, -0.77459666924148, 0, 0.77459666924148, -0.77459666924148, 0, 0.77459666924148];
	second_gauss_coord=[ -0.77459666924148, -0.77459666924148, -0.77459666924148, 0, 0, 0, 0.77459666924148, 0.77459666924148, 0.77459666924148];

else
	error('GaussQuad error message: order not supported yet');
end
