function [Sigma] = Sigma_Compute(b,c,d,q,hfrac_vec)
% Compute Sigma for "b,c,d" model for q cosine transforms and forecast
% horizon hfrac. Note hfrac_vec is a vector that contains the forecast
% horizons in fractions of the sample size
%
% Calculations for appendix 8.4.3  (copy of UM Fortran program)
%

d0 = d;
c0 = c;

% -- Sigma Appendix Method;
    mkSig_T = floor(1/.0002);        % Numerical intergration grid given appendix
    nhor = ceil(hfrac_vec*mkSig_T);
    Tmax = mkSig_T + max(nhor) + 2;
    v = zeros(Tmax,q+size(hfrac_vec,1));
    vqd = NaN(Tmax,q+size(hfrac_vec,1));  
   
    d = d0;
    if (d < 0.5); d = d+1; end;
    if (d==0.5); d = 0.51; end;
    c = c0;
    if (c0 ==0); c = 0.01; end;
   
    % Compute cosine transforms and weights 
    fvec = pi*(1:1:q);
    tvec = (0.5:1:mkSig_T)'/mkSig_T;
    v(1:mkSig_T,1:q) = (sqrt(2))*cos(tvec*fvec);
    v(1:mkSig_T,q+1:end) = -1;
    for ih = 1:size(hfrac_vec,1);
      v(mkSig_T+1:mkSig_T+nhor(ih),q+ih) = ones(nhor(ih),1)/hfrac_vec(ih);
    end;
    
    % Summation by parts modification
    if (d0 < 0.5);
      vqd(1:Tmax-1,:)=(c/mkSig_T)*v(1:Tmax-1,:)-(v(2:Tmax,:)-v(1:Tmax-1,:));
      vqd(1,:) = vqd(1,:)-v(1,:);
      vqd(Tmax,:) = 0;
      v = mkSig_T*vqd;
    end;
    
    ga_0 = c^(1-2*d)*sqrt(pi)*gamma(-0.5+d)/gamma(d);
    x = (1/mkSig_T)*(1:Tmax-1)';
    cx = c*x;
    bseval = besselk(d-0.5,cx);
    ga = ((2^(1.5-d))*sqrt(pi)/gamma(d))*((x/c).^(d-0.5)).*bseval;
    vqd(1,:)=[ga_0 ga']*v;
    for t = 1:Tmax-2;
        vqd(t+1,:)= [flipud(ga(1:t))' ga_0 ga(1:Tmax-t-1)']*v;
    end;
    vqd(Tmax,:) = [flipud(ga)' ga_0]*v;
    Sig_cd = v'*vqd/(mkSig_T^2);
    Sig_cd = Sig_cd/(2*pi);
    
    % -- b component
    Sig_b = zeros(q+size(hfrac_vec,1),q+size(hfrac_vec,1));
    Sig_b(1:q,1:q) = eye(q);
    for j = 1:size(hfrac_vec,1);
        Sig_b(q+j,q+j:end)=1+1/hfrac_vec(j);
        Sig_b(q+j:end,q+j)=1+1/hfrac_vec(j);
    end;
    
    Sig_b = (b^2)*Sig_b;
    Sig_bcd = Sig_b + Sig_cd;
    
    % Scaled so first element is unity
    Sigma = Sig_bcd/Sig_bcd(1,1);    

end

