function track=sybslab(x,ngauss)
% produce a Gauss tracking in 1D slab geometry
% function track=sybslab(x,ngauss)
% (c) 2016 Alain Hebert, Ecole Polytechnique de Montreal
  nreg=size(x,2)-1 ;
  vol=zeros(1,nreg) ;
  for ik=1:nreg
    vol(ik)=x(ik+1)-x(ik) ;
  end
  %[mu0, w0]=lgwt(ngauss, -1.0, 1.0) ;
  if ngauss == 2
    mu0=[0.577350269189626, -0.577350269189626] ; w0=[1.0, 1.0] ;
  elseif ngauss == 4
    mu0=[0.861136311594052, 0.339981043584856] ;
    w0=[0.347854845137454, 0.652145154862546] ;
  elseif ngauss == 6
    mu0=[0.932469514203152, 0.661209386466264,  0.238619186083197] ;
    w0=[0.171324492379170, 0.360761573048139, 0.467913934572691] ;
  elseif ngauss == 8
    mu0=[0.960289856497536, 0.796666477413627, 0.525532409916329, 0.183434642495650] ;
    w0=[0.101228536290377, 0.222381034453374, 0.313706645877887, 0.362683783378362] ;
  elseif ngauss == 10
    mu0=[0.973906528517172, 0.865063366688985, 0.679409568299024, 0.433395394129247, ...
         0.148874338981631] ;
    w0=[0.066671344308688, 0.149451349150581, 0.219086362515982, 0.269266719309996, ...
        0.295524224714753] ;
  elseif ngauss == 14
    mu0=[0.986283808696812, 0.928434883663574, 0.827201315069765, 0.687292904811685, ...
         0.515248636358154, 0.319112368927890, 0.108054948707344] ;
    w0=[0.035119460331752, 0.080158087159760, 0.121518570687903, 0.157203167158194, ...
        0.185538397477938, 0.205198463721296, 0.215263853463158] ;
  elseif ngauss == 18
    mu0=[0.991565168420931, 0.955823949571398, 0.892602466497556, 0.803704958972523, ...
         0.691687043060353, 0.559770831073948, 0.411751161462843, 0.251886225691506, ...
         0.084775013041735] ;
    w0=[0.021616013526483, 0.049714548894969, 0.076425730254889, 0.100942044106287, ...
        0.122555206711478, 0.140642914670651, 0.154684675126265, 0.164276483745833, ...
        0.169142382963144] ;
  else
    error('invalid number of Gauss-Legendre points.')
  end
  eta=zeros(1,ngauss) ; xi=eta ; weight=eta ;
  for igauss=1:ngauss/2
    eta(igauss)=mu0(igauss) ;
    xi(igauss)=sqrt(1.0-mu0(igauss)^2);
    weight(igauss)=w0(igauss)*mu0(igauss)  ;
    eta(ngauss/2+igauss)=mu0(igauss) ;
    xi(ngauss/2+igauss)=-sqrt(1.0-mu0(igauss)^2);
    weight(ngauss/2+igauss)=weight(igauss)  ;
  end
  nbtrack=7+nreg+2*ngauss+ngauss*(5+2*nreg) ;
  track=zeros(1,nbtrack) ;
  nbtrack=7+nreg+2*ngauss ;
  track(1:nbtrack)=[2, nreg, ngauss, ngauss, 1.0, 1.0, 1.0, vol, eta, xi] ;
  for iangle=1:ngauss
    track(nbtrack+1:nbtrack+5+2*nreg)=[iangle, 1, 2, weight(iangle), nreg, ...
       1:nreg, vol/eta(iangle)] ;
    nbtrack=nbtrack+5+2*nreg ;
  end