	subroutine lens2( xs, ys, xp, yp, qp, amp )
* gravitational lens - 2 point masses
* Input:
*	xs, ys	r4 source coordinates (origin=M1, unit=R_E(M1+M2) )
*	xp, yp	r4 planet coordinates (origin=M1, unit=R_E(M1+M2) )
*	qp	r4 planet/star mass ratio (M2/M1)
* Output:
*	amp	r4 amplification
* Output to common blocks:
*	nimage	i4 number of images
*	x_i()	r4 image x coordinate
*	y_i()	r4 image y coordinate
*	a_i()	r4 image amplification
* 2002 Mar Keith Horne @ St.Andrews
* 2008 Jul KDH @ St.A - fix error when all images rejected.
	parameter ( maxim = 5 )
	common/nimage/	nimage
	common/x_i/	x_i( maxim )
	common/y_i/	y_i( maxim )
	common/a_i/	a_i( maxim )
	common/err_i/	err_i( maxim )
* rotate planet from (xp,yp) to (0,d)
	d = sqrt( xp*xp + yp*yp )
	cth = xp / d
	sth = yp / d
	xsp = cth * xs + sth * ys
	ysp = - sth * xs + cth * ys
* find images and sum amplifications
	call lens2a( qp, d, xsp, ysp, amp )
* de-rotate images
      do i=1,nimage
	xi = x_i(i)
	yi = y_i(i)
	x_i(i) = cth * xi - sth * yi
	y_i(i) = sth * xi + cth * yi
      end do
* report
c	write(*,*) 'images:', nimage, ' amp:', amp
c      do i=1,nimage
c	write(*,*) 'xya:', i, x_i(i), y_i(i), a_i(i) 
c      end do

* report failures
	amplo = 0.99
      if( nimage .le. 0 .or. amp .lt. amplo ) then
	write(*,*) '** ERROR lens2 A=', amp, ' <', amplo, ' :((((('
	write(*,*) '** source xy', xs, ys, ' planet xy', xp, yp, ' q', qp
	write(*,*) '** lens2a(', qp, d, xsp, ysp, ' )=', amp
	write(*,*) '** images:', nimage
      do i=1,nimage
	x = x_i(i)
	y = y_i(i)
	a = a_i(i)
	err = err_i(i)
	write(*,*) '** image', i, 'x,y', x, y, ' err', err, ' A', a
      end do
      end if

	return
	end
*----------------------------------------------------
      subroutine lens2d( q, sep, xs, ys, delta )
* Calculates amplification excess for the binary lens
* Input:
*	q	r4 planet/star mass ratio m2/m1
*	sep	r4 binary separation: 1 = R_E(m1+m2)
*	xs,ys	r4 source coords: (0,0) = m1, (sep,0)=m2
* Output:
*	delta	r4 amplification(m1+m2) / amplification(m1)
* 2002 Apr Keith Horne @ St-Andrews
	call lens2a( q, sep, xs, ys, amp )
	us = sqrt( xs*xs + ys*ys )
	delta = amp / atotl( us )
	return
	end
*----------------------------------------------------
	subroutine lens2a( qin, sin, xin, yin, amp )
* image positions and amplifications for the binary lens
* Input:
*	qin	r4 planet/star mass ratio m2/m1
*	sin	r4 binary separation: 1 = R_E(m1+m2)
*	xin,yin	r4 source coords: (0,0) = m1, (sep,0)=m2
* Output:
*	amp	r4 amplification
* Output to common:
*	nimage	i4 number of images
*	x_i(*)	r4 image x coords
*	y_i(*)	r4 image y coords
*	a_i(*)	r4 image amplifications
*	err_i(*) r4 error in predicted image position
* uses ZROOTS and LAGUER from Numerical Recipes.
* 1998 written by R.Carson
* 2001 modified by Y.Tsapras and converted to IDL
* 2002 Mar Keith Horne @ St-Andrews - simplify for speed
* 2002 Apr KDH @ St-And - simplify inputs
* 2008 Jul KDH @ St-And - keep best image and others almost as close
* 2008 Jul KDH @ SAAO - i_ir(i) to rank and keep(i) to keep images
* 2008 Jul KDH @ SAAO - amp = max( 1., amp )
* 2008 Aug KDH @ StA - try implicit real*8
* 2008 Aug KDH @ StA - use dzroots
* 2008 Aug KDH @ StA - keep 3 or 5 roots
	implicit double complex ( z )
	implicit real*8 (a-h,o-y)
	real*4 qin, sin, xin, yin, amp
	real*4 x_i, y_i, a_i, err_i
	parameter ( maxim = 5 )
* image positions
	common/nimage/	nimage
	common/x_i/	x_i( maxim )
	common/y_i/	y_i( maxim )
* image magnifications (can be negative)
	common/a_i/	a_i( maxim )
* error in predicted source position
	common/err_i/	err_i( maxim )
* ranked image numbers
	integer*4 i_ir( maxim )
* keep images
	logical keep( maxim )
* Complex coefficients of 5th order polynomial
	dimension zc( maxim + 1 )
* Complex image positions
	dimension zi( maxim )
	real*8 sum, time1, time2, gm1, gm2

* double the input arguments
	q = dble( qin )
	xs = dble( xin )
	ys = dble( yin )
	sep = dble( sin )

* Criterion for checking the validity of solutions:
* Map images back to source and keep if error less than ERRMX
	ERRMX = 1.d-4

* masses (m1+m2=1)
	gm1 = 1.d0 / ( 1.d0 + q )
	gm2 = q * gm1
* half binary separation
	d = sep * 0.5d0
* lens configuration
	Z1 = dCMPLX( -d, 0.d0 )
	Z2 = dCMPLX( +d, 0.d0 )
* source position
	ZS = dCMPLX( xs - d, ys )
* shift origin to m1
c	z0 = z1
c	z1 = z1 - z0
c	z2 = z2 - z0
c	zs = zs - z0
* complex conjugates
	Z1C = Z1
	Z2C = Z2
	ZSC = dCMPLX( dreal( zs ), - dimag( zs ) )

* complex coefficients of 5th order polynomial
	Z1_Z1 = Z1 * Z1
	ZSC_ZSC = ZSC * ZSC
	ZC( 6 ) = Z1_Z1 - ZSC_ZSC

	SM = GM1 + GM2
	DM = GM2 - GM1
	ZA = DM + Z1 * ZS
	Z1_ZA = Z1 * ZA
	ZF =  ZS * ZSC_ZSC - Z1_ZA  
	ZSC_SM =  ZSC * SM
	ZC( 5 ) = ZF - ZSC_SM

	SM_2 = SM + SM
	ZS_ZSC = ZS * ZSC
	Z1_Z6 = Z1 * ZC( 6 )
	ZB = ( DM + DM ) * ZSC
	ZC( 4 ) = SM_2 * ZS_ZSC + Z1 * ( ZB - Z1_Z6 - Z1_Z6 )

	SDM = SM * DM
	DDM = DM * DM
	SSM = SM * SM
	ZS_ZB = ZS * ZB
	Z1_ZF = Z1 * ZF
	ZC( 3 ) = SSM * ZS + Z1 * ( SDM - ZS_ZB - Z1_ZF - Z1_ZF )

	ZC( 2 ) = -Z1 * ( (SDM + SDM ) * ZS
     &	+ Z1 * ( DDM + SSM + SM_2 * ZS_ZSC
     &	+ Z1 * ( ZB - Z1_Z6 ) ) )

	ZC( 1 ) = Z1_Z1 * ( DDM * ZS + Z1 * ( SDM + ZS_ZB
     &	+ Z1 * ( ZF + ZSC_SM ) ) )

* order of polynomial
	M = 5
      if( zc( 6 ) .ne. 0. ) then
	m = 5
      else if( zc( 5 ) .ne. 0. ) then
	m = 4
      else if( zc( 4 ) .ne. 0. ) then
	m = 3
      else if( zc( 3 ) .ne. 0. ) then
	m = 2
      else if( zc( 2 ) .ne. 0. ) then
	m = 1
      else
	m = 4
      end if

* report coeficients
      if( m .ne. 5 ) then
	m = 5
	write(*,*) 'Poly Coefficients:', m+1, ' q', qin
	write(*,*) 'zs:', zs, ' zs_zsc:', zs_zsc
	write(*,*) 'z1:', z1, ' z2:', z2
	write(*,*) 'za:', za, ' zb:', zb
	write(*,*) 'zf:', zf, ' zsc_sm:', zsc_sm
      do i=1,m+1
	write(*,*) 'zc(', i, '):', zc(i)
      end do
      end if
* find roots
c	CALL ZROOTS( ZC, M, ZI, .TRUE. )
	CALL dZROOTS( ZC, M, ZI, .TRUE. )

* shoot back to source positions and compute magnifications 
      DO I = 1, M
	XI = dREAL( ZI( I ) )
	YI = dIMAG( ZI( I ) )
	ZIC = dCMPLX( XI, -YI )
	ZDC1 = ZIC - Z1C
	ZDC2 = ZIC - Z2C
	ZB = GM1 / ZDC1 + GM2 / ZDC2
* ZE map back position of source
	ZE = ZI( I ) - ZB
* reject invalid image
	ERR = CdABS( ZE - ZS )
* magnification
	ZD = GM1 / ( ZDC1 * ZDC1 ) + GM2 / ( ZDC2 * ZDC2 )
	CJ = CdABS( ZD )
	CJ = CJ * CJ
	VJ = 1.D0 - CJ
	AJ = 1.D0 / VJ
* stow
	x_i( i ) = xi + d
	y_i( i ) = yi
	a_i( i ) = aj
	err_i( i ) = err
* next image
      END DO

* sort the images by error in predicted source position
	call quicksort( m, err_i, i_ir )

* decide which images to keep
* Keep 3 best images, or 2 more if they are both close
	sum = 0.d0
	nim = 3
	i = i_ir( 5 )
	if( abs( err_i( i ) ) .le. errmx ) nim = 5
      do ir = 1, m
	i = i_ir( ir )
	err = abs( err_i( i ) )
c	keep(i) = err .le. errmx .or. ir .le. nim
	keep(i) = ir .le. nim
      if( keep(i) ) then
	sum = sum + abs( a_i(i) )
      end if
      end do
* summed magnification
	amp = sum

* report error if magnification < 1 (and try to recover)
	amplo = 0.99
      if( amp .lt. amplo ) then
	write(*,*) '** ERROR in lens2a. A=', amp, ' <', amplo, ' :(((.'
	us = sqrt( xs * xs + ys * ys )
	write(*,*) 'Planet u=', sep, ' Source xy=', xs, ys, ' u=', us
* Keep the image closest to A=1.
	k = 0
	best = 1000.
      do i = 1, m
	test = abs( a_i(i) - 1. )
      if( test .lt. best ) then
	k = i
	best = test
      end if
      end do
	if( k .gt. 0 ) keep(k) = .true.
* report images
      if( .false. ) then
      do ir = 1, m
	i = i_ir(ir)
	a = a_i(i)
	err = err_i(i)
	write(*,*) '** Image', i, ' rank', ir, ' err', err
     &		, ' A', a, ' keep', keep(i)
      end do
      end if
      end if

* trim the image list
	n = 0
      do i = 1, m
      if( keep(i) ) then
	n = n + 1
	x_i( n ) = x_i(i)
	y_i( n ) = y_i(i)
	a_i( n ) = abs( a_i(i) )
	err_i( n ) = err_i(i)
      end if
      end do
	nimage = n

* report error in number of images
      if( nimage .ne. 3 .and. nimage .ne. 5 ) then
	write(*,*) '** WARNING: nimage', nimage, ' not 3 or 5 :((((.'
	write(*,*) '** q', qin, ' sep', sin, ' xs,ys', xin, yin
	sum = 0.d0
      do ir=1,nimage	
	i = i_ir(ir)
	a = a_i(i)
	sum = sum + dabs( a )
	err = err_i(i)
	write(*,*) '**', i, ' rank', ir, ' keep', keep(i) 
     &		, ' err', err, ' Ai', a, ' A', sum
      end do
      end if

* lensing should never yield total magnification below 1.
	amp = max( amp, 1. )
* fini
	RETURN
	END
