PFV4 quadature rules

From CFD-Wiki

function [ xa,ya,wt,nq ] = GQuad2( ord )
%GQuad2
% Returns arrays of integration points and weights for specified degree
% Usage:
%   [xa,ya,wt,nq] = GQuad2( ord )
% where
%   ord - highest degree of term integrated exactly on interval [-1,1]
% 
  if  ord<1,     error('integration order < 1');
  elseif ord>15, error('max quad data order = 15');
  end
  nt = fix(ord/2)+1;
  
    switch nt;       
       case {1, 2}  % degree 3
   Aq=[-0.5773502691896257, 0.5773502691896257];
   Hq=[ 1.0,                1.0];  
        case 3    % degree 5
   Aq=[-0.7745966692414833, 0.000000000000000, 0.7745966692414833];
   Hq=[ 0.5555555555555558, 0.8888888888888888,0.5555555555555558];
        case 4    % degree 7
   Aq=[-0.8611363115940526,-0.3399810435848563, 0.3399810435848563, ...
        0.8611363115940526];
   Hq=[ 0.3478548451374538, 0.6521451548625461, 0.6521451548625461, ...
        0.3478548451374538];
        case 5    % degree 9
   Aq=[-0.9061798459386640,-0.5384693101056831, 0.0000000000000000, ...
        0.5384693101056831, 0.9061798459386640];
   Hq=[ .23692688505618909, .47862867049936647, .56888888888888889, ...
        .47862867049936647, .23692688505618909];
        case 6    % degree 11
   Aq=[-.932469514203152,-.661209386466265,-.238619186083197, ...
        .238619186083197, .661209386466265, .932469514203152];
   Hq=[ .171324492379170, .360761573048139, .467913934572691, ...
        .467913934572691, .360761573048139, .171324492379170];
        case 7   % degree 13
   Aq=[-0.9491079123427585,-0.7415311855993945,-0.4058451513773972, ...
        0.0000000000000000, 0.4058451513773972, 0.7415311855993945, ...
        0.9491079123427585];
   Hq=[ 0.1294849661688699, 0.2797053914892767, 0.3818300505051190, ...
        0.4179591836734694, 0.3818300505051190, 0.2797053914892767, ...
        0.1294849661688699];
        case 8   % degree 15
   Aq=[-0.9602898564975363,-0.7966664774136267,-0.5255324099163290, ...
       -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, ...
        0.7966664774136267, 0.9602898564975363];
   Hq=[ 0.1012285362903761, 0.2223810344533745, 0.3137066458778873, ...
        0.3626837833783620, 0.3626837833783620, 0.3137066458778873, ...
        0.2223810344533745, 0.1012285362903761];

    end
    
    nt=size(Aq,2);
    xa = repmat(Aq,nt,1);
    ya = xa';
    wt = repmat(Hq,nt,1);
    wt = wt.*wt';
    nq = nt*nt;
    return;