subroutine kardan (a,rr,ri,icode,ier)
c*
c***********************************************************************
c*
c*          "Copyright [c] CNES 98 - tous droits reserves"
c*          **********************************************
c*
c*PRO MAGLIB
c*
c*VER 99.03.31 - V 1.0
c*VER 01.06.05 - V 2.0
c*VER 03.01.06 - V 2.1
c*
c*AUT spec. CNES - Mme Soubry
c*AUT spec. CNES - JC KOSIK - janvier 1991
c*AUT port. CISI 
c*
c*ROL Theme : Calculs mathematiques
c*ROL         Calcul des racines reelles ou imaginaires d'une equation
c*ROL         du 3 eme degre.
c*
c*PAR a     (I) : tableau des 4 coefficients de l'equation
c*PAR           : du 3eme degre
c*
c*PAR rr    (O) : tableau des 3 racines reelles
c*PAR ri    (O) : tableau des 3 racines imaginaires
c*
c*PAR icode (O) : nombre de racines reelles trouvees (de 0 a 3)
c*
c*PAR ier   (O) : code de retour
c*
c*NOT ier       : sans objet
c*
c*NOT common    : util
c*
c*INF utilise   : sans objet
c*
c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP
c*HST version 2.0 - 01.06.05 - correction de commentaires de code
c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77
c*
c***********************************************************************
c*
      implicit none
c
c     ---------------------------------
c*FON Declaration identificateur rcs_id
c     ---------------------------------
c
      character rcs_id*100
c
c     --------------------------
c*FON Declaration des parametres
c     --------------------------
c
      double precision a(4)
      double precision rr(3), ri(3)
      integer icode
      integer ier
c
c     ----------------------------------
c*FON Declaration des variables communes
c     ----------------------------------
c
      double precision pi,dpi,rad,deg,pid,xmu,rayt
c
c*COM pi   : constante pi (obtenue a partir de acos(-1.))
c*COM dpi  : constante 2 * pi
c*COM pid  : constante pi / 2
c*COM rad  : facteur de conversion degres  ----> radians
c*COM deg  : facteur de conversion radians ----> degres
c*COM xmu  : constante de gravitation terrestre (km**3/sec**2)
c*COM rayt : rayon equatorial terrestre (km)
c
      common/util/pi,dpi,rad,deg,pid,xmu,rayt
c
c     ---------------------------------
c*FON Declaration des variables locales
c     ---------------------------------
c
      double precision delta
c*LOC delta : discriminant
c
      double precision c1,c2,c3,s1,s2,s3,coef
      double precision fp,fq,rho3,rho,dzeta,fu,fv,fup,fvp
c*LOC Variables de travail intermediaires
c
      SAVE
c
c     ---------------------------------
c*FON Affectation identificateur rcs_id
c     ---------------------------------
c
      data rcs_id /"
     >$Id$"/
c
c     ******************
c     Debut de programme
c     ******************
c
      ier   = 0
      icode = 0
c
c     --------------------------------------------------------------
c*FON Algorithme de calcul des racines si le degre du polynome est 3
c     --------------------------------------------------------------
c
      if (a(1) .ne. 0.0d0) then
c
         c1    = cos(2.d0 * pi / 3.d0)
         c2    = 2.d0 * c1 * c1 - 1.d0
         s1    = sin(2.d0 * pi / 3.d0)
         s2    = 2.d0 * s1 * c1
         coef  = -a(2) / (3.d0 * a(1))
         fp    = (-a(2)**2 + 3.d0 * a(3) * a(1)) /
     >           (3.d0 * a(1)**2)
         fq    = (2.d0 * a(2)**3 - 9.d0 * a(2) *
     >           a(1) * a(3) + 27.d0 * a(1) * a(4) * a(1))
     >           / (27.d0 * a(1)**3)
         delta = fq * fq + 4.d0 * fp**3 / 27.d0
c
c     -------------------------------------------------------------
c*FON Modification de Kardan indiquant le nombre de racines reelles
c     -------------------------------------------------------------
c
         if (delta .lt. 0.d0) then
c
            icode = 3
c
c     ----------------------
c*FON Calculs intermediaires
c     ----------------------
c
            delta = abs(delta)
            rho3  = (fq * fq + delta) / 4.d0
            rho3  = sqrt(rho3)
            rho   = rho3**(1.d0 / 3.d0)
            delta = sqrt(delta)
            dzeta = atan2(delta,-fq)
            dzeta = dzeta / 3.d0
            c3    = cos(dzeta)
            s3    = sin(dzeta)
c
c     ------------------
c*FON Calcul des racines
c     ------------------
c
            rr(1) = 2.d0 * rho * c3 + coef
            ri(1) = 0.d0
            rr(2) = 2.d0 * rho * (c3 * c1 - s3 * s1)
     >              + coef
            ri(2) = 0.d0
            rr(3) = 2.d0 * rho * (c3 * c2 - s3 * s2)
     >              + coef
            ri(3) = 0.d0
c
         else
c
            icode = 1
c
c     ----------------------
c*FON Calculs intermediaires
c     ----------------------
c
            delta = sqrt(delta)
            fu    = (-fq + delta) / 2.d0
            fv    = (-fq - delta) / 2.d0
            fup   = abs(fu)
            fvp   = abs(fv)
            fup   = fup**(1.d0 / 3.d0)
            fvp   = fvp**(1.d0 / 3.d0)
            fu    = sign(fup,fu)
            fv    = sign(fvp,fv)
c
c     ------------------
c*FON Calcul des racines
c     ------------------
c
            rr(1) = fu + fv + coef
            ri(1) = 0.d0
            rr(2) = fu * c1 + fv * c2 + coef
            ri(2) = fu * s1 + fv * s2
            rr(3) = fu * c2 + fv * c1 + coef
            ri(3) = fu * s2 + fv * s1
c
         endif
c
      endif
c
c     ****************
c     Fin de programme
c     ****************
c
      return
      end