	function dzwstd( q, u, smag0, b0 )
* detection zone area d(theta)d(lnu) for standard telescope 
* Input:
*	q	r4 mass ratio
*	u	r4 source-lens distance in Einstein ring radii
*	smag0	r4 baseline magnitude
*	b0	r4 blend flux / source flux
* Output:
*	dzwstd	r4 detection zone area for standard telescope
* 2008 Keith Horne @ SAAO
	data qold, uold, sold, bold/ -1., -1., -100., -1000. /
	cm2 = 1.e4
	bw = 1.e3
	t = 100.
	fwhm = 1.
	skymag = 19.
* use old value
      if( q .eq. qold .and. u .eq. uold
     &	.and. smag0 .eq. sold .and. b0 .eq. bold ) then
	dzwstd = dzwold
* recompute
      else
	dzwstd = dzwtel( q, u, smag0, b0, cm2, bw, t, fwhm, skymag )
      end if
* stow old values
	qold = q
	uold = u
	sold = smag0
	bold = b0
	dzwold = dzwstd
	return
	end
*--------------------------------------------------------------------
	function dzwtel( q, u, smag0, b0, cm2, bw, t, fwhm, skymag )
* detection zone area d(theta)d(lnu) for specified telescope 
* Input:
*	q	r4 mass ratio
*	u	r4 source-lens distance in Einstein ring radii
*	smag0	r4 baseline magnitude
*	b0	r4 blend flux / source flux
*	cm2	r4 effective area of telescope (cm2)
*	bw	r4 bandwidth (Angstroms)
*	fwhm	r4 seeing fwhm (arcsec)
*	skymag	r4 sky magnitude (per sq arcsec)
* output
*	dzwtel	r4 detection zone area
* 2008 Keith Horne * SAAO
	logical firstcall/.true./
      if( firstcall ) then
	pi = 4. * atan2( 1., 1. )
* gaussian psf has effective area = gpix * fwhm^2
	gpix = pi / ( 2. * alog( 2. ) )
	firstcall = .false.
      end if
* magnification
	a = atotl( u )
* Vega flux 500 photons/cm2/s/A for I band
	fvega = 500. * cm2 * bw
* baseline flux
	f0 = fvega * 10.**( -0.4 * smag0 )
* source flux
	fs = abs( f0 / ( 1. + b0 ) )
* blend flux
	fb = abs( fs * b0 )
* square arcseconds
	sqas = gpix * fwhm * fwhm
* sky flux
	fsky = fvega * 10.**( -0.4 * skymag ) * sqas
* s/n ratio
	fnoise = fs * a + fb + fsky
	snr = fs * a * sqrt( t / fnoise )
	snrdet = 10.
* monte-carlo trials (increase for more accuracy)
	m = 5
* detection zone area
c	dzwtel = dzw( q, u, snrdet/snr, m )
* fast version
	dzwtel = dzwfast( q, u, 1./snr )

* report
      if( .false. ) then
	write(*,*) '** dzwtel(', q, u, snrdet/snr, ' )=', dzwtel
	write(*,*) '** A', a, ' m0', smag0, ' b0', b0
	write(*,*) '** fs', fs, ' fs A', fs*a, ' fb', fb, ' fsky', fsky
	write(*,*) '**', cm2/1e4, 'm^2', bw, 'A', t, 's  S/N', snr
      end if
	return
	end


*=================================================
	function dzwfast( q, u, dm )
* fast approx to detection area d(theta)d(lnu) on lens plane
* Input:
*	q	r4 planet/star mass ratio
*	u	r4 source-lens distance in lens plane (Re=1)
*	dm	r4 magnitude change for detection
* Output:
*	dzw	r4 planet detection zone area d(theta)d(lnu).
* 2002 Apr Keith Horne @ St.Andrews
* 2008 Jun KDH @ St.And - d(theta)d(lnu) metric
* 2008 Aug KDH @ St.And - fudge by factor pi/3.
	logical firstcall /.true./
      if( firstcall ) then
	pi = 4. * atan2(1.,1.)
	fudge = pi / 3.
	firstcall = .false.
      end if

	uu = u * u
	t = uu + 2.
	b = sqrt( uu * ( uu + 4. ) )

	fp = ( t + t ) / ( b * ( t + b ) )
	fm = 2. / b
	f = fp + fm

	dzwfast = fudge * q * f / dm
	return
	end
*=================================================
	function dzwplus( q, u, dm )
	twopi = 8. * atan2( 1., 1. )
	uu = u * u
	t = uu + 2.
	b = sqrt( uu * ( uu + 4. ) )
	f = t / ( b * ( t + b ) )
	dzwplus = ( q * twopi * f ) / ( 3. * dm )
	return
	end
*=================================================
	function dzwminus( q, u, dm )
	twopi = 8. * atan2( 1., 1. )
	uu = u * u
	t = uu + 2.
	b = sqrt( uu * ( uu + 4. ) )
	f = 1. / b
	dzwminus = ( q * twopi * f ) / ( 3. * dm )
	return
	end
*=================================================
	function dzw( q, u, dm, m )
* planet detection area d(theta)d(lnu) on lens plane
* Input:
*	q	r4 planet/star mass ratio
*	u	r4 source-lens distance in lens plane (Re=1)
*	dm	r4 magnitude change for detection
*	m	i4 number of monte-carlo trials (0=pixel centre)
* Output:
*	dzw	r4 planet detection zone area d(theta)d(lnu)
* 2002 Apr Keith Horne @ St.Andrews
* 2008 Jun KDH @ St.And - d(theta)d(lnu) metric
	dzw = dzw2( q, 0., u, dm, 0., 0., -1., m )
	return
	end
*=================================================
	function dzw1( q, xs, ys, dm, m )
* planet detection area d(theta)d(lnu) on lens plane
* Input:
*	q	r4 planet/star mass ratio
*	xs,ys	r4 source-lens coords in lens plane (Re=1)
*	dm	r4 magnitude change for detection
*	m	i4 number of monte-carlo trials (0=pixel centre)
* Output:
*	dzw1	r4 planet detection zone area d(theta)d(lnu).
* 2002 Apr Keith Horne @ St.Andrews
* 2008 Jun KDH @ St.And - d(theta)d(lnu) metric
	dzw1 = dzw2( q, xs, ys, dm, 0., 0., -1., m )
	return
	end
*=================================================
	function dzw2( q, xs, ys, dm, xsp, ysp, dmp, monte )
* planet detection area d(theta)d(lnu) on lens plane
* Input:
*	q	r4 planet/star mass ratio
*	xs,ys	r4 source-lens coords in lens plane (Re=1)
*	dm	r4 magnitude change for detection
*	xsp,ysp	r4 prior source-lens coords in lens plane (Re=1)
*	dmp	r4 prior accuracy (mag, <0 if no prior)
*	monte	i4 number of monte-carlo trials (0=pixel centre)
* Output:
*       dzw2     r4 change in detection zone area d(theta)d(lnu).
*
* - for q not small, dzw too big since 2 boxes overlap
* - for us or dm small, dzw too small since boxes too small
* 2002 Apr Keith Horne @ St.Andrews
* 2002 Apr KDH @ St.And - exploit symmetry
* 2002 Apr KDH @ St.And - Monte-Carlo sub-pixel sampling
* 2008 Jun KDH @ St.And - d(theta)d(lnu) metric
	wp = dzw2i( 1, q, xs, ys, dm, xsp, ysp, dmp, monte )
	wm = dzw2i( 2, q, xs, ys, dm, xsp, ysp, dmp, monte )
	dzw2 = wm + wp
	return
	end
*=================================================
	function dzw2i( im, q, xs, ys, dm, xsp, ysp, dmp, monte )
* planet detection area d(theta)d(lnu) on lens plane
* Input:
*	im	i4 image number (1=major, 2=minor)
*	q	r4 planet/star mass ratio
*	xs,ys	r4 source-lens coords in lens plane (Re=1)
*	dm	r4 magnitude change for detection
*	xsp,ysp	r4 prior source-lens coords in lens plane (Re=1)
*	dmp	r4 prior accuracy (mag, <0 if no prior)
*	monte	i4 number of monte-carlo trials (0=pixel centre)
* Output:
*       dzw2     r4 change in detection zone area d(theta)d(lnu).
*
* - for q not small, dzw too big since 2 boxes overlap
* - for us or dm small, dzw too small since boxes too small
* 2002 Apr Keith Horne @ St.Andrews
* 2002 Apr KDH @ St.And - exploit symmetry
* 2002 Apr KDH @ St.And - Monte-Carlo sub-pixel sampling
* 2008 Jun KDH @ St.And - d(theta)d(lnu) metric

	real*8 sum, sump, suma
	character*75 xlabel, ylabel, title
	parameter ( maxxy = 1e5 )
	real*4 pic( maxxy ), tr(6)
	logical firstcall / .true. /
      if( firstcall ) then
	iseed = 107176
	firstcall = .false.
      end if
	write(*,*) 'dzw2i   q', q, ' monte', monte, ' image', im
	write(*,*) 'New data sigma', dm, 'at x,y=', xs, ys
	write(*,*) 'Old data sigma', dmp, 'at x,y=', xsp, ysp

* sanity checks
	dzw2i = 0.
	if( q .le. 0. ) return
	if( dm .le. 0. ) return

* plots?
	iplot = 1
	iplot = 0

* planet einstein ring radius
	rp = sqrt( q )
* lens-source separation
	us = sqrt( xs*xs + ys*ys )
* lens-image distances
	u1 = uplus( us )
	u2 = uminus( us )
* image magnifications
	a1 = aplus( us )
	a2 = aminus( us )
* total magnification
	amp1 = a1 + a2

* prior
	prior = 0.

      if( dmp .gt. 0. ) then
	usp = sqrt( xsp*xsp + ysp*ysp )
	amp1p = atotl( usp )
	write(*,*) 'prior magnification:', amp1p
	u1p = uplus( usp )
	u2p = uminus( usp )
	delx = xsp - xs
	dely = ysp - ys
	cc = xs / us
	ss = ys / us
	delu = delx * cc + dely * ss
	delv = delx * ss - dely * cc
      end if

* sum over images
	nim = 2
	suma = 0.d0
c      do im=1,nim

* open plot
      if( iplot .gt. 0 ) then
      if( im .eq. 1 ) then
	call pgstart(1,2,new )
	csize = 1.5
	call pgsch( csize )
      end if
      end if


* lens-image separation and magnification
      if( im .eq. 1 ) then
	aim = a1
	ipm = 1
      else
	aim = a2
	ipm = -1
      end if

* detection zone grid
	nx = maxxy
	ny = 1
	call dzwgrid( q, us, dm, ipm, nx, ny, tr )
      if( nx .gt. 0 .and. ny .gt. 0 ) then
	dxdy = tr(2) * tr(6) - tr(3)*tr(5)
* why ?
	dxdy = dxdy * 2.

* detection zone map
	call dzwmap( q, us, nx, ny, tr, pic )

	nsum = 0
* planet search positions
	sump = 0.d0
	ixy = 0
      do iy=1,ny
      do ix=1,nx
	ixy = ixy + 1
	xp0 = tr(1) + tr(2) * ix + tr(3) * iy
	yp0 = tr(4) + tr(5) * ix + tr(6) * iy
* convert from ln-polar to cartesian
	radian = xp0
	r = exp( yp0 )
	xp0 = r * sin( radian )
	yp0 = r * cos( radian )

	pdet = 0.

* pixel centre
      if( monte .le. 0 ) then
	xp = xp0
	yp = yp0

* check for prior data
      if( dmp .gt. 0. ) then
	call lens2( delv, us + delu, xp, yp, q, amp2 )
	change =  2.5 * alog10( amp2 / amp1p ) 
	prior = abs( change / dmp )
      end if

      if( prior .lt. 1. ) then
c	call lens2( 0., us, xp, yp, q, amp2 )
c	change = 2.5 * alog10( amp2 / amp1 )
	change = pic( ixy )
	if( abs(change) .ge. dm ) pdet = 1.
      end if

* monte-carlo average over pixel
      else

	ntest = max( 1, monte )
	sum = 0.d0
	ndet = 0
	px0 = ix - 0.5
	py0 = iy - 0.5
      do i=1,ntest
	px = px0 + ran( iseed )
	py = py0 + ran( iseed )
	xp = tr(1) + tr(2) * px + tr(3) * py
	yp = tr(4) + tr(5) * px + tr(6) * py
* lnu,theta -> x,y
	radian = xp
	r = exp( yp )
	xp = r * sin( radian )
	yp = r * cos( radian )

	call lens2( 0., us, xp, yp, q, amp2 )
	change = 2.5 * alog10( amp2 / amp1 )
	sum = sum + abs( change )
	if( abs(change) .ge. dm ) ndet = ndet + 1
      end do
	pdet = float( ndet ) / ntest
	change = sum / ntest
      end if

* sum detections
	sump = sump + pdet

* fix pic for plot
      if( iplot .gt. 0 ) then
	pic(ixy) = abs( change )
      end if

* next ix,iy
      end do
      end do
* accumulate detection area
	add = sump * dxdy
	suma = suma + add
c	write(*,*) 'Image:', im, ' Ntest:', nx*ny,
c     &	' Nsum:', nsum, ' Area:', add
      end if

* plot
      if( iplot .gt. 0 ) then 
* window
	xmx = tr(1) + (nx+0.5) * tr(2)
	xmn = tr(1) - (nx-0.5) * tr(2)
	ymn = tr(4) + 0.5 * tr(6)
	ymx = ymn + ny * tr(6)
	write(*,*) 'pgenv(', xmn, xmx, ymn, ymx, ' )'
	call pgsci( 5 )
	call pgenv( xmn, xmx, ymn, ymx, 1, 0 )

* labels
      if( im .eq. 1 ) then
	title = 'Planet Detection Zones'
	xlabel = ' '
	ylabel = 'Major Image'
      else if( im .eq. 2 ) then
	ylabel = 'Minor Image'
	write( title, '( a, f7.5, a, f6.3, a, f6.2 )' )
     &	'q=', q,
     &	' U=', us,
     &	' A=', amp1

	dzw2i = suma
	write( xlabel, '( a, f6.1, a, f8.5 )' )
     &	' \gs:', dm*100., '%  DZW=', dzw2i
      end if	

	call pgsci(7)
	call pglabel( xlabel, ylabel, title )
	write( ylabel, '( A, F6.2, A, F8.5 )' )
     &	'A=', aim, '  DZW=', add
	call pgmtext( 'R', 1.5, 0.5, 0.5, ylabel )
* greyscale
	call mnmx( nx*ny, pic, fg, bg )
	write(*,*) fg, bg
	fg = 0.
	bg = 2. * dm
      do i=1,2
	call pggray( pic, nx,ny, 1,nx, 1,ny, bg, fg, tr)
	call pgcont( pic, nx,ny, 1,nx, 1,ny, dm, 1, tr )
	tr(2) = - tr(2)
	tr(1) = tr(1) - tr(2)
      end do
	call pgpoint( 1, xim, yim, 2 )
	call pgsci( 7 )
      end if

      if( iplot .gt. 0 ) then
      if( im .eq. nim ) then
	call pgstop
      end if
      end if

* next image
c      end do

* area where planet was detected
	dzw2i = suma
	return
	end

*=================================================
	subroutine dzwgrid( q, u, dm, ipm, nx, ny, tr )
* detection zone grid
* Input:
*	q	r4 planet/star mass ratio (>0)
*	u	r4 source-lens separation (1=Re)
*	dm	r4 mag change (>0)
*	ipm	i4 >0 for major image, <0 for minor image
*	nx,ny	i4 max pixels (x=theta, y=lnu)
* Output:
*	nx,ny	i4 grid dimensions (0 if failed)
*	tr(6)	r4 grid transformation matrix
* 2002 May Keith Horne @ St.Andrews
* 2008 Jun KDH @ St.And - d(theta)d(lnu) grid
	real*4 tr(6)
	real*4 uold/0./
	real*8 uu, top, bot
	pi = 4. * atan2(1.,1.)
	maxx = nx
	maxy = ny
	maxxy = maxx * maxy
	if( maxxy .le. 0 ) return
	nx = 0
	ny = 0
	if( dm .le. 0. ) return
	if( q .le. 0. ) return
* planet einstein ring radius
	rp = sqrt( q )
      if( u .ne. uold ) then
	uold = u
	uu = u * u
	top = uu + 2.d0
	bot = dsqrt( uu * ( uu + 4.d0 ) )
	uu1 = ( top + bot ) * 0.5d0
	uu2 = ( top - bot ) * 0.5d0
* image magnifications
	a1 = uu1 / bot
	a2 = uu2 / bot
* lens-image distances
	u1 = sqrt( uu1 )
	u2 = sqrt( uu2 )
* total magnification
	a = a1 + a2
	write(*,*) 'dzwgrid: u', u, ' u+', u1, ' u-', u2
	write(*,*) 'dzwgrid: A', a, ' A+', a1, ' A-', a2
      end if

	qdm = q / dm

* image position and magnification
      if( ipm .gt. 0 ) then
	ui = u1
	ai = a1
	oi = ( a1 + a1 - 1. ) * qdm / uu1
      else
	ui = u2
	ai = a2
	oi = ( a2 + a2 ) * qdm / uu2
      end if
	if( ai .le. 0. ) return

* box size
	up = sqrt( ai * max( 1., 1./dm ) )
	rtest = rp * up
	omega = ( a + a - 1. ) * qdm
	size = sqrt( omega )
	si = sqrt( oi )
	rtest = size
	write(*,*) 'dzwgrid im', ipm, ' u', ui, ' a', ai, ' omega', oi

      if( ipm .gt. 0 ) then
	fudge = sqrt( 1. + sqrt( dm/0.1 ) )
	xpand = 0.8 * fudge
	ypand = 0.8 * fudge + 0.2 * ( fudge - 1. )
      else
	fudge = sqrt( ( oi + oi ) * u2 / omega )
	xpand = 0.7 * fudge * max( 1., sqrt( dm/0.1 ) )
	ypand = 0.6 * fudge
      end if
	xbox = xpand * rtest
	ybox = ypand * rtest

	factor = exp( abs( ybox / ui ) )
	yout = ui * factor
	yin = ui / factor

	if( ui * yin .lt. 0. ) yin = 0.
	xmx = max( xbox, 0.5 * abs( yout - ui ) )
	ymn = min( yout, yin )
	ymx = max( yout, yin )
	
* change variables to x -> theta and y -> lnu
* range in y = lnu
	
	add = abs( ybox / ui )
	uln = alog( ui )
	ymn = uln - add
	ymx = uln + add
* half-range in x = theta
	xmn = 0.
	if( ipm .lt. 0 ) xmn = pi
	xmx = min( pi*0.5, max( xbox / ui, add ) )

* split box into pixels ( x-> theta, y-> lnu )
c	stepx = min( rp, rtest ) / 10.
c	stepx = min( rp, rtest, si ) / 10.
	div = 50. * ui
	step1 = rp / div
	step2 = rtest / div
	step3 = si / div
	write(*,*) 'im', im, ' size/rp', 1., step2 / step1, step3 / step1
	stepx = min( rp / 2./ ui, step2, step3 )
	stepy = stepx

	nx = nint( xmx / stepx )
	ny = nint( ( ymx - ymn ) / stepy )
	nxy = nx * ny
	factor = float( maxxy ) / max( 1, nxy )
	maxx = nx * factor
	maxy = ny * factor
	minx = 10
	miny = 10
	nx = max( minx, min( maxx, nx ) )
	ny = max( miny, min( maxy, ny ) )
	stepx = xmx / nx
	stepy = ( ymx - ymn ) / ny

* transform matrix
	x0 = xmn - 0.5 * stepx
	y0 = ymn - 0.5 * stepy
	tr(6) = stepy
	tr(5) = 0.
	tr(4) = y0
	tr(3) = 0.
	tr(2) = stepx
	tr(1) = x0
* report
	write(*,*) 'dzwgrid: u', u, ' delta', dm, ' q', q, ' im', ipm
	write(*,*) 'dzwgrid: x=theta', nx, stepx,
     &		tr(1)+tr(2), tr(1)+nx*tr(2)
	write(*,*) 'dzwgrid: y=ln(u)', ny, stepy,
     &		tr(4)+tr(6), ui, tr(4)+ny*tr(6)
	return
	end


*=================================================
	subroutine dzwmap( q, u, nx, ny, tr, dm )
* detection zone map on (theta,lnu) grid of planet positions
* Input:
*	q	r4 planet/star mass ratio (>0)
*	u	r4 source-lens separation (1=Re)
*	nx,ny	i4 grid dimensions (x=theta,y=lnu)
*	tr(6)	r4 grid transformation matrix
* Output:
*	dm(nx*ny) r4 mag changes
* 2002 May Keith Horne @ St.Andrews
* 2008 Jun KDH @ St-And x,y > theta, lnu
	real*4 dm(*)
	real*4 tr(6)
	a = atotl( u )
* planet search positions
	ixy = 0
      do iy=1,ny
      do ix=1,nx
	ixy = ixy + 1
	radian = tr(1) + tr(2) * ix + tr(3) * iy
	rln = tr(4) + tr(5) * ix + tr(6) * iy
	r = exp( rln )
	x = r * sin( radian )
	y = r * cos( radian )
	call lens2( 0., u, x, y, q, a2 )
	dm( ixy ) = 2.5 * alog10( a2 / a )
      end do
      end do
	return
	end
