;------------------------------------------------------------------
;+
; NAME:
;	NMOSAIC
; PURPOSE:
;
;       This program adds an image to a mosaic of BBSO VMG
;       images, and also adds the corresponding J image to
;       a J mosaic.  It is designed for polar images.
;
; CATEGORY:
;	POLAR
; CALLING SEQUENCE:
;	nmosaic,stackimage,newimage,nhdr,imgscl,jstack,jimage
; INPUTS:
;	stackimage = image containing composite map.
;	             A "blank" stack can be created with the command
;	             newstack = REPLICATE(-32000,xsize,ysize)
;                    xsize = 4000, ysize=2000
;	newimage   = image to add to stack.
;       nhdr       = FITS header of new image.
;       imgscl     = image scale (pix/arcsec).
;       jstack     = J mosaic image (xsize=4000, ysize=2000)
;       jimage     = J image to add.
; KEYWORD PARAMETERS:
; OUTPUTS:
;	The stackimage and jstack are modified.
; COMMON BLOCKS:
; NOTES:
; MODIFICATION HISTORY:
;	J. Varsik,	08 March 1999
;       J. Varsik,      19 April 1999  Changed mosaic size to
;       3000x1500.
;       J. Varsik,      06 May   1999  Changed back to 4000x2000.
;-
;-------------------------------------------------------------------

pro nmosaic,stack,new,nhdr,imgscl,jstack,jimage,help=help

    ; Display IDL header if help is required.
    IF (KEYWORD_SET(help)) THEN BEGIN
      GET_IDLHDR,'nmosaic.pro'
      GOTO,finishup
    ENDIF

    TRUE = 1
    FALSE = 0

    ; Obtain initial guess of position of new image in mosaic.
    ew = FXPAR(nhdr,'CRVAL1')
    ns = FXPAR(nhdr,'CRVAL2')
    h = FIX(imgscl * ew + 256)
    k = FIX((-1.0 * imgscl * ns) + 192)
;   h and k are the position of the center of the Sun in picture coordinates.
    IF (k LT 192) THEN npole = TRUE
    IF (k GE 192) THEN npole = FALSE

    goodnew = new[16:505,4:375]
    goodj   = jimage[16:505,4:375]
    ; edges of images are bad from flatfield and VMG system
    mxmin = 2000 - h + 16 
    mxmax = 2000 - h + 505
     ; Find the initial coordinates of the new picture in the mosaic
    IF (npole) THEN BEGIN
        PRINT,'North pole image'
        mymin = -k + 4
        mymax = -k + 375
    ENDIF ELSE BEGIN
        PRINT,'South pole image'
        mymin = 2000 - k + 4
        mymax = 2000 - k + 375
    ENDELSE
    ; Have we put any pictures on this mosaic?
    stindex = WHERE(stack GT -32000, ncount)
    IF (ncount EQ 0) THEN BEGIN
        ; If this is the first picture, just place it and quit.
        stack[mxmin:mxmax, mymin:mymax] = goodnew
        jstack[mxmin:mxmax, mymin:mymax] = goodj
    ENDIF ELSE BEGIN
        ; Otherwise,
        ; Extract the current image from the mosaic that is where we think
        ; the new picture will go.
        PRINT,mxmin,mxmax,mymin,mymax
        mextract = stack[mxmin:mxmax,mymin:mymax]
        ; Find the estimated overlap area.
            ; Steps to do this:
            ; a. find the largest rectangle in mextract that is not -32000
            ; ** Simplest way to do this now might be to display
            ; mextract and have user click on
            ; lower left and upper right
            ; of an appropriate rectangle
        WINDOW,/FREE,XSIZE=489,YSIZE=371
        mywin = !D.WINDOW  ; save the index of my new window
        TV,BYTSCL(mextract,MIN=-30,MAX=30)
        CURSOR,xlo,ylo,/UP,/DEVICE
        PLOTS,xlo,ylo,PSYM=1,/DEVICE
        CURSOR,xhi,yhi,/UP,/DEVICE
        PLOTS,xhi,yhi,PSYM=1,/DEVICE
        PRINT,xlo,xhi,ylo,yhi
        ; Draw the rectangle
        PLOTS,[xlo,xhi],[ylo,ylo],/DEVICE
        PLOTS,[xhi,xhi],[ylo,yhi],/DEVICE
        PLOTS,[xhi,xlo],[yhi,yhi],/DEVICE
        PLOTS,[xlo,xlo],[yhi,ylo],/DEVICE
            ; b. get mx2a (which is that rectangle from mextract) and
            ;    mx2b (which is that rectangle from goodnew)
        mx2a = mextract[xlo:xhi,ylo:yhi]
        mx2b = goodnew[xlo:xhi,ylo:yhi]
        ; Use align.pro to find the correction to the initial position.
        algnres = ALIGN(mx2a,mx2b)
        PRINT,algnres
        ; Apply the correction and place the new picture in the mosaic.
        mxmin = mxmin + ROUND(algnres[0])
        mxmax = mxmin + 489
        mymin = mymin + ROUND(algnres[1])
        mymax = mymin + 371
        stack[mxmin:mxmax, mymin:mymax] = goodnew
        ; Now we need to adjust the new jimage to match the existing
        ; J mosaic.
        ;rtest = ROUND(850.0 * imgscl) ;select points at 850 arcsec radius.
        ; Maybe we don't really need to do this.
        jstack[mxmin:mxmax, mymin:mymax] = goodj
        WDELETE,mywin
    ENDELSE

    finishup:

    return
end