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