pro chi2_profile,visfit,visc

;number of sources fitted
n = n_elements(visfit)

type = visfit[0].srctype
  
if n gt 3 then begin
  print, 'CANNOT FIT MORE THAN THREE SOURCES, RETURNING'
  return
endif

xyoffset = visc[0].xyoffset
visin = visfit

;make chi^2 plots of source parameters (SRCX, SRCY, Flux, FWHM)
;code built upon program of Ed Schmahl

;initialize parameter arrays
;11 is for 5 data points on each side of the fitted parameter
fwhm = dblarr(n,11)
flux = dblarr(n,11)
srcx = dblarr(n,11)
srcy = dblarr(n,11)
eccen = dblarr(n,11)
srcpa = dblarr(n,11)
loop = dblarr(n,11)

;fitted parameters, calculated in hsi_vis_mpfit
fit_fwhm = dblarr(n)
fit_flux = dblarr(n)
fit_srcx = dblarr(n)
fit_srcy = dblarr(n)
fit_eccen = dblarr(n)
fit_srcpa = dblarr(n)
fit_loop = dblarr(n)

;create "chi^2" space arrays and assign fitted parameters
for ctr=0,n-1 do begin
  fwhm[ctr,*]=visin[ctr].srcfwhm+(findgen(11)-5.)
  flux[ctr,*]=visin[ctr].srcflux+(findgen(11)-5.)
  srcx[ctr,*]=visin[ctr].srcx+(findgen(11)-5.)
  srcy[ctr,*]=visin[ctr].srcy+(findgen(11)-5.)
  eccen[ctr,*]=visin[ctr].eccen+((findgen(11)-5.)/100.)
  srcpa[ctr,*]=visin[ctr].srcpa+(findgen(11)-5.)
  loop[ctr,*]=visin[ctr].loop_angle+(findgen(11)-5.)
  
  fit_fwhm[ctr] = visin[ctr].srcfwhm
  fit_flux[ctr] = visin[ctr].srcflux
  fit_srcx[ctr] = visin[ctr].srcx
  fit_srcy[ctr] = visin[ctr].srcy
  fit_eccen[ctr] = visin[ctr].eccen
  fit_srcpa[ctr] = visin[ctr].srcpa
  fit_loop[ctr] = visin[ctr].loop_angle
end

;"chi^2" parameter...reduced chi^2 arrays for each of the 11 data points
;in chi^2 space, TBD in hsi_vis_mpfit_plotfit
c_fwhm = dblarr(n,11)
c_flux = dblarr(n,11)
c_srcx = dblarr(n,11)
c_srcy = dblarr(n,11)
c_eccen = dblarr(n,11)
c_srcpa = dblarr(n,11)
c_loop = dblarr(n,11)

;initialize reduced chi^2?
rchisq=-1

;check device name/change to PS
save_dev = !d.name
set_plot,'ps'

dname = !d.name
set_plot,'Z' ; send plots to z buffer instead of screen or printer

;calculate reduced chi^2 for every data point in chi^2 space

;fluxes
visin = visfit
for ctr1=0,n-1 do begin ;loop over number of sources
  for ctr2=0,10 do begin  ;loop over data points in chi^2 space
    visin[ctr1].srcflux = flux[ctr1,ctr2]
    hsi_vis_mpfit_plotfit,visc,visin,xyoffset,rchisq
    c_flux[ctr1,ctr2] = rchisq
  endfor
endfor

;sizes
visin = visfit
for ctr1=0,n-1 do begin ;loop over number of sources
  for ctr2=0,10 do begin  ;loop over data points in chi^2 space
    visin[ctr1].srcfwhm = fwhm[ctr1,ctr2]
    hsi_vis_mpfit_plotfit,visc,visin,xyoffset,rchisq
    c_fwhm[ctr1,ctr2] = rchisq
  endfor
endfor

;x position
visin = visfit
for ctr1=0,n-1 do begin ;loop over number of sources
  for ctr2=0,10 do begin  ;loop over data points in chi^2 space
    visin[ctr1].srcx = srcx[ctr1,ctr2]
    hsi_vis_mpfit_plotfit,visc,visin,xyoffset,rchisq
    c_srcx[ctr1,ctr2] = rchisq
  endfor
endfor

;y position
visin = visfit
for ctr1=0,n-1 do begin ;loop over number of sources
  for ctr2=0,10 do begin  ;loop over data points in chi^2 space
    visin[ctr1].srcy = srcy[ctr1,ctr2]
    hsi_vis_mpfit_plotfit,visc,visin,xyoffset,rchisq
    c_srcy[ctr1,ctr2] = rchisq
  endfor
endfor

if type ne 'circle' then begin
  ;eccentricity
  for ctr1=0,n-1 do begin ;loop over number of sources
    for ctr2=0,10 do begin  ;loop over data points in chi^2 space
      visin[ctr1].eccen = eccen[ctr1,ctr2]
      hsi_vis_mpfit_plotfit,visc,visin,xyoffset,rchisq
      c_eccen[ctr1,ctr2] = rchisq
    endfor
  endfor

  ;position angle (degrees East of North)
  for ctr1=0,n-1 do begin ;loop over number of sources
    for ctr2=0,10 do begin  ;loop over data points in chi^2 space
      visin[ctr1].srcpa = srcpa[ctr1,ctr2]
      hsi_vis_mpfit_plotfit,visc,visin,xyoffset,rchisq
      c_srcpa[ctr1,ctr2] = rchisq
    endfor
  endfor
  
  if type eq 'loop' then begin
    ;loop angle
    for ctr1=0,n-1 do begin ;loop over number of sources
      for ctr2=0,10 do begin  ;loop over data points in chi^2 space
        visin[ctr1].loop_angle = loop[ctr1,ctr2]
        hsi_vis_mpfit_plotfit,visc,visin,xyoffset,rchisq
        c_loop[ctr1,ctr2] = rchisq
      endfor
    endfor
  endif

endif

time=visc[0].trange
date=anytim(time(0),/ccsds,/date)
set_plot,dname   ; restore plot device
if dname eq 'X' then window,/free,xsize=512,ysize=768

energy = visc[0].erange
en1 = string(fix(energy[0]))
en2 = string(fix(energy[1]))
energy_string = en1+'-'+en2
energy_string=strcompress(energy_string,/rem)
psname='chisq_profiles'+date+'_E'+energy_string+'.ps
if !d.name eq 'PS' then device,/port,xsize=16,ysize=24,yoffset=1,file=psname

if type eq 'circle' then nparam = 4 else begin
  if type eq 'ellipse' then nparam = 6 else nparam = 7
endelse
!p.multi=[0,n,nparam,0,1]
datim=anytim(time[0],/ecs)
en_lbl = string(energy[0],energy[1],format='(f5.1,1x,"-",1x,f5.1,1x,"keV",1x)')
date_en = datim+' '+en_lbl

!x.style=1 & !y.style=1
!p.charsize=1.5
yt='CHISQUARE'

for ctr=0,n-1 do begin
  t1 = '!7v!e2!n!3 PROFILE FOR '
  t2 = 'FWHM'+strcompress(ctr+1,/rem)
  xt = t2+' (arcsec)'
  plot,fwhm[ctr,*],c_fwhm[ctr,*],psym=-6,title=t1+t2,ytitle=yt,xtitle=xt
  legend,['fit_fwhm='+strcompress(fit_fwhm[ctr],/rem)],charsize=.75

  t2 = 'FLUX'+strcompress(ctr+1,/rem)
  xt = t2+' (photon cm!e-2!n s!e-1!n)'
  plot,flux[ctr,*],c_flux[ctr,*],psym=-6,title=t1+t2,ytitle=yt,xtitle=xt
  legend,['fit_flux='+strcompress(fit_flux[ctr],/rem)],charsize=.75

  t2 = 'SRCX'+strcompress(ctr+1,/rem)
  xt = t2+' (arcsec)'
  plot,srcx[ctr,*],c_srcx[ctr,*],psym=-6,title=t1+t2,ytitle=yt,xtitle=xt
  legend,['fit_srcx='+strcompress(fit_srcx[ctr],/rem)],charsize=.75

  t2 = 'SRCY'+strcompress(ctr+1,/rem)
  xt = t2+' (arcsec)'
  plot,srcy[ctr,*],c_srcy[ctr,*],psym=-6,title=t1+t2,ytitle=yt,xtitle=xt
  legend,['fit_srcy='+strcompress(fit_srcy[ctr],/rem)],charsize=.75
  
  if type ne 'circle' then begin
    t2 = 'Eccentricity'
    xt = t2+' (arcsec)'
    plot,eccen[ctr,*],c_eccen[ctr,*],psym=-6,title=t1+t2,ytitle=yt,xtitle=xt
    legend,['fit_eccen='+strcompress(fit_eccen[ctr],/rem)],charsize=.75

    t2 = 'Position Angle'
    xt = t2+' (Degrees East of North)'
    plot,srcpa[ctr,*],c_srcpa[ctr,*],psym=-6,title=t1+t2,ytitle=yt,xtitle=xt
    legend,['fit_srcpa='+strcompress(fit_srcpa[ctr],/rem)],charsize=.75
    
    if type eq 'loop' then begin
      t2 = 'Loop Angle'
      xt = t2+' (Degrees)'
      plot,loop[ctr,*],c_loop[ctr,*],psym=-6,title=t1+t2,ytitle=yt,xtitle=xt
      legend,['fit_loop='+strcompress(fit_loop[ctr],/rem)],charsize=.75
    endif
  endif
endfor

;flare data as a "universal" y-axis title
xyouts,0.015,0.5,date_en,orient=90,align=0.5,/norm,charsize=1

if !d.name eq 'PS' then device,/close
if !d.name eq 'PS' then spawn,'gv '+psname

!p.multi=0
set_plot,save_dev
end