%
%     SUBROUTINE CUBIC
%
%     PURPOSE
%        SUBROUTINE FINDS THE REAL ROOT, SX, OF THE CUBI% EQUATION,
%        SX**3+SA*SX**2+SBR*SX+SCR=0, FOR WHICH THE FOURTH DEGREE
%        POLYNOMIAL, SX**4/4.+SA*SX**3/3.+SBR*SX**2/2.+SCR*SX,  IS
%        A MINIMUM.
%
%     USAGE
%        CALL CUBIC(SA,SBR,SCR,SX,E)
%
%     DESCRIPTION OF PARAMETERS
%        SA    - SINGLE PRECISION COEFFICIENT OF SX**2. INPUT PARAMETER.
%        SBR   - SINGLE PRECISION COEFFICIENT OF SX. INPUT PARAMETER.
%        SCR   - SINGLE PRECISION CONSTANT TERM, INPUT PARAMETER.
%        SX    - SINGLE PRECISION COMPUTED ROOT. OUTPUT PARAMETER.
%        E     - SINGLE PRECISION COEFFICIENT OF SX**3. HOWEVER ON INPUT
%                TO THE SUBROUTINE THE COEFFIENT OF SX**3. IS ASSUMED
%                TO BE 1.THEREFORE E IS THE NUMBER BY WHICH THE EQUATION
%                WAS DIVIDED BEFORE CALLING CUBIC, IN ORDER TO PRODUCE
%                COEFFICIENTS 1,SA,SBR,SCR.
%
%     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
%        NONE
%
%     METHOD
%        CARDAN'S FORMULAE FOR FINDING THE ROOTS OF A CUBI% EQUATION,
%        AS DESCRIBED IN "THE THEORY OF EQUATIONS" BY J.V. USPENSKY,
%        CHAPTER 5(MCGRAW-HILL BOOK CO, INC., 1948)
%
%     NOTE
%        AFTER THE DESIRED ROOT HAS BEEN OBTAINED, THE SUBROUTINE
%        TESTS TO SEE WHETHER THE FOURTH DEGREE POLYNOMIAL HAS IN
%        FACT BEEN MINIMISED, I.E. WHETHER ITS SECOND DERIVATIVE
%        IS POSITIVE. IF IT IS <=0, AN ERROR MESSAGE IS PRINTED
%        AND EXECUTION TERMINATES.
%

function SX = cubic(SA,SBR,SCR,E)

      EPS = 1; THIRD = 0.3333333333;
      A=SA;
      BR=SBR;
      CR=SCR;
      FE=E;
      S=A/3;
      T=S*A;
      B=0.5*(S*(T/1.5-BR)+CR);
      T=(T-BR)/3;
      C=T^3;
      D=B*B-C;
      if ((abs(D) - EPS) < 0),
         D=(sqrt(D)+abs(B))^(THIRD);
         if (D ~= 0),
             if (B > 0),
                 B=(-1)*D;
             else
                 B=D;
             end;
         end;
         C=T/B;
         D=sqrt(0.75)*(B-C);
         P2=D;
         B=B+C;
         C=(-0.5)*B-S;
         P=C;
         if (B > 0) & (S <= 0),
             Q=C;
             Q2= (-1)*D;
             X = B-S;
             X2= 0;
         else
             Q = B-S;
             Q2= 0;
             X = C;
             X2=(-1)*D;
         end;
         if (((abs(X2)-.001) <0) & ((abs(P2)-.001) <0) & ((abs(Q2)-.001) < 0)),
            R1=(X^4/4+A*X^3/3+BR*X^2/2+CR*X)*FE;
            Q1=(Q^4/4+A*Q^3/3+BR*Q^2/2+CR*Q)*FE;
            P1=(P^4/4+A*P^3/3+BR*P^2/2+CR*P)*FE;
            if ((P1-Q1) <= 0),
               if ((P1-R1) <= 0), X = P; end;
            else
               if ((Q1-R1) <= 0), X = Q; end;
            end;
         else,
            if (P2 == 0), X=P; end;
            if (Q2 == 0), X=Q; end;
         end;
      else
         if (D>=0),
            D=(sqrt(D)+abs(B))^(THIRD);
            if (D ~= 0),
                if (B > 0),
                    B=(-1)*D;
                else
                    B=D;
                end;
            end;
            C=T/B;
            D=sqrt(0.75)*(B-C);
            P2=D;
            B=B+C;
            C=(-0.5)*B-S;
            P=C;
            if (B > 0) & (S <= 0),
                Q=C;
                Q2= (-1)*D;
                X = B-S;
                X2= 0;
            else
                Q = B-S;
                Q2= 0;
                X = C;
                X2=(-1)*D;
            end;
            if (((abs(X2)-.001) <0) & ((abs(P2)-.001) <0) & ((abs(Q2)-.001) < 0)),
               R1=(X^4/4+A*X^3/3+BR*X^2/2+CR*X)*FE;
               Q1=(Q^4/4+A*Q^3/3+BR*Q^2/2+CR*Q)*FE;
               P1=(P^4/4+A*P^3/3+BR*P^2/2+CR*P)*FE;
               if ((P1-Q1) <= 0),
                  if ((P1-R1) <= 0), X = P; end;
               else
                  if ((Q1-R1) <= 0), X = Q; end;
               end;
            else,
               if (P2 == 0), X=P; end;
               if (Q2 == 0), X=Q; end;
            end;
         else

            if ( B ~= 0),
               SQ=sqrt((-1)*D)/abs(B);
               D=atan(SQ)/3;
            else
               D=atan(1)/1.5;
            end;
            if (B < 0),
               B=2*sqrt(T);
            else
               B=(-2)*sqrt(T);
            end;
            C=cos(D)*B;
            T=(-1)*sqrt(0.75)*sin(D)*B-0.5;
            D=-T-C-S;
            C=C-S;
            T=T-S;
            if (abs(C) <= abs(T)),
               X = T;
               T = C;
            else
               X = C;
            end;
            if (abs(D) <= abs(T)),
               Q=T;
               T=D;
            else
               Q=D;
            end;
            R1=(X^4/4+A*X^3/3+BR*X^2/2+CR*X)*FE;
            Q1=(Q^4/4+A*Q^3/3+BR*Q^2/2+CR*Q)*FE;
            P1=(P^4/4+A*P^3/3+BR*P^2/2+CR*P)*FE;
            if ((P1-Q1) <= 0),
                if ((P1-R1) <= 0), X = P; end;
            else
                if ((Q1-R1) <= 0), X = Q; end;
            end;
         end;
      end;
      SD=(3*X*X+2*A*X+BR)*FE;
      SX = X;

%[SD SX]
%pause
