!*****************************************************
!*    Program to demonstrate the complex domain      *
!*            Mueller's subroutine                   *
!* ------------------------------------------------- *
!* Reference: BASIC Scientific Subroutines, Vol. II  *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
!*                                                   *
!*                F90 version by J-P Moreau, Paris.  *
!*                        (www.jpmoreau.fr)          *
!* ------------------------------------------------- *
!* Example:   Find a complex root of zē + 1 = 0      *
!*                                                   *
!* SAMPLE RUN:                                       *
!*                                                   *
!* Write the initial guesses and their bounds:       *
!*                                                   *
!*    X0 = 0                                         *
!*    Bond on X0 = 3                                 *
!*    Y0 = 0                                         *
!*    Bond on Y0 = 3                                 *
!*                                                   *
!* Convergence criterion: 0                          *
!*                                                   *
!* Maximum number of iterations: 100                 *
!*                                                   *
!* The calculated zero is (X,Y) =                    *
!* -0.000063  -1.000283                              *
!*                                                   *
!* The associated Z value is (U,V) =                 *
!* -0.000566   0.000126                              *
!*                                                   *
!* The number of steps was: 100                      *
!*****************************************************
PROGRAM DEMO_ZMUELLER

integer k,n
real*8  b1,b2,e,u,v,x,y,x0,y0

  print *,' '
  print *,'Write the initial guesses and their bounds:'
  print *,' '
  write(*,"('    X0 = ')",advance='no')
  read *, x0
  write(*,"('    Bond on X0 = ')",advance='no')
  read *, b1
  print *,' '
  write(*,"('    Y0 = ')",advance='no')
  read *, y0
  write(*,"('    Bond on Y0 = ')",advance='no')
  read *, b2
  print *,' '
  write(*,"(' Convergence criterion: ')",advance='no')
  read *, e
  write(*,"(' Maximum number of iterations: ')",advance='no')
  read *, n

  call ZMueller(x0,y0,b1,b2,e,n,x,y,k)      !Call complex domain Mueller subroutine

  write(*,50) x,y
  call Z(x,y,u,v)
  write(*,60) u,v
  write(*,70) k

  stop

50 format(//' The calculated zero is (X,Y) = ',F9.6,' ',F9.6) 
60 format(' The associated Z value is (U,V) = ',F9.6,' ',F9.6) 
70 format(' The number of steps was: ',i4//)

end

!**********************************************
!  Function subroutine 
! ---------------------------------------------               
Subroutine Z(x,y,u,v) 
  real*8 x,y,u,v
  u = x*x-y*y+1.d0
  v = 2.d0*x*y
  return
end
!**********************************************


!*************************************************
!* Mueller's method for complex roots subroutine *
!* --------------------------------------------- *
!* This routine uses the parabolic fitting tech- *
!* nique associated with Mueller"s method, but   *
!* does it in the complex domain.                *
!* --------------------------------------------- *
!* INPUTS:                                       *
!*   x0, y0 - The initial guess                  *
!*   b1, b2 - A bound on the error in this guess *
!*   e - The convergence criteria                *
!*   n - The maximum number of iterations        *
!* OUTPUTS:                                      *
!*   x,y - The value of the complex root found   *
!*   k - The number of iterations performed,     *
!*************************************************
Subroutine ZMueller(x0,y0,b1,b2,e,n,x,y,k)
  !Labels: 100,200,300
  integer k, n
  real*8  x0,y0,b1,b2,e,x,y
  real*8  d,d1,d2,x1,x2,x3,y1,y2,y3;
  real*8  u1,u2,u3,v1,v2,v3,xl1,xl2;
  real*8  a,a1,a2,b,c1,c2,e1,e2;
  ! Start calculations
  k=1
  x3=x0; y3=y0; x1=x3-b1; y1=y3-b2; x2=x3+b1; y2=y3+b2
100 d=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)
  ! Avoid divide by zero
  if (d.eq.d0) d=1.d-7;
  xl1=(x3-x2)*(x2-x1)+(y3-y2)*(y2-y1);
  xl1=xl1/d;
  xl2=(x2-x1)*(y3-y2)+(x3-x2)*(y2-y1);
  xl2=xl2/d;
  d1=(x3-x1)*(x2-x1)+(y3-y1)*(y2-y1);
  d1=d1/d;
  d2=(x2-x1)*(y3-y1)+(x3-x1)*(y2-y1);
  d2=d2/d;
  ! Get function values
  call Z(x1,y1,u1,v1);
  call Z(x2,y2,u2,v2);
  call Z(x3,y3,u3,v3);
  ! Calculate Mueller parameters
  e1=u1*(xl1*xl1-xl2*xl2)-2.d0*v1*xl1*xl2-u2*(d1*d1-d2*d2);
  e1=e1+2.d0*v2*d1*d2+u3*(xl1+d1)-v3*(xl2+d2);
  e2=2.d0*xl1*xl2*u1+v1*(xl1*xl1-xl2*xl2)-2.d0*d1*d2*u2-v2*(d1*d1-d2*d2);
  e2=e2+u3*(xl2+d2)+v3*(xl1+d1);
  c1=xl1*xl1*u1-xl1*xl2*v1-d1*xl1*u2+xl1*d2*v2+u3*xl1;
  c1=c1-u1*xl2*xl2-v1*xl1*xl2+u2*xl2*d2+v2*d1*xl2-v3*xl2;
  c2=xl1*xl2*u1+xl1*xl1*v1-d2*xl1*u2-xl1*d1*v2+v3*xl1;
  c2=c2+u1*xl1*xl2-v1*xl2*xl2-u2*xl2*d1+v2*d2*xl2+u3*xl2;
  b1=e1*e1-e2*e2-4.d0*(u3*d1*c1-u3*d2*c2-v3*d2*c1-v3*d1*c2);
  b2=2.d0*e1*e2-4.d0*(u3*d2*c1+u3*d1*c2+v3*d1*c1-v3*d2*c2);
  ! Guard against divide by zero
  if (b1.eq.d0) b1=1.d-7;
  a=datan(b2/b1); a=a/2.d0;
  b=dsqrt(dsqrt(b1*b1+b2*b2));
  b1=b*dcos(a); b2=b*dsin(a);
  a1=(e1+b1)*(e1+b1)+(e2+b2)*(e2+b2);
  a2=(e1-b1)*(e1-b1)+(e2-b2)*(e2-b2);
  if (a1>a2) goto 200
  a1=e1-b1; a2=e2-b2; goto 300
200 a1=e1+b1; a2=e2+b2;
300 a=a1*a1+a2*a2;
  xl1=a1*d1*u3-a1*d2*v3+a2*u3*d2+a2*v3*d1;
  ! Guard against divide by zero
  if (a.eq.d0) a=1.d-7;
  xl1=-2.d0*xl1/a;
  xl2=-d1*u3*a2+d2*v3*a2+a1*u3*d2+a1*v3*d1;
  xl2=-2.d0*xl2/a;
  ! Calculate new estimate
  x=x3+xl1*(x3-x2)-xl2*(y3-y2);
  y=y3+xl2*(x3-x2)+xl1*(y3-y2);
  ! Test for convergence
  if (dabs(x-x0)+dabs(y-y0)<e) return;
  ! Test for number of iterations
  if (k>=n) return;
  ! Continue
  k=k+1; x0=x; y0=y; x1=x2; y1=y2; x2=x3; y2=y3; x3=x; y3=y;
  goto 100
end ! ZMueller()

! End of file zmueller.f90