pro addradcomps, img, radprof, outimg,  $
                 inimg = inimg, $
                 nsigma = NSIGMA, pa = PA,  $
                 minrad = minrad, maxrad = maxrad, drad = drad,  $
                 review = review

;+
; given an unsaturates image and its radial noise profile (as
; generated by AORADNORM.PRO), insert artificial companions
; at multiple radii.  flux scaling is based on peak pixel of the input
; image (admittedly very crude).
; 
; this routine provides a method for verifying the results 
; companion limits extracted from AORADNORM.PRO.
;
; star doesn't have to be centered on input image for this to work.
; however, if 'inimg' is passed (used for saturated images), then it
; has to be aligned with the input image.
;
; 09/28/01 M. Liu
;-


if not(keyword_set(NSIGMA)) then NSIGMA = 10.0  ; detection limit 
if not(keyword_set(minrad)) then minrad = 10  ; min radius for comps (pix)
if not(keyword_set(maxrad)) then maxrad = 90  ; max radius for comps (pix)
if not(keyword_set(drad)) then drad = 10      ; spacing in radius
if not(keyword_set(PA)) then PA = 30.  ; CCW from +x axis, in degrees


if n_params() ne 3 then begin
    print, 'pro addradcomps, img, radprof, outimg,'
    print, '                 [nsigma=', strc(NSIGMA), '[pa=', strc(PA), '], '
    print, '                 [minrad=', strc(minrad), ', [maxrad=', strc(maxrad), '], '+ $
      '[dra = ', strc(drad), ']'
    return
endif


; initialize
outimg = img
rr = makex(minrad, maxrad, drad)

; assign image to use for adding companions
; (normally this is just the input image.  but for saturated images
; you can pass an appropriately scale unsaturated image)
if not(keyword_set(inimg)) then inimg = img

; compute flux ratio curve for companion
dflux = max(inimg)/radprof(*, 2)/NSIGMA

; add companions
for i = 0, n_elements(rr)-1 do begin
    df = interpol(dflux, radprof(*, 0), rr(i))
    x = rr(i)*cos(PA*!dtor)
    y = rr(i)*sin(PA*!dtor)
    outimg = outimg + shift(inimg/df, x, y)
    print, i, rr(i), x, y, df
endfor


if keyword_set(review) then begin
    ; normalize (assume star is centered for convenience)
    aoradnorm, outimg, norm, /nodisp,  $
      x = (size(outimg))(1)/2., y = (size(outimg))(2)/2.

    loadct, 15
    win, 0
    display2, img, tit = 'input image'
    win, 1
    display2, outimg, tit = 'output image'
    win, 2
    display2, norm > (-3) < 5, /tv, tit = 'noise-norm'
    if not(getyn('ADDRADCOMPS: satified with added companions?')) then stop
    wdelete, 0
    wdelete, 1
    wdelete, 2    
endif

end