pro visfit,visin,srcparm=srcparm,event=event,$
      verbose=verbose,$
       spawn_gs=spawn_gs,$
       ZBUFFER=zbuffer

; Performs a vis fwdfit and outputs a parameter file named using time & energy

; srcparm examples are given below
; event can be specified instead of srcparm
; input srcparm positions are relative to sun center
; they are converted to xoffset-relative by this program

; Example of finding a double-Gaussian fit of double footpoint vis 
; OPTIONAL KEYWORDS
; set verbose to get screen printouts of chisq slices
; Set spawn_gs to show ps plots by spwaning ghostscript
 
; Create the starting values for [x1,y1,x2,y2,flux1,flux1,fwhm1,fwhm2]

; example visin for 25-50 keV is found in
; /home/schmahl/3jun17/vis_savefiles/vis3jun17_22:53:00.sav'
; for 30-40 keV:


default,event,13
default,verbose,0
default,spawn_gs,0

nsrc=2
srcparm=fltarr(10,nsrc)
; Same array as used in gh vispkg, HSI_VIS_FWDFIT_ARRAY2STRUCTURE
; srcparm(*,0)=[index,srcflux,srcx-mapctr,srcy-mapctr,srcfwhm,ecmsr1,ecsrm2,loop_angle,albedo_ratio,srcheight]

vise=hsi_vis_edit(visin)
message,'min(vise.sigamp)='+string(min(vise.sigamp)),/info
visc=hsi_vis_combine(vise,/conjugate)
message,'min(visc.sigamp)='+string(min(visc.sigamp)),/info


  ; starting values

case event of
1:srcparm=[[1, 4.38, -839.0,-409.0, 6.0, 0, 0, 0, 0, 0],$
           [1, 11.0, -831.0,-423.0, 7.0, 0, 0, 0, 0, 0]]
2:srcparm=[[1, 67.,   430.0, 128.0, 1.0, 0, 0, 0, 0, 0],$
           [1, 68,5., 430.5, 152.5,15.0, 0, 0, 0, 0, 0]]
3:srcparm=[[1, 2.5,  -345.5, 149.7, 6.0, 0, 0, 0, 0, 0],$
           [1, 7.2,  -330.9, 107.1, 7.0, 0, 0, 0, 0, 0]]
4:srcparm=[[1, 29.,   790.6,-251.7, 5.5, 0, 0, 0, 0, 0],$
           [1, 27.5,  814.0,-225.0, 8.0, 0, 0, 0, 0, 0]]
5:srcparm=[[1, 8.8,  -664.0,-247.5, 7.0, 0, 0, 0, 0, 0],$
           [1, 4.5,  -654.2,-267.5, 7.0, 0, 0, 0, 0, 0]]
6:srcparm=[[1, 20.5,  636.0,-269.5, 6.0, 0, 0, 0, 0, 0],$
           [1, 15.,   611.0,-266.0, 7.0, 0, 0, 0, 0, 0]]
7:srcparm=[[1, 32.8, -912.0,-199.0, 9.0, 0, 0, 0, 0, 0],$
           [1, 11.0, -918.0,-177.0, 7.0, 0, 0, 0, 0, 0]]
8:srcparm=[[1, 48.0, -798.0,-148.0, 5.0, 0, 0, 0, 0, 0],$ ; jun 17, 2003
           [1, 36.0, -782.0,-132.0, 5.0, 0, 0, 0, 0, 0]]
9:srcparm=[[1, 33.0, -912.0,-199.0, 5.0, 0, 0, 0, 0, 0],$  ; sep 8, 2002
           [1, 11.0, -918.0,-177.0, 5.0, 0, 0, 0, 0, 0]]
10:srcparm=[[1, 33.0,-873.0,-244.0, 3.0, 0, 0, 0, 0, 0],$  ; jul 23, 2002
           [1, 11.0, -871.0,-216.0, 3.0, 0, 0, 0, 0, 0]]
11:srcparm=[[1, 33.0, -96.0, -20.0, 3.0, 0, 0, 0, 0, 0],$  ; nov 19, 2003
           [1, 11.0,  -90.0, -10.0, 3.0, 0, 0, 0, 0, 0]]
12:srcparm=[[1, 127.0,-827.0, 421.0, 6.0, 0, 0, 0, 0, 0],$  ; 2005, aug 25  043715
           [1,  54.0, -821.0, 429.0, 3.0, 0, 0, 0, 0, 0]]
13:srcparm=[[1, 14.6,-96.2, -21.0, 4.0, 0, 0, 0, 0, 0],$  ; 2003, dec 01 035928
           [1,  12.5,-89.7, -10.7, 3.0, 0, 0, 0, 0, 0]]
14:srcparm=[[1, 21.0,  114.0, 298.0, 2.0, 0, 0, 0, 0, 0],$  ; 2005, jan 15  223230
           [1,  34.7,  141.0, 309.0, 2.0, 0, 0, 0, 0, 0]]
endcase

xyoffset=visin(0).xyoffset
;   [(srcparm(2,0)+srcparm(2,1))/2.,(srcparm(3,0)+srcparm(3,1))/2.]
srcparm(2,0:1)=srcparm(2,0:1)-xyoffset[0]
srcparm(3,0:1)=srcparm(3,0:1)-xyoffset[1]

;stop

; Try a fit with starting values

src_struc=hsi_vis_fwdfit_array2structure(srcparm,xyoffset)
if verbose then print,'srcparm starting values:
if verbose then print,srcparm

hsi_vis_fwdfit,visc,/multi,srcout=srcfit1,srcin=src_struc,maxiter=1500
  ; As a check, use last outputs as new inputs
srcparm1=hsi_vis_fwdfit_structure2array(srcfit1,xyoffset)

if verbose then print,'New srcparm starting values:
if verbose then print,srcparm1
fitstddev=1
hsi_vis_fwdfit,visc,/multi,srcout=srcfit2,srcin=srcfit1,maxiter=1500,FITSTDDEV=fitstddev

eflux1=fitstddev[0].srcflux & eflux2=fitstddev[1].srcflux 
efwhm1=fitstddev[0].srcfwhm & efwhm2=fitstddev[1].srcfwhm 
esrcy1=fitstddev[0].srcy    & esrcy2=fitstddev[1].srcy    
esrcx1=fitstddev[0].srcx    & esrcx2=fitstddev[1].srcx    

;stop

fwhm1 = srcfit2[0].srcfwhm + (findgen(11)-5.)
fwhm2 = srcfit2[1].srcfwhm + (findgen(11)-5.)
srcx1 = srcfit2[0].srcx + (findgen(11)-5) 
srcy1 = srcfit2[0].srcy + (findgen(11)-5)
srcx2 = srcfit2[1].srcx + (findgen(11)-5)
srcy2 = srcfit2[1].srcy + (findgen(11)-5)
flux1 = srcfit2[0].srcflux + (findgen(11)-5)/2
flux2 = srcfit2[1].srcflux + (findgen(11)-5)/2

c_fwhm1 = fltarr(11)
c_fwhm2 = fltarr(11)
c_srcx1 = fltarr(11)
c_srcy1 = fltarr(11)
c_srcx2 = fltarr(11)
c_srcy2 = fltarr(11)
c_flux1 = fltarr(11)
c_flux2 = fltarr(11)



srcin=srcfit2 ; outputs of previous run become inputs
rchisq=-1
;stop

set_plot,'ps'

dname=!d.name
if keyword_set(zbuffer) then set_plot,'Z' ; send plots to z buffer instead of screen or printer

hfwhm1=srcfit2[0].srcfwhm 
hfwhm2=srcfit2[1].srcfwhm 
hflux1=srcfit2[0].srcflux 
hflux2=srcfit2[1].srcflux 
hsrcx1=srcfit2[0].srcx    
hsrcy1=srcfit2[0].srcy    
hsrcx2=srcfit2[1].srcx    
hsrcy2=srcfit2[1].srcy    

if verbose then print,'COMPONENT 1 FWHM AND CHISQ'
for j=0,10 do begin
  srcin[0].srcfwhm=fwhm1[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[0].srcfwhm=', srcin[0].srcfwhm,' rchisq=',rchisq
  c_fwhm1(j)=rchisq
endfor

if verbose then print,'COMPONENT 1 FWHM AND CHISQ'

srcin=srcfit2
for j=0,10 do begin
  srcin[1].srcfwhm=fwhm2[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[1].srcfwhm=', srcin[1].srcfwhm,' rchisq=',rchisq
  c_fwhm2(j)=rchisq
endfor
if verbose then print,'vis.xyoffset=',visin(0).xyoffset(0),visin(0).xyoffset(1)
if verbose then print,'xyoffset=',xyoffset(0),xyoffset(1)
;stop


if verbose then print,'COMPONENT 1 FLUX AND CHISQ'

srcin=srcfit2
for j=0,10 do begin
  srcin[0].srcflux=flux1[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[0].srcflux=', srcin[0].srcflux,' rchisq=',rchisq
  c_flux1(j)=rchisq
endfor

if verbose then print,'COMPONENT 2 FLUX AND CHISQ'

srcin=srcfit2
for j=0,10 do begin
  srcin[1].srcflux=flux2[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[1].srcflux=', srcin[1].srcflux,' rchisq=',rchisq
  c_flux2(j)=rchisq
endfor


if verbose then print,'COMPONENT 1 X POSITION AND CHISQ'

srcin=srcfit2
for j=0,10 do begin
  srcin[0].srcx=srcx1[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[0].srcx=', srcin[0].srcx,' rchisq=',rchisq
  c_srcx1(j)=rchisq
endfor

if verbose then print,'COMPONENT 1 Y POSITION AND CHISQ'

srcin=srcfit2
for j=0,10 do begin
  srcin[0].srcy = srcy1[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[0].srcy=', srcin[0].srcy,' rchisq=',rchisq
  c_srcy1(j)=rchisq
endfor

if verbose then print,'COMPONENT 2 X POSITION AND CHISQ'

srcin=srcfit2
for j=0,10 do begin
  srcin[1].srcx = srcx2[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[1].srcx=', srcin[1].srcx,' rchisq=',rchisq
  c_srcx2(j)=rchisq
endfor

if verbose then print,'COMPONENT 2 Y POSITION AND CHISQ'

srcin=srcfit2
for j=0,10 do begin
  srcin[1].srcy=srcy2[j]
  hsi_vis_fwdfit_plotfit, visc, srcin, xyoffset,rchisq
  if verbose then print,' srcin[1].srcy=', srcin[1].srcy,' rchisq=',rchisq
  c_srcy2(j)=rchisq
endfor

time=visc[0].trange(0)
date=anytim(time(0),/ccsds,/date)
tim=str_replace(anytim(visc[0].trange[0],/ecs,/time,/trunc),':','')
set_plot,dname   ; restore plot device
if dname eq 'X' then window,/free,xsize=512,ysize=768
erange=visc[0].erange
estr=strcompress(fix(erange(0)),/rem)+'_'+strcompress(fix(erange(1)),/rem)

psname='chisq_slices'+date+'_'+tim+'E'+estr+'.ps
if !d.name eq 'PS' then device,/port,xsize=16,ysize=24,yoffset=1,file=psname

!p.multi=[0,2,4]
datim=anytim(time[0],/ecs)

!x.style=1 & !y.style=1
!p.charsize=1.5
yt='CHISQUARE'
xt='FWHM1 (arcsec)'
plot,fwhm1,c_fwhm1,psym=-6,tit='PROFILE FOR FWHM1',ytit=yt,xtit=xt
  legend,['hfwhm1='+strcompress(hfwhm1,/rem)],charsize=.75

xt='FWHM2 (arcsec)'
plot,fwhm2,c_fwhm2,psym=-6,tit='!7v!e2!n!3 PROFILE FOR FWHM2',ytit=yt,xtit=xt
  legend,['hfwhm2='+strcompress(hfwhm2,/rem)],charsize=.75
xyouts,0.025,0.5,datim,orient=90,align=0.5,/norm,charsize=1

xt='FLUX1 (photon cm!e-2!n s!e-1!n)'
plot,flux1,c_flux1,psym=-6,tit='!7v!e2!n!3 PROFILE FOR FLUX1',ytit=yt,xtit=xt
  legend,['hflux1='+strcompress(hflux1,/rem)],charsize=.75

xt='FLUX2 (photon cm!e-2!n s!e-1!n)'
plot,flux2,c_flux2,psym=-6,tit='!7v!e2!n!3 PROFILE FOR FLUX2',ytit=yt,xtit=xt
  legend,['hflux2='+strcompress(hflux2,/rem)],charsize=.75

xt='X1 (arcsec)'
plot,srcx1,c_srcx1,psym=-6,tit='!7v!e2!n!3 PROFILE FOR SRCX1',/XST,/YST,ytit=yt,xtit=x
  legend,['hsrcx1='+strcompress(hsrcx1,/rem)],charsize=.75

xt='Y1 (arcsec)'
plot,srcy1,c_srcy1,psym=-6,tit='!7v!e2!n!3 PROFILE FOR SRCY1',/XST,/YST,ytit=yt,xtit=xt
  legend,['hsrcy1='+strcompress(hsrcy1,/rem)],charsize=.75

xt='X2 (arcsec)'
plot,srcx2,c_srcx2,psym=-6,tit='!7v!e2!n!3 PROFILE FOR SRCX2',/XST,/YST,ytit=yt,xtit=xt
  legend,['hsrcx2='+strcompress(hsrcx2,/rem)],charsize=.75

xt='Y2 (arcsec)'
plot,srcy2,c_srcy2,psym=-6,tit='!7v!e2!n!3 PROFILE FOR SRCY2',/XST,/YST,ytit=yt,xtit=xt
  legend,['hsrcy2='+strcompress(hsrcy2,/rem)],charsize=.75


if !d.name eq 'PS' then begin
  device,/close
  if keyword_set(spawn_gs) then spawn,'gv '+psname
endif

; This is the same format as accepted by HSI_VIS_FWDFIT_ARRAY2STRUCTURE
;         [typ,flux1,srcx1,srcy1,fwhm1,loop_angl1,alb_ratio1,src_ht1,ecmsr1,eccen1]
;         [typ,flux2,srcx2,srcy2,fwhm2,loop_angl2,alb_ratio2,src_ht2,ecmsr2,eccen2]
osrcparm=[[1, 48.0, -798.0,-148.0, 9.0, 0, 0, 0, 0, 0],$
           [1, 36.0, -783.0,-128.0, 7.0, 0, 0, 0, 0, 0]]

close,1
fn='fit_data'+date+'_'+tim+'E'+estr+'.dat'

openw,1,fn
printf,1,anytim(visc[0].trange,/ecs),' Energy range = ',estr
;         01234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789
printf,1,'type    hflux     hsrcx     hsrcy     hfwhm loop_angl alb_ratio    src_ht     ecmsr     eccen'
printf,1,1.,   hflux1,hsrcx1, hsrcy1, hfwhm1,0.,0.,0.,0.,0., format='(i3,9f10.4)'
printf,1,1.,   hflux2,hsrcx2, hsrcy2, hfwhm2,0.,0.,0.,0.,0. ,format='(i3,9f10.4)'
printf,1,'type    eflux     esrcx     esrcy     efwhm loop_angl ealb_ratio   esrc_ht    e_ecmsr   e_eccen'
printf,1,'SD ',eflux1,esrcx1, esrcy1, efwhm1,0.,0.,0.,0.,0., format='(a3,9f10.4)'
printf,1,'SD ',eflux2,esrcx2, esrcy2, efwhm2,0.,0.,0.,0.,0. ,format='(a3,9f10.4)'

close,1

print,'type    eflux     esrcx     esrcy     efwhm loop_angl ealb_ratio   esrc_ht    e_ecmsr   e_eccen'
print,'SD ',eflux1,esrcx1, esrcy1, efwhm1,0.,0.,0.,0.,0., format='(a3,9f10.4)'
print,'SD ',eflux2,esrcx2, esrcy2, efwhm2,0.,0.,0.,0.,0. ,format='(a3,9f10.4)'

set_plot,dname

if !d.name ne 'PS' then begin
 tt=tvrd()
 write_png,'chisq_profiles'+datim+estr+'.png',255b-tt
 message,'saved '+'chisq_slices'+date+'.png',/info
endif


end

;savefiles=['vis3jun17_225240E30-40.sav',$
;'vis3jun17_225240E40-60.sav',$
;'vis3jun17_225310E30-40.sav',$
;'vis3jun17_225310E40-60.sav',$
;'vis3jun17_225340E30-40.sav',$
;'vis3jun17_225340E40-60.sav',$
;'vis3jun17_225410E30-40.sav',$
;'vis3jun17_225410E40-60.sav']

;savefiles=['vis2sep08_014220E25-40.sav','vis2sep08_014220E40-70.sav']

;savefiles=['vis3jul23_0041E40-80.sav','vis3jul23_0041E80-160.sav']
;savefiles=['vis3nov19_035928E25-50.sav']
;savefiles=['vis5aug25_043715E50-100.sav','vis5aug25_043716E50-100.sav']
;savefiles=['vis5jan15_223230E25-50.sav']

;  for j=0,7 do begin
;  for j=0,1 do begin
;for j=0,0 do begin
;  restore,savefiles(j),/verb
;  visfit,vis,srcparm=srcparm,event=8   ; jun 17, 2003
;  visfit,vis,srcparm=srcparm,event=9   ; sep08
;  visfit,vis,srcparm=srcparm,event=10  ; july 23, 2002
;  visfit,vis,srcparm=srcparm,event=11  ; nov 19, 2003
; visfit,vis,srcparm=srcparm,event=12   ; aug 25, 2005
; visfit,vis,srcparm=srcparm,event=13   ; jan 15, 2005
;endfor

;end