PRO hsi_vis_mpfit_plotfit, visin, srcstr, mapcenter, rchisq, $
   MAP_RESID=map_resid, EBPLOT=ebplot, TIME=time,$
   PLOT=plot,_EXTRA=extra
;
; Generates comparison plot between input visibility data and model visibilities
;
; visin is an array of visibility structures
; srcstr is an array of source component structures
; /MAP_RESID also generates back projection map of residuals.
;
; 10-Nov-05 	Initial working version, broken out from hsi_vis_fwdfit. (ghurford@ssl.berkeley.edu)
; 13-Nov-05 gh	Adapt to revised definition of srcstr.xyoffset
; 14-Nov-05 gh	Add MAP_RESID keyword
;  9-Dec-05 gh	Minor changes.
; 21-Dec-05 gh	Add EBPLOT keyword to optionally plot error bars in blue.
; 22-dec-05 gh	Suppress amplitude and residual points for which sigamp > amplitude
; 26-Jan-06 ejs Added optional rchisq variable (reduced chi^2)
; 01-mar-06 ejs Added plot command (default=1) to prevent plotting (plot=0)
; 07-mar-06 ejs Added time keyword to put time on map; added _EXTRA keyword
; 28-Mar-06 gh  Reimplement fit plot colors with thicker lines
; june 2007 ejs  Revised name and calls for mpfit compatibility

DEFAULT, plot, 1
DEFAULT, ebplot, 1
if keyword_set(PLOT) then !p.multi 	= [0,1,1]
nvis    	= N_ELEMENTS(visin)
npt     	= 2*nvis
jdum 		= FINDGEN(npt)                    		     ; dummy 'x' values used in fitting routine
srcparm   	= hsi_vis_fwdfit_structure2array(srcstr, mapcenter)
ampobs 		= ABS(visin.obsvis)
;
			; normal
    visx    = FLOAT(visin.obsvis)
    visy    = IMAGINARY(visin.obsvis)
; Will the following clobber the subsequent plot commands?
IF KEYWORD_SET(nophase) THEN BEGIN                       		; Set phase to zero if /NOPHASE is set
    visx    = ABS(visin.obsvis)
    visy    = FLTARR(nvis)
ENDIF
;
dummy 		= hsi_vis_select (visin, PAOUT=paout)				; calculate paout
scpa        = visin.isc+1 + (paout/180. MOD 1) 		       		; (sc# + pa/180) is used as abscissa for plots
i 			= SORT(scpa)
j 			= WHERE(ampobs GT visin.sigamp, nok)				; cases where amplitude > its error
;
; Default option is to generate fit plots.
IF keyword_set(PLOT) then begin
    IF nok LE 1 THEN j = i											; Use j to just plot good amplitudes
    PLOT,  scpa[j], ampobs[j], XTITLE='SC + PA/180.',      XRANGE=[1,10], XSTYLE=1, $
                     YTITLE='ph / cm2 / sec',    /YLOG,  PSYM=7,  SYMSIZE=0.6, THICK=2, $
                      TITLE="Observed (x), Fitted (-) and Difference (!9B!3) Amplitudes"
    for kgrid=1,9 do oplot,[kgrid,kgrid],[0.1,1000],line=1 ; plot grid separators

    IF KEYWORD_SET(ebplot) NE 0 THEN BEGIN
        ylim = 1.02*10.^!Y.CRANGE[0]			; a minimum value that will just fit on plot
	    LOADCT, 6, /SILENT ; loadct,6 produces a blank screen (ejs)
	    ERRPLOT, scpa[i], (ampobs[i]-visin[i].sigamp > ylim), ampobs[i]+visin[i].sigamp, width=0, COLOR=192
    ENDIF
loadct,6
; Overplot 'fitted' amplitudes from source model
    x=[[visin.u],[visin.v]]
    y=fltarr(npt)
    err=y+1.0
    visxyfit    = hsi_vis_mpfit_func(srcparm,x=x,Y=y,ERR=err,pa=pa,mapcenter=mapcenter)
    ampfit      = SQRT(visxyfit[0:nvis-1]^2 + visxyfit[nvis:*]^2)
    IF KEYWORD_SET(ebplot) NE 0 then OPLOT, scpa[i], ampfit[i], PSYM=10,THICK=2,COLOR=64
message,'',/info
 vfitx=visxyfit[0:nvis-1]
 vfity=visxyfit[nvis:*]
 ;window,0 & plot,vfitx,vfity,tit='vfitx vs vfity',/iso
 ;window,1 & plot,visx,visy,tit='visx vs visy',/iso
loadct,0
; Overplot amplitude of difference between fitted and observed visibilities
    visdiff     = SQRT((visx-visxyfit[0:nvis-1])^2 + (visy-visxyfit[nvis:*])^2)
    IF KEYWORD_SET(ebplot) NE 0 then OPLOT, scpa[j], visdiff[j], PSYM=6, SYMSIZE=0.4, THICK=2, COLOR=128
    LOADCT, 0, /SILENT
ENDIF

message,'Plotting finished',/info



;
IF N_ELEMENTS(RCHISQ) NE 0 THEN rchisq=total((visdiff/visin.sigamp)^2)/(npt-8) ; reduced chi^2
;
; Optionally plot residual map
IF KEYWORD_SET(map_resid) EQ 0 THEN RETURN
	visdiff 		= visin
	visfit 			= COMPLEX(visxyfit[0:nvis-1], visxyfit[nvis:*])
	visdiff.obsvis 	= visin.obsvis - visfit
	hsi_vis_bpmap, visdiff,time=time,_EXTRA=extra


RETURN
END