!******************************************************
!*   ROOTS OF A SECOND DEGREE EQUATION WITH COMPLEX   *
!*   COEFFICIENTS  AZ2 + BZ + C = 0                   *
!* -------------------------------------------------- *
!* SAMPLE RUN:                                        *
!* (Find complex roots of equation:                   *
!* (-4+i)z2 + (2-3i)z +(5-i) = 0).                    *
!*                                                    *
!* SOLVING A COMPLEX 2ND DEGREE EQUATION              *
!*     az2 + bz + c = 0                               *
!*                                                    *
!*  a real part      = -4                             *
!*  a imaginary part = 1                              *
!*  b real part      = 2                              *
!*  b imaginary part = -3                             *
!*  c real part      = 5                              *
!*  c imaginary part = -1                             *
!*                                                    *
!*  Root 1 =   1.444644 -  0.352759 i                 *
!*  Root 2 =  -0.797586 -  0.235476 i                 *
!*                                                    *
!* -------------------------------------------------- *
!* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp  *
!* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"   *
!* [BIBLI 05].                                        *
!*                                                    *
!*                  F90 Release By J-P Moreau, Paris. *
!*                          (www.jpmoreau.fr)         *
!******************************************************
Program TEQUA2C

  COMPLEX  a,b,c,s1,s2

  print *,' '
  print *,' SOLVING A COMPLEX 2ND DEGREE EQUATION'
  print *,'      az2 + bz + c := 0'
  print *,' '
  write(*,10,advance='no'); read *, ax
  write(*,11,advance='no'); read *, ay
  a = CMPLX(ax,ay)
  write(*,20,advance='no'); read *, bx
  write(*,21,advance='no'); read *, by
  b = CMPLX(bx,by)
  write(*,30,advance='no'); read *, cx
  write(*,31,advance='no'); read *, cy
  c = CMPLX(cx,cy)

  call equa2c(a,b,c,s1,s2)

  print *,' '
  write(*,40,advance='no') REAL(s1)
  if (IMAG(s1) > 0.) then 
    write(*,50) abs(IMAG(s1))
  else 
    write(*,51) abs(IMAG(s1))
  end if
  write(*,41,advance='no') REAL(s2)
  if (IMAG(s2) > 0.) then 
    write(*,50) abs(IMAG(s2))
  else 
    write(*,51) abs(IMAG(s2))
  end if
  print *,' '
  stop
10 format('    a      real part = ')
11 format('    a imaginary part = ')
20 format('    b      real part = ')
21 format('    b imaginary part = ')
30 format('    c      real part = ')
31 format('    c imaginary part = ')
40 format(' Root 1 = ',F10.6)
41 format(' Root 2 = ',F10.6)
50 format(' + ',F10.6,' i')
51 format(' - ',F10.6,' i')
END

Subroutine equa2c(a,b,c,s1,s2)
COMPLEX a,b,c,s1,s2
REAL*8 r,sx,sy,ux,uy,vx,vy,wx,wy
  ux = REAL(b)*REAL(b) - IMAG(b)*IMAG(b) - 4.d0 * REAL(a)*REAL(c) + 4.d0 * IMAG(a)*IMAG(c)
  uy = 2.d0*REAL(b)*IMAG(b) - 4.d0 * REAL(a)*IMAG(c) - 4.d0*IMAG(a)*REAL(c)
  r = dsqrt(ux*ux + uy*uy)
  vx = dsqrt((r+ux)/2.d0)
  vy = dsqrt((r-ux)/2.d0)
  if (uy<0.)  vy = -vy
  wx = (-REAL(b) - vx)/2.d0
  wy = (-IMAG(b) - vy)/2.d0
  ux = (-REAL(b) + vx)/2.d0
  uy = (-IMAG(b) + vy)/2.d0
  r = REAL(a)*REAL(a) + IMAG(a)*IMAG(a)
  sx = (REAL(a)*wx + IMAG(a)*wy)/r
  sy = (REAL(a)*wy - IMAG(a)*wx)/r
  s1 = CMPLX(sx,sy)
  sx = (REAL(a)*ux + IMAG(a)*uy)/r
  sy = (REAL(a)*uy - IMAG(a)*ux)/r
  s2 = CMPLX(sx,sy)
  return
End

!end of file tequa2.f90