function[alpha]= Split_Integration(N_alpha,kappa)

% Gauss quadrature nodes and weights on [-1,1]:
smallInterval=2^(-N_alpha);

if kappa ==0

    [x1,w1] = lgwt(N_alpha,-0.5,-smallInterval);
    x1 = x1(end:-1:1);
    w1 = w1(end:-1:1);

    alpha.nodes = [x1;-smallInterval;0;smallInterval;-x1(end:-1:1)];
    alpha.weights =[w1;smallInterval/2;smallInterval; smallInterval/2;w1(end:-1:1)];

elseif kappa ==0.5 || kappa ==-0.5

    alpha.nodes = zeros(N_alpha+4,1);
    alpha.weights = alpha.nodes;

    [x1,w1] = lgwt(N_alpha,-0.5+smallInterval,0.5-smallInterval);

    x1 = x1(end:-1:1);
    w1 = w1(end:-1:1);

    alpha.nodes(1:2) = [-1/2;-1/2+smallInterval];
    alpha.nodes(end-1:end) = [1/2-smallInterval;1/2];

    alpha.nodes(3:end-2)=x1;

    alpha.weights(3:end-2) = w1;
    alpha.weights(1:2) = smallInterval/2;
    alpha.weights(end-1:end) = smallInterval/2;

else

    alpha.nodes = zeros(3*N_alpha+6,1);
    alpha.weights = alpha.nodes;

    [x1,w1] = lgwt(N_alpha,-1,1);
    x1 = x1(end:-1:1);
    w1 = w1(end:-1:1);

    intervals=[-0.5,-kappa-smallInterval;-kappa+smallInterval,kappa-smallInterval;kappa+smallInterval,0.5];

    for int = 0:2
        index = (N_alpha+3)*int+1:(N_alpha+3)*int+N_alpha;
        [alpha.nodes(index),alpha.weights(index)] = linearMap(x1,w1,intervals(int+1,1),intervals(int+1,2));
    end

    for int = [-1,1]
        index=(N_alpha+3)*0.5*(int+1)+N_alpha+(1:3);
        alpha.nodes(index) = kappa*int+[-smallInterval;0;smallInterval];
        alpha.weights(index) = [smallInterval/2;smallInterval;smallInterval/2];
    end  
end
end

function [x,w]=linearMap(y,w0,a,b)
% Linear map from[-1,1] to [a,b]
x = (a*(1-y)+b*(1+y))/2;      
% Compute the weights
w = 0.5*(b-a)*w0;
end
