      subroutine cubic(coeff,xreal,ximage)
      implicit none
      double precision coeff(4),xreal(3),ximage(3),a,b,c,r,q,disc
      double precision r2,q3,theta,arg,third,pi,barg
      third=1.0d0/3.0d0
      pi=dacos(-1.0d0)
      a=coeff(2)/coeff(1)
      b=coeff(3)/coeff(1)
      c=coeff(4)/coeff(1)
      q=(a*a-3.0d0*b)/9.0d0
      r=(2.0d0*a**3-9.0d0*a*b+2.7d1*c)/5.4d1
      r2=r*r
      q3=q**3
      disc=r2-q3
      if(disc.lt.0.0d0) then
        arg=r/dsqrt(q3)
        theta=dacos(arg)*third
        arg=-2.0d0*dsqrt(q)
        xreal(1)=arg*dcos(theta)-third*a
        xreal(2)=arg*dcos(theta+2.0d0*third*pi)-third*a
        xreal(3)=arg*dcos(theta-2.0d0*third*pi)-third*a
        ximage(1)=0.0d0
        ximage(2)=0.0d0
        ximage(3)=0.0d0
      else if (disc.ge.0.0d0) then
        if(r.lt.0.0d0) then
           arg=(dabs(r)+dsqrt(disc))**third
        else
           arg=-(dabs(r)+dsqrt(disc))**third
        endif
        if(arg.eq.0.0d0) then
           barg=0.0d0
        else
           barg=q/arg
        endif
        xreal(1)=arg+barg-third*a
        ximage(1)=0.0d0
        xreal(2)=-0.5d0*(arg+barg)-third*a
        ximage(2)=0.5d0*dsqrt(3.0d0)*(arg-barg)
        xreal(3)=xreal(2)
        ximage(3)=-ximage(2)
      endif
      return
      end   
      
      
