function [Reslim, Rerlim, fcub, sigwfrlim, frough, psigmasmooth, psigmarough, a2, a1, a0, nstar] = ...
    getfrictioncoeffs(roughness, lower, upper, nterms)

% Polynomial  coefficients for unsteady friction weighting terms
a2all = [0.0196466 0.0433595 0.0314749 0.0247278 0.0126470 0.0043957 0.0008864 0.0001254 -0.0000727 -0.0001712 0.0007275 0.0001389];
a1all = [0.150494 -0.094614 -0.036212 -0.066953 -0.037386 -0.015141 -0.002871 -0.000666 0.0010975 -0.0086700 0.011337 0.0015409];
a0all = [-0.25198 0.26256 0.79916 1.33203 1.79049 2.25970 2.73528 3.19317 3.75574 3.98899 4.47219 5.39664];
nstarall = [20 60 360*9.^[0:9]];
if nterms > 12
    error ('Number of friction weighting terms can''t exceed 12')
end
a2 = a2all(1:nterms);
a1 = a1all(1:nterms);
a0 = a0all(1:nterms);
nstar = nstarall(1:nterms);

% Polynomial coefficients for sigma
psigmasmooth = [0.1309   -0.1119];
psigmarough = [-0.7025   -2.9936   -5.4931   -3.5609];

% Find upper limit of Re for fully smooth flow
Rerootf=lower*sqrt(2)./roughness;
factor = 2*log10(Rerootf*2) - 0.8;
if factor <= 0 % out of range. 
    emsg = ['Lower bound of transition region is too low for this roughness. ',...
        'You may need to increase the lower transition limit, or reduce the relative roughness.'];
    error (emsg)
end
rootf = 1./(2*factor);
fslim = rootf^2;
Reslim = Rerootf ./ rootf

dlfdres = -1/(log(10)/(8*sqrt(fslim))+0.5);

% Transition - rough point. 
% find f for fully rough flow (independent of Re)
rootfrough = 1./(2*(1.74 - 2*log10(2*roughness))); 
frough = (rootfrough.^2);
frlim = frough;
% Find lower limit of Re for rough flow
Rerlim = upper./(roughness .* rootfrough/sqrt(2)); 
dlfdrer = 0;

% Find sigmas at this point.
% sigmacw
sigcwrlim = polyval(psigmarough,log10(frlim));

% sigmacf
sigcfrlim = log10(0.065*Rerlim/2*sqrt(frlim/2));

% find sigmawf from these.
sigwfrlim = (sigcfrlim - sigcwrlim)

% convert to cubic coefficients
dxdlr = log10(Rerlim) - log10(Reslim);
fcub=ccoeff(log10(fslim), dlfdres*dxdlr, log10(frlim), dlfdrer*dxdlr);



function u = ccoeff(x, dx, y, dy)
%note: terms have been flipped to match normal Matlab polynomial convention
u(3) = dx;
u(4) = x;
u(1) = 2*u(4) + u(3) - 2*y + dy;
u(2) = y - u(4) - u(3) - u(1);