c..this file contains root finding routines:   
c..  
c..quadratic, cubic and quartic roots: 
c..routine sroot finds the complex roots of a quadratic with real coefficients 
c..routine csroot does the above with complex coefficients 
c..routine croot finds the complex roots of a cubic with real coefficients 
c..routine ccroot does the above with complex coefficients 
c..routine q4oot finds the complex roots of a quartic with real coefficients 
c..routine cq4root does the above with complex coefficients 



      subroutine sroot(a,b,c,r1,r2)  
      implicit none
      save

c..this routine finds the roots of the quadratic equation 
c..           a*x**2 + b*x + c = 0  
c..input are the real coefficients a,b and c.    
c..output are the complex roots r1 and r2.   

c..declare the pass
      double precision  a,b,c
      double complex    r1,r2 

c..local variables
      double precision  tiny
      parameter         (tiny = 1.0d-14)  
      double complex    aa,bb,cc

      if (abs(a) .le. tiny) stop 'coef a is zero in routine sroot'  
      aa  = dcmplx(a,0.0d0) 
      bb  = dcmplx(b,0.0d0) 
      cc  = dcmplx(c,0.0d0) 
      call csroot(aa,bb,cc,r1,r2) 

      return 
      end    



      subroutine csroot(a,b,c,r1,r2)  
      implicit none
      save

c..this routine finds the roots of the quadratic equation 
c..           a*x**2 + b*x + c = 0  
c..input are the complex coefficients a,b and c.    
c..output are the complex roots r1 and r2.   

c..declare the pass
      double complex    a,b,c,bb,cc,r1,r2

c..local variables
      double precision  asgn,tiny 
      parameter         (tiny  = 1.0d-14)  
      double complex    fac 

      if (abs(a) .le. tiny) stop 'a is zero in routine csroot'  
      bb   = conjg(b) 
      cc   = sqrt(b*b - 4.0d0*a*c) 
      asgn = sign(1.0d0,dreal(bb * cc)) 
      fac  = -0.5d0*(b + asgn*cc) 
      r1   = fac/a 
      r2   = c/fac 

      return 
      end    





      subroutine croot(a,b,c,d,r1,r2,r3)  
      implicit none
      save

c..this routine finds the roots of the cubic equation  
c..          a*x**3 + b*x**2 + c*x + d = 0  
c..input are the real coefficients a,b,c, and d. 
c..output are the complex roots r1, r2, and r3. 

c..declare the pass
      double precision  a,b,c,d
      double complex    r1,r2,r3 

c..local variables
      double precision  tiny
      parameter         (tiny = 1.0d-14) 
      double complex    aa,bb,cc,dd

      if (abs(a) .le. tiny) stop 'coef a is zero in routine croot'  
      aa  = dcmplx(a,0.0d0) 
      bb  = dcmplx(b,0.0d0) 
      cc  = dcmplx(c,0.0d0) 
      dd  = dcmplx(d,0.0d0) 
      call ccroot(aa,bb,cc,dd,r1,r2,r3) 

      return 
      end 





      subroutine ccroot(a,b,c,d,r1,r2,r3)   
      implicit none
      save

c..this routine finds the roots of the cubic equation 
c..          a*x**3 + b*x**2 + c*x + d = 0  
c..input are the complex coefficients a,b,c, and d. 
c..output are the complex roots r1, r2, and r3. 
c.. 

c..declare the pass
      double complex     a,b,c,d,r1,r2,r3

c..local variables
c..the constant fac1 below is sqrt(3)/2 

      double precision   third,ninth,five4,asgn,fac1,tiny 
      parameter          (third = 1.0d0/3.0d0, 
     1                    ninth = 1.0d0/9.0d0, 
     2                    five4 = 1.0d0/54.0d0, 
     3                    fac1  = 0.8660254037844386d0, 
     4                    tiny  = 1.0d-14)  
      double complex     p,q,r,bigq,bigr,cc,rr,biga,bigb,fac


c..put it in standard form 
      if (abs(a) .le. tiny) stop 'a is zero in routine ccroot'  
      p    = b/a 
      q    = c/a 
      r    = d/a 

      bigq = ninth * (p*p - 3.0d0*q)    
      bigr = five4 * ( p*(2.0d0*p*p - 9.0d0*q) + 27.0d0*r) 
      cc   = sqrt(bigr*bigr - bigq*bigq*bigq) 
      rr   = conjg(bigr) 
      asgn = sign(1.0d0,dreal(rr*cc)) 

      biga = bigr + asgn*cc 
      if (abs(biga) .le. tiny) then 
       biga = dcmplx(0.0d0,0.0d0) 
       bigb = dcmplx(0.0d0,0.0d0) 
      else 
       biga = -1.0d0 * biga**third 
       bigb = bigq/biga 
      end if 

      fac = dcmplx(0.0d0,fac1)  
      cc   = -0.5d0 * (biga + bigb) - third*p  
      rr   = biga - bigb 

      r1   = biga + bigb - third*p 
      r2   = cc + fac * rr 
      r3   = cc - fac * rr 
      return 
      end    




      subroutine q4root(a,b,c,d,e,r1,r2,r3,r4)  
      implicit none

c..this routine finds the roots of the quartic equation  
c..       a*x**4 + b*x**3 + c*x**2 + d*x + e = 0  
c..input are the real coefficients a, b, c, d, and e 
c..output are the complex roots r1, r2, r3, and r4. 

c..declare the pass
      double precision  a,b,c,d,e
      double complex    r1,r2,r3,r4 

c..local variables
      double precision  tiny 
      parameter         (tiny = 1.0d-14)
      double complex    aa,bb,cc,dd,ee

      if (abs(a) .le. tiny) stop 'coef a is zero in routine q4root'  
      aa  = dcmplx(a,0.0d0) 
      bb  = dcmplx(b,0.0d0) 
      cc  = dcmplx(c,0.0d0) 
      dd  = dcmplx(d,0.0d0) 
      ee  = dcmplx(e,0.0d0) 
      call cq4root(aa,bb,cc,dd,ee,r1,r2,r3,r4) 

      return 
      end 




      subroutine cq4root(a,b,c,d,e,r1,r2,r3,r4)   
      implicit none
      save

c..this routine finds the roots of the quartic equation  
c..       a*x**4 + b*x**3 + c*x**2 + d*x + e = 0  
c..input are the complex coefficients a, b, c, d, and e. 
c..output are the complex roots r1, r2, r3, and r4. 

c..declare the pass
      double complex     a,b,c,d,e,r1,r2,r3,r4

c..local variables
      double precision   x1,x2,x3,x4,tiny 
      parameter          (tiny = 1.0d-14) 
      double complex     a3,a2,a1,a0,t1,t2,t3,t4,biga,bigb,bigc,bigd 

c..put it in standard form 
      if (abs(a) .le. tiny) stop 'a is zero in routine cq4root'  
      a3   = b/a 
      a2   = c/a 
      a1   = d/a 
      a0   = e/a 
      biga = 0.5d0*a3 

c..solve the cubic 
      t1   = (1.0d0,0.0d0) 
      t2   = -a2 
      t3   = a1*a3 - 4.0d0*a0 
      t4   = a0*(4.0d0*a2 - a3*a3) - a1*a1 
      call ccroot(t1,t2,t3,t4,r1,r2,r3) 

c..find the largest real root 
      x1 = dreal(r1) 
      x2 = dreal(r2) 
      x3 = dreal(r3) 
      x4 = max(x1,x2,x3) 
      if (x4 .eq. x1) then 
       bigb = 0.5d0*r1 
      else if (x4 .eq. x2) then 
       bigb = 0.5d0*r2 
      else if (x4 .eq. x3) then 
       bigb = 0.5d0*r3 
      end if 

c..get the other stuff 
      t1 = bigb*bigb - a0 
      if (abs(t1) .le. tiny) then 
        bigd = (0.0d0,0.0d0) 
        bigc = sqrt(biga*biga + 2.0d0*bigb - a2) 
      else 
        bigd = sqrt(t1) 
        bigc = (biga*bigb - 0.5d0*a1)/bigd 
      end if 

c..solve the two quadratics 
      t1 = (1.0d0,0.0d0) 
      t2 = biga - bigc 
      t3 = bigb - bigd 
      call csroot(t1,t2,t3,r1,r2) 

      t1 = (1.0d0,0.0d0) 
      t2 = biga + bigc 
      t3 = bigb + bigd 
      call csroot(t1,t2,t3,r3,r4) 

      return 
      end