*-------------------------
	function asp( xs, ys, npl, xpl, ypl, qpl )
* evaluate magnification for lensing by star + npl planets
* Input:
*	xs,ys	r4 source position
*	npl	i4 number of planets
*	xpl(npl) r4 planet x coordinate
*	ypl(npl) r4 planet y coordinate
*	qpl(npl) r4 planet/star mass ratio
* Output:
*	asp	r4 total amplification
* Output to common:
*	nimage	i4 number of images
*	x_i(i)	r4 x coordinate of image i
*	y_i(i)	r4 y coordinate of image i
*	a_i(i)	r4 amplification of image i
* 1998 Sep Keith Horne @ St-Andrews
* 2002 Mar KDH @ St.And -- lens2
* 2002 Apr KDH @ St.And -- lens1
* 2003 Jan KDH @ St.And -- thresh 1e-5->1e-6
* 2003 Jan KDH @ St.And -- image-planet distance in theta and ln(u)
* 2010 Apr KDH @ St.And -- simplify code
	real*4 xpl(*), ypl(*), qpl(*)
	character*1 method
	common/nimage/	nimage
	common/x_i/	x_i(2)
	common/y_i/	y_i(2)
	common/a_i/	a_i(2)
	data kold/-1/
	pi = 4. * atan2( 1., 1. )
	twopi = 2. * pi

* threshold for including a planet
c	thresh = 1.e-5
c	thresh = 1.e-6
	thresh = 1.e-7

c	method = 'G'
	method = 'B'

* gould & loeb model
c      if( method .eq. 'G' ) then
c	call gould_loeb_92( xs, ys, npl, xpl, ypl, qpl, val )
* binary lens model
c      else

* no planet model
	call lens1( xs, ys, val )
	a1 = val
	asp = val

* trap error
	alo = 0.999
      if( a1 .lt. alo ) then
	write(*,*) '** ERROR in asp xs=', xs, ' ys=', ys, ' npl', npl
	write(*,*) '** PSPL magnification', a1,' <', alo, ' :((((('
      do i=1,npl
	write(*,*) '** planet', i, ' q=', qpl(i)
     &		, ' xp=', xpl(i), ' yp=', xpl(i)
      end do
	stop 
      end if

	if( npl .eq. 0 ) return
	call mnmx( npl, qpl, qmn, qmx )
	if( qmx .le. 0. ) return

* one or more planets
* find 1 or 2 most important planets
* SHOULD DO ALL PLANETS IN ORDER OF DECREASING IMPORTANCE
	k1 = 0
	k2 = 0
	delta = 0.
	best = 0.
* planet loop
      do k = 1, npl
	q = qpl( k )
      if( q .gt. 0. ) then
	xp = xpl( k )
	yp = ypl( k )
	uup = xp * xp + yp * yp
	azp = atan2( yp, xp )
	test = 0.
	delta = 0.
* image loop
      do i = 1, 2
* x-y separation of planet and image
      if( .false. ) then
	xip = x_i( i ) - xp
	yip = y_i( i ) - yp
	uu = ( xip * xip + yip * yip ) / q
* image-planet separation in azimuth and ln(u)
      else
	ai = a_i( i )
	xi = x_i( i )
	yi = y_i( i )
	uui = xi * xi + yi * yi
	dlnu = 0.5 * alog( uui / uup )
	azi = atan2( yi, xi )
	daz = mod( azi - azp, twopi )
	if( daz .gt. pi ) daz = daz - twopi
	if( daz .lt. -pi ) daz = daz + twopi
	uu = ( dlnu * dlnu + daz * daz ) * uui / q
      end if
	uip = sqrt( uu )
	delta =  ( atotl( uip ) - 1. ) * ai / a1
	test = max( test, abs( delta ) )
* next image
      end do
* stow best
      if( k1 .eq. 0 ) then
	best = test
	k1 = k
      else if( test .gt. best ) then
	second = best
	best = test
	k2 = k1
	k1 = k
      else if( test .gt. second ) then
	second = test
	k2 = k
      end if

* next planet
      end if
      end do

* first planet
      if( k1 .gt. 0 ) then
      if( abs( best ) .ge. thresh ) then
	call lens2( xs, ys, xpl(k1), ypl(k1), qpl(k1), val )
	asp = val
* second planet
      if( k2 .gt. 0 ) then
      if( abs( second ) .gt. thresh ) then
	call lens2( xs, ys, xpl(k2), ypl(k2), qpl(k2), val )
	asp = asp + ( val - a1 )
      end if
      end if
      end if
      end if

* report errors and changes in lens model
      if( k1 .ne. kold .or. a1 .lt. alo .or. asp .lt. alo ) then
	u = sqrt( xs*xs + ys*ys )
	write(*,*) 'PSPL model xs,ys', xs, ys,' A(', u, ')=', a1
	delta = asp / a1 - 1.
	write(*,*) 'Planets:', npl, ' del', delta, ' A', asp
	write(*,*) 'Planet', k1, ' del', best,
     &		' q', qpl(k1), ' x,y', xpl(k1), ypl(k1)
	if( k2 .gt. 0 )
     &	write(*,*) 'Planet', k2, ' del', second,
     &		' q', qpl(k2), ' x,y', xpl(k2), ypl(k2)
	kold = k1
      end if

c      end if


	return
	end

