*-------------------------
        subroutine lens1( xs, ys, amp )
* gravitational lens - 1 point mass
* Input:
*       xs,ys   r4 source coordinates (origin=M1, unit=R_E)
* Output:
*       amp     r4 total magnification
* 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 magnification of image i
* 2002 Apr Keith Horne @ St-Andrews
        common/nimage/  nimage
	common/x_i/     x_i( 2 )
	common/y_i/     y_i( 2 )
	common/a_i/     a_i( 2 )
* unlensed star position
	us = sqrt( xs * xs + ys * ys )
c	write(*,*) 'star at x,y,u:', xs, ys, us
	cth = xs / us
	sth = ys / us
	nimage = 2
	a_i( 1 ) = aplus( us )
	u = uplus( us )
	x_i( 1 ) = u * cth
	y_i( 1 ) = u * sth
	a_i( 2 ) = aminus( us )
	u = uminus( us )
	x_i( 2 ) = u * cth
	y_i( 2 ) = u * sth
	amp = a_i( 1 ) + a_i( 2 )
	return
	end
*--------------------
	function uplus( u )
* Input:	u	r4 separation / R_E of source
* Output:	uplus	r4 separation / R_E of major image
* 1998 Aug Keith Horne @ St.Andrews
	half = 0.5 * u
	uplus = half + sqrt( half*half + 1. )
	return
	end
*--------------------
	function uminus( u )
* Input:	u	r4 separation / R_E of source
* Output:	uminus	r4 separation / R_E of minor image
* 1998 Aug Keith Horne @ St.Andrews
	a = 0.5 * u
	root = sqrt( a * a + 1. )
	uminus = a - root
	return
	end
*--------------------
	function duplus( u )
* Input:	u	r4 separation / R_E of source
* Output:	duplus	r4 radial magnification of major image
* 1998 Aug Keith Horne @ St.Andrews
	a = 0.5 * u
	root =  sqrt( a * a + 1. )
	duplus = 0.5 * ( 1. + a / root )
	return
	end
*--------------------
	function duminus( u )
* Input:	u	r4 separation / R_E of source
* Output:	duminus	r4 radial magnification of minor image
* 1998 Aug Keith Horne @ St.Andrews
	a = 0.5 * u
	root =  sqrt( a * a + 1. )
	duminus = 0.5 * ( 1. - a / root )
	return
	end
*--------------------
	function aplus( u )
* Input:	u	r4 separation / R_E of source
* Output:	aplus	r4 magnification of major image
* 1998 Aug Keith Horne @ St.Andrews
	aplus = abs( duplus( u ) * uplus( u ) / u )
	return
	end
*--------------------
	function aminus( u )
* Input:	u	r4 separation / R_E of source
* Output:	aminus	r4 magnification of minor image
* 1998 Aug Keith Horne @ St.Andrews
	aminus = abs( duminus( u ) * uminus( u ) / u ) 
	return
	end
*--------------------
	function atotl( u )
* Input:	u	r4 separation / R_E of source
* Output:	atotl	r4 total magnification
* 1998 Aug Keith Horne @ St.Andrews
	usq = u * u
      if( u .gt. 1. ) then
	atotl = ( usq + 2. ) / ( abs(u) * sqrt( usq + 4. ) )
      else
	bot = usq * ( usq + 4. )
	atotl = sqrt( 1. + 4. / bot )
      end if
	return
	end
*--------------------
	function ainv( a )
* Input:	a	r4 magnification
* Output:	ainv	r4 separation / R_E of source
* 1998 Aug Keith Horne @ St.Andrews
	ainv = -1.
	if( a .lt. 1. ) return
	root = sqrt( 1. - 1. / ( a * a ) )
	ainv = sqrt( 2. / root - 2. )

	return
	end
*--------------------
	function thalf( a0 )
* time from event centre to half-peak magnification
* in units of t_e
* Input:	a0	r4 peak magnification
* Output:	thalf_a	r4 time to half peak / t_E
* 2008 Jun Keith Horne @ St.Andrews
	thalf = -1.
	if( a0 .lt. 1. ) return
	u = ainv( ( 1. + a0 ) * 0.5 )
	u0 = ainv( a0 )
	thalf = sqrt( u*u - u0*u0 )
	return
	end
*--------------------
	function thalf_u0( u0 )
* time from event centre to half-peak magnification
* in units of t_e
* Input:	u0	r4 impact parameter in R_E units
* Output:	thalf_a	r4 time to half peak / t_E
* 2009 Jan Keith Horne @ St.Andrews
	thalf_u0 = 0.
	if( u0 .eq. 0. ) return
	a0 = atotl( u0 )
	u = ainv( ( 1. + a0 ) * 0.5 )
	thalf_u0 = sqrt( u*u - u0*u0 )
	return
	end

*---------------------
	function pspl( u )
* PSPL magnification
* 2008 Jul Keith Horne @ SAAO
	real*8 uu, au
	uu = u * u
	uu = max( uu, 1.d-40 )
	au = dsqrt( uu )
	pspl = ( uu + 2.d0 ) / ( au * dsqrt( uu + 4.d0 ) )
	return
	end
*---------------------
	function dpspl( u )
* 1st derivitive of PSPL magnification
* 2008 Jul Keith Horne @ SAAO
	real*8 uu, uu4
	uu = u * u
	uu = max( uu, 1.d-40 )
	uu4 = uu + 4.d0
	dpspl = -8.d0 / ( uu * uu4 * dsqrt( uu4 ) )
	return
	end
*---------------------
	function ddpspl( u )
* 2nd derivitive of PSPL magnification
* 2008 Jul Keith Horne @ SAAO
	real*8 uu, uu4, bb
	uu = u * u
	uu = max( uu, 1.d-40 )
	uu4 = uu + 4.d0
	bb = uu * uu4
	ddpspl = 8.d0 * ( 8.d0 - 5.d0 * uu ) / ( bb * dsqrt( bb ) )
	return
	end

