function [Fx, Fy, Fz, Tx, Ty, Tz, Fxp, Fyp, Fzp, Txp, Typ, Tzp] = SBT_helical_Lighthill_1(a_lambda, theta, nlength, res, v, 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

    % phin = 1, phim = 2
    phimin=0; phimax=nlength*(2*pi);
    xi=@(phi) sqrt(2 - 2*cos(phi(2)-phi(1))+(phi(2)-phi(1))^2/(tan(theta)^2)); % norm of r(phi2)-r(phi1)) - normalized by R, its dimensionless
    fx12=@(phi) [cos(phi(1))-cos(phi(2)), sin(phi(1))-sin(phi(2)), (phi(1)-phi(2))/tan(theta)]; % (r(phi2)-r(phi1)) - normalized by R, its dimensionless
    px12=@(phi) (rotmatrix(-phi(1))*(fx12(phi)'*fx12(phi))*rotmatrix(phi(2))/xi(phi)^3+rotmatrix(phi(2)-phi(1))/xi(phi)); % non- dimensionalized by 1/R
    Xnorm = @(phi) sqrt(R^2 + (R*phi/tan(theta)+b+h)^2); % has apt dimensions
    Xj = @(phi) [R*cos(phi); R*sin(phi); R*phi/tan(theta)+b+h]; % has apt dimensions
    
    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); % units taken care of here. diagonal terms = 0
        % vphi(i) = phiN, vphi(j) = phiM
    end
    end
    
    vp=zeros(1, 3*nres);
    for i = 1:nres
        A = [v(i);v(nres+i);v(2*nres+i)];
        
        Trans = 3/4*b*(eye(3)/Xnorm(vphi(i)) + Xj(vphi(i))*Xj(vphi(i))'/Xnorm(vphi(i))^3) - b^3/4*(-eye(3)/Xnorm(vphi(i))^3 + 3*Xj(vphi(i))*Xj(vphi(i))'/Xnorm(vphi(i))^5);
%         Trans = 0;
%         Trans = zeros(3,3);
        A = rotmatrix(-vphi(i))*(eye(3)-Trans)*A;
        vp(i) = A(1);
        vp(nres+i) = A(2);
        vp(2*nres+i) = A(3);
    end
    
    fp=mattr\vp';

    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);
    end

    % d = (R cos(phi), R sin(phi), (R phi cot(theta)+b+h))
    Fx = sum(f(1:nres)); Fy = sum(f(nres+1:2*nres)); Fz = sum(f(2*nres+1:3*nres));
    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;
       
       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 = -cons;
%         cons = 1;
       Txp = Txp + t(i)*cons;
       Typ = Typ + t(nres+i)*cons;
       Tzp = Tzp + t(2*nres+i)*cons;
  end

    
end