function tw=trigauss_kronrod(n,a0,b0)

n0 = 2*n+1;
omega=(b0-a0)/2;

%%%
% Calculation of moments for 2n+1 Gaussian nodes
%%%

% modified Chebyshev moments by recurrence 
z(1)=2*omega;
z(n0+1)=quadgk(@(t)cos(2*n0*acos(sin(t/2)/sin(omega/2))),-omega,omega,'MaxIntervalCount',5000);
temp=(2:2:2*n0-1);
dl=1/4-1./(4*(temp-1));
dc=1/2-1/sin(omega/2)^2-1./(2*(temp.^2-1));
du=1/4+1./(4*(temp+1));
d=4*cos(omega/2)/sin(omega/2)./(temp.^2-1)';
d(n0-1)=d(n0-1)-du(n0-1)*z(n0+1);
z(2:n0)=tridisolve(dl(2:n0-1),dc(1:n0-1),du(1:n0-2),d(1:n0-1));
mom=zeros(1,2*n0+2);
mom(1:2:2*n0+1)=z(1:n0+1);

% normalization of the moments (monic polynomials)
k=(3:length(mom));
mom(3:end)=exp((2-k)*log(2)).*mom(3:end);

% recurrence coeffs of the monic Chebyshev polynomials
abm(:,1)=zeros(2*n0+1,1);
abm(:,2)=0.25*ones(2*n0+1,1); abm(1,2)=pi; abm(2,2)=0.5;

% recurrence coeffs for the monic OPS w.r.t. the weight function
% w(x)=2*sin(omega/2)/sqrt(1-sin^2(omega/2)*x^2)
% by the modified Chebyshev algorithm 
[ab,normsq]=chebyshev(n0,mom,abm);

%%%
% Start of Laurie's algorithm - modifying the recurrence coefficients
%%%

% Offsets so that we can use same notation as in Laurie's paper
% despite indices starting from 1 here
w_off = 2;
i_off = 1;

% Initing first part equal to the gaussian one
a = zeros(2*n+i_off,1);
b = zeros(2*n+i_off,1);
a(0+i_off:floor(3*n/2)+i_off) = ab(0+i_off:floor(3*n/2)+i_off,1);
b(0+i_off:ceil(3*n/2)+i_off) = ab(0+i_off:ceil(3*n/2)+i_off,2);

% Temp vectors
s = zeros(floor(n/2)+w_off,1);
t = zeros(floor(n/2)+w_off,1);


t(0+w_off) = b(n+1+i_off);

for m=0:n-2
    u = 0;
    for k = floor((m+1)/2):-1:0
        l = m-k;
        u = u + (a(k+n+1+i_off) - a(l+i_off))*t(k+w_off) + ...
            b(k+n+1+i_off)*s(k-1+w_off) - ...
            b(l+i_off)*s(k+w_off);
        s(k+w_off) = u;
    end
    temp = t;
    t = s;
    s = temp;
end
for j=floor(n/2):-1:0
    s(j+w_off) = s(j-1+w_off);
end
for m=n-1:2*n-3
    u = 0;
    for k = m+1-n:floor((m-1)/2)
        l = m-k;
        j = n-1-l;
        u = u - (a(k+n+1+i_off)-a(l+i_off))*t(j+w_off) - ...
            b(k+n+1+i_off)*s(j+w_off) + ...
            b(l+i_off)*s(j+1+w_off);
        s(j+w_off) = u;
    end
    if (mod(m,2) == 0)
        k = m/2;
        a(k+n+1+i_off) = a(k+i_off) + ...
            (s(j+w_off) - b(k+n+1+i_off)*s(j+1+w_off))/ ...
            t(j+1+w_off);
    else
        k = (m+1)/2;
        b(k+n+1+i_off) = s(j+w_off)/s(j+1+w_off);
    end
    temp = t;
    t = s;
    s = temp;
end
a(2*n+i_off) = a(n-1+i_off) - b(2*n+i_off)*s(0+w_off)/ ...
    t(0+w_off);

%%% 
% End of Laurie's algorithm for recurrence coefficients
%%%

% Gaussian formula 
xw=gauss(2*n+1,[a';b']');

% angles and weights for the trigonometric gaussian formula
tw(:,1)=2*asin(sin(omega/2)*xw(:,1))+(b0+a0)/2;
tw(:,2)=xw(:,2);

end

% Helper function

function x = tridisolve(a,b,c,d)    
%   TRIDISOLVE  Solve tridiagonal system of equations.
% From Cleve Moler's Matlab suite 
% http://www.mathworks.it/moler/ncmfilelist.html
 
%     x = TRIDISOLVE(a,b,c,d) solves the system of linear equations
%     b(1)*x(1) + c(1)*x(2) = d(1),
%     a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,
%     a(n-1)*x(n-1) + b(n)*x(n) = d(n).
%
%   The algorithm does not use pivoting, so the results might
%   be inaccurate if abs(b) is much smaller than abs(a)+abs(c).
%   More robust, but slower, alternatives with pivoting are:
%     x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
%     x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n)

x = d;
n = length(x);
for j = 1:n-1
   mu = a(j)/b(j);
   b(j+1) = b(j+1) - mu*c(j);
   x(j+1) = x(j+1) - mu*x(j);
end
x(n) = x(n)/b(n);
for j = n-1:-1:1
   x(j) = (x(j)-c(j)*x(j+1))/b(j);
end
end

