Quadrature_1D.m

function [zgp, wgp] = Quadrature_1D(ngaus)
%
% [zgp, wgp] = Quadrature_1D(ngaus)
% Integration points (zgp) and weights (wgp) for a reference element [-1,1].
% ngaus is the number of integration points to be used.

if ngaus == 1
    zgp = 0;
    wgp = 2;
elseif ngaus == 2
   pos1 = 1/sqrt(3); 
   zgp = [-pos1; pos1];
   wgp = [ 1 1 ]; 
elseif ngaus == 3
   pos1 = sqrt(3/5); 
   zgp = [-pos1; 0; pos1];
   wgp = [ 5/9   8/9   5/9 ];   
elseif ngaus == 4
   pos1 = sqrt(525+70*sqrt(30))/35;
   pos2 = sqrt(525-70*sqrt(30))/35;
   zgp = [-pos1; -pos2; pos2; pos1];       
   w1 = sqrt(30)*(3*sqrt(30)-5)/180;
   w2 = sqrt(30)*(3*sqrt(30)+5)/180;
   wgp = [w1   w2   w2   w1];
elseif ngaus == 5
     r70 = sqrt(70);
     pos1 = sqrt(245+14*r70)/21;
     pos2 = sqrt(245-14*r70)/21;
     zgp = [-pos1; - pos2; 0; pos2; pos1];
     w1 = (7+5*r70)*3*r70/(100*(35+2*r70)); 
     w2 = -(-7+5*r70)*3*r70/(100*(-35+2*r70));
     w0 = 128/225;
     wgp = [w1,w2,w0,w2,w1];
elseif ngaus == 6
    zgp = [ 0.23861918608319690863050172168066
              0.66120938646626451366139959501991
              0.93246951420315202781230155449406
             -0.23861918608319690863050172168066
             -0.66120938646626451366139959501991
             -0.93246951420315202781230155449406  ]; 
    wgp = [ 0.46791393457269104738987034398801, ...
            0.36076157304813860756983351383812, ...
            0.17132449237917034504029614217260, ...
            0.46791393457269104738987034398891, ...
            0.36076157304813860756983351383816, ...
            0.17132449237917034504029614217271 ];
elseif ngaus == 7
    zgp = [ 0.
              0.40584515137739716690660641207692
              0.74153118559939443986386477328078
              0.94910791234275852452618968404784
             -0.40584515137739716690660641207692
             -0.74153118559939443986386477328078
             -0.94910791234275852452618968404784 ]; 
    wgp = [ 0.4179591836734693877551020408166, ...
            0.38183005050511894495036977548841, ...
            0.27970539148927666790146777142377, ...
            0.12948496616886969327061143267904, ...
            0.38183005050511894495036977548964, ...
            0.27970539148927666790146777142336, ...
            0.12948496616886969327061143267912 ];    
else
   error('Nonavailable 1D quadrature');   
end