function [Fx, Fy, Fz, Tx, Ty, Tz, Fxp, Fyp, Fzp, Txp, Typ, Tzp, Q] = SBT_helical_Lighthill_2(a_lambda, theta, nlength, res, omega, R,  b, h)
% Note that these output forces are normalized by 4 pi mu!  
    
    function m = rotmatrix(angle)
        m=[cos(angle), -sin(angle), 0; sin(angle), cos(angle), 0; 0, 0, 1];
    end

    phimin=0; phimax=nlength*(2*pi);
    xi=@(phi) sqrt(2 - 2*cos(phi(2)-phi(1))+(phi(2)-phi(1))^2/(tan(theta)^2));
    fx12=@(phi) [cos(phi(1))-cos(phi(2)), sin(phi(1))-sin(phi(2)), (phi(1)-phi(2))/tan(theta)];
    px12=@(phi) (rotmatrix(-phi(1))*(fx12(phi)'*fx12(phi))*rotmatrix(phi(2))/xi(phi)^3+rotmatrix(phi(2)-phi(1))/xi(phi)); 
    Xnorm = @(phi) sqrt(R^2 + (R*phi/tan(theta)+b+h)^2);
    Xj = @(phi) [R*cos(phi); R*sin(phi); R*phi/tan(theta)+b+h];
    
    nres=floor(nlength*res+.5);
    dphi=(phimax-phimin)/nres;
    vphi=phimin+.5*dphi:dphi:phimax-.5*dphi;
    mattr=zeros(3*nres, 3*nres);
    cutofflen=a_lambda*pi*exp(0.5)*cos(theta);
    cl=log(.5*dphi/cutofflen);  % Approximation!

    for i=1:nres

    mattr(i,i)=1+cl;
    mattr([nres+i, 2*nres+i], [nres+i, 2*nres+i])=[cos(theta)^2+(1+sin(theta)^2)*cl, (cl-1)*sin(theta)*cos(theta); (cl-1)*cos(theta)*sin(theta), sin(theta)^2+(1+cos(theta)^2)*cl];
    end

    for i=1:nres
    for j=1:nres
        if (abs(vphi(i)-vphi(j))<cutofflen+.5*dphi)
        continue;
        end
        mattr([i, nres+i, 2*nres+i], [j, nres+j, 2*nres+j])=px12([vphi(i), vphi(j)])*dphi*.5/sin(theta); % Order of passing phi?
        % vphi(i) = phiN, vphi(j) = phiM
    end
    end
    
    v=zeros(1, 3*nres);
    vp=zeros(1, 3*nres);
    Q=zeros(3*nres,1);
    
    % d = (R cos(phi), R sin(phi), (R phi cot(theta)+b+h))
    % v = omega cross r(position vector)
    for i = 1:nres
        v(i) = (R*vphi(i)/tan(theta)+b+h)*omega(2) - R*sin(vphi(i))*omega(3);
        v(nres+i) = R*cos(vphi(i))*omega(3) - (R*vphi(i)/tan(theta)+b+h)*omega(1);
        v(2*nres+i) = R*sin(vphi(i))*omega(1) - R*cos(vphi(i))*omega(2);
        
        Rot = (1- (b/Xnorm(vphi(i)))^3);
%         Rot = 1;
        
        A = [v(i);v(nres+i);v(2*nres+i)];
        A = rotmatrix(-vphi(i))*Rot*A;
        vp(i) = A(1);
        vp(nres+i) = A(2);
        vp(2*nres+i) = A(3);
        
        Q(i) = -R*sin(vphi(i));
        Q(nres+i) =  R*cos(vphi(i));
        Q(2*nres+i) = 0;
%         Q(i) = 0; Q(nres+i) = 0; Q(2*nres+i) = 0;
        
        B = [Q(i);Q(nres+i);Q(2*nres+i)];
        B = rotmatrix(-vphi(i))*B;
        Q(i) = B(1);
        Q(nres+i) = B(2);
        Q(2*nres+i) = B(3);
        
    end
    
    fp=mattr\vp';
    Q = mattr\Q;

    f = zeros(3*nres,1);
    t = zeros(3*nres,1);
    for i = 1:nres
        A = [fp(i);fp(nres+i);fp(2*nres+i)];
        A = rotmatrix(vphi(i))*A;
        f(i) = A(1);
        f(nres+i) = A(2);
        f(2*nres+i) = A(3);
        
        B = [Q(i);Q(nres+i);Q(2*nres+i)];
        B = rotmatrix(vphi(i))*B;
        Q(i) = B(1);
        Q(nres+i) = B(2);
        Q(2*nres+i) = B(3);
    end

    Qp = zeros(3,1);
    Fx = sum(f(1:nres)); Fy = sum(f(nres+1:2*nres)); Fz = sum(f(2*nres+1:3*nres));
    Qp(1) = sum(Q(1:nres)); Qp(2) = sum(Q(nres+1:2*nres)); Qp(3) = sum(Q(2*nres+1:3*nres));
    Q = Qp; clear Qp;
    
    Tx = 0; Ty = 0; Tz = 0;
    for i = 1:nres

        t(i) = (R*sin(vphi(i))*f(2*nres+i) - (R*vphi(i)/tan(theta)+b+h)*f(nres+i));
        t(nres+i) = ((R*vphi(i)/tan(theta)+b+h)*f(i) - R*cos(vphi(i))*f(2*nres+i));
        t(2*nres+i) = (R*cos(vphi(i))*f(nres+i) - R*sin(vphi(i))*f(i));

    end

    Tx = sum(t(1:nres)); Ty = sum(t(nres+1:2*nres)); Tz = sum(t(2*nres+1:3*nres));
    
   Fxp = 0; Fyp = 0; Fzp = 0; Txp = 0; Typ = 0; Tzp = 0; 
   for i = 1:nres
       cr = -3/2*b/Xnorm(vphi(i)) + 1/2*(b/Xnorm(vphi(i)))^3;
       ct = -3/4*b/Xnorm(vphi(i)) - 1/4*(b/Xnorm(vphi(i)))^3;

%        cr = 0; ct = 0;
%        cr = -cr; ct = -ct;
       D = Xj(vphi(i))*Xj(vphi(i))'/Xnorm(vphi(i))^2;
       N = D*[f(i);f(nres+i);f(2*nres+i)];
       Fxp = Fxp + f(i)*(1+ct) + (cr-ct)*N(1);
       Fyp = Fyp + f(nres+i)*(1+ct) + (cr-ct)*N(2);
       Fzp = Fzp + f(2*nres+i)*(1+ct) + (cr-ct)*N(3);
   end
   
  for i = 1:nres
       cons = 1-(b/Xnorm(vphi(i)))^3;
%         cons = 1;
%        cons = -cons;
       Txp = Txp + t(i)*cons;
       Typ = Typ + t(nres+i)*cons;
       Tzp = Tzp + t(2*nres+i)*cons;
   end

end