function [REALrts, ALLrts] = cardanos(p,q)
% CARDANOS(p,q) root finding algorithm for depressed cubic polynomial with real valued p and q
% reduced form of CardanRoots for surface renewal vectors p and q to determine ramp Amplitudes
% polynomial should be of form A^3 + p*A + q = 0 with p and q real valued (no complex values)
% RETURNS only positive real valued solutions, with complex and negative solutions replaced with NaN
% refs: <a href="matlab:web('https://en.wikipedia.org/wiki/Cubic_function#Cardano.27s_method','-browser')">Wiki</a>
% SEE ALSO: cardanroots.m

% c3*x^3+c2*x^2+c1*x+c0=0 attempts at full form
%{
% % <a href="matlab:web('https://www.mathworks.com/matlabcentral/newsreader/view_thread/165013?requestedDomain=www.mathworks.com','-browser')">Source Code</a>
% &prev cardanos(b,c,d)

%  bs = b.^2;
%   p = c-bs./3;
%   q = ((-2/27)*bs+(1/3)*c).*b - d;
%    Q = p./3;
%    R = q./2;
%   D = Q.^3 + R.^2;


% k = find(D>=0);
%    S = zeros(length(D),1);
%    T = zeros(length(D),1);
% S(k) = nthroot(R(k)+sqrt(D(k)),3);
% T(k) = nthroot(R(k)-sqrt(D(k)),3);
% 
% root      = zeros(length(D),1);
% root(k,1) = -b./3+S(k)+T(k);
% root(k,2) = -b./3+(1/2).*(S(k)+T(k))+(1/2).*1i.*sqrt(3).*(S(k)-T(k));
% root(k,3) = -b./3+(1/2).*(S(k)+T(k))-(1/2).*1i.*sqrt(3).*(S(k)-T(k));
% 
% k=find(D<0);
%    phi    = zeros(length(D),3);
%    phi(k) = acos(R(k)./sqrt(-Q(k).^3));
% root(k,1) = 2.*sqrt(-Q(k)).*cos(phi(k)./3)-b(k)./3;
% root(k,2) = 2.*sqrt(-Q(k)).*cos((phi(k)+2*pi)./3)-b(k)./3;
% root(k,3) = 2.*sqrt(-Q(k)).*cos((phi(k)+4*pi)./3)-b(k)./3;
%}

D = q.^2 + (4/27)*p.^3;         % the discriminant

                                     %deltanull = abs(delta)<tol; % = 0
 Dneg = D<0;
 Dpos = ~Dneg;                       %~(deltanull | deltaneg);

    n = size(D,1);
 rts = zeros(n, 3);
    %root(Dneg,:) = CardanNeg(q(Dneg), D(Dneg));
        a = -q(Dneg);
        b = sqrt(-D(Dneg));
        r2 = a.^2-D(Dneg);
        rho = (4^(1/3))*exp(log(r2)/6);
        theta = atan2(b,a)/3;
        a = rho.*cos(theta);
        b = rho.*sin(theta);
        S1 = a;
        x = (-0.5)*a;
        y = (sqrt(3)/2)*b;
        S2 = x-y;
        S3 = x+y;
        rts(Dneg,1:3) = [S1 S2 S3];
    %rts(Dpos,:) = CardanPos(q(Dneg), D(Dneg));
        E = sqrt(D(Dpos));
        u3 = (-q(Dpos)+E)/2;
        v3 = (-q(Dpos)-E)/2;
        u = sign(u3).*exp(log(abs(u3))/3);  % Cubic roots of u3 and v3
        v = sign(v3).*exp(log(abs(v3))/3);
        S1 = u+v;
        j = complex(-0.5,sqrt(3)/2);        % Complex solutions
        j2 = complex(-0.5,-sqrt(3)/2);
        S2 = j*u+j2*v;
        S3 = conj(S2);
        rts(Dpos,1:3) = [S1 S2 S3];

%roots = bsxfun(@minus, roots, b/3);        %not necessary for reduced form

ALLrts = rts;

rts(imag(rts)~=0) = NaN;
%rts(rts<=0) = NaN;

%REALPOSrts = rts;
REALrts = rts;

end