pro phase_vis,ut,vt,vist,mvist,svist,circles,dates,err_bars=err_bars, $
	phase=phase,help=help

if keyword_set(help) then begin
	print,'PURPOSE:'
	print,'	Plots RHESSI Visibilities and Model (MEM) Visibilities'
	print,'	vs. Theta for each UV circle used, for two flares'
	print,'	Keyword "phase" plots Phase and Model Phase vs. Theta'
	print,'HELP:'
	print,'	pro phase_vis,u,v,vis,mvis,svis,circles,dates,/err_bars,/phase'
	print,'	inputs: u,v (Fourier coordinates array)'
	print,'		vis (RHESSI Visibilities array,nflares x nvis)'
	print,'		mvis (Model Visibilities array, nflares x nvis)'
	print,'		svis (Sigma array, nflares x nvis)'
	print,'		circles (integer array, starting and ending UV circles used '
	print,'			in mem_map.pro, must be equal for both flares)'
	print,'		dates (string array of flare dates for title)'
	print,'		err_bars (default:no, /err_bars plots sigma error bars)'
	print,'		phase (default:no, /phase plots phase vs. theta)'
	print,'		*All inputs are individual flare arrays combined together*'
	print,'	outputs: PS file 2flares_vis.ps or 2flares_phase.ps'
	print,'AUTHOR:'
	print,'	Rick Pernak, Goddard Space Flight Center'
	return
endif

if keyword_set(phase) then fp='2flares_phase.ps' else fp='2flares_vis.ps'

set_plot,'ps'
device,filename=fp

nflares = size(vist,/n_dimensions)

multiplot,/def
for fctr=0,nflares-1 do begin ;loop for each flare
	if (nflares eq 1) then begin
		vis = vist
		svis = svist
		mvis = mvist
		u = ut
		v = vt
	endif else begin
		vis = vist(fctr,*)
		svis = svist(fctr,*)
		mvis = mvist(fctr,*)
		u = ut(fctr,*)
		v = vt(fctr,*)
	endelse
	
	vis_size = size(vis,/n_elements)
	mvis_size = size(mvis,/n_elements)
	if (vis_size ne mvis_size) then begin
		print,'Vis and Mvis sizes do not match, returning'
		return
	endif
	
	n = vis_size ;number of visibilities
	start_circle = circles[0]
	end_circle = circles[1]
	ncircles = (end_circle-start_circle)+1 ;number of circles
	vis_per = n/ncircles ;visibilities per UV circle
	
	for ctr=0,(ncircles-1) do begin ;loop over uv circles
		min = ctr*vis_per
		max = (ctr+1)*vis_per-1
		;reconstruct theta array so there's no wrap around
		;i.e. go from negative to positive
		theta = atan(v[min:max],u[min:max]) * !radeg;roll angle
		temp_circle = ctr+start_circle
		temp = sort(theta)
		theta = theta(temp)

		;arg: vis argument, marg: model vis argument, for plotting
		if keyword_set(phase) then begin
			;arguments for observed phase and model phase
			arg = atan(imaginary(vis[min:max]),float(vis[min:max]))*!radeg
			marg = atan(imaginary(mvis[min:max]),float(mvis[min:max]))*!radeg
			yr=[-180,180]
			
			;shift phases 360 degrees so observed and model have same sign
			for nvis=0,vis_per-1 do begin
				if (marg[nvis] - arg[nvis] gt 180) then arg[nvis]=arg[nvis]+360.
				if (marg[nvis] - arg[nvis] lt -180) then arg[nvis]=arg[nvis]-360.
			endfor

			print,mean(marg-arg)

		endif else begin
			;arguments for observed and model amplitudes
			arg = abs(vis[min:max])
			marg = abs(mvis[min:max])
		endelse
		
		;assign left side y title
		if (fctr eq 0) then uv = string(temp_circle,format='("UV Circle",1x,i1)') $
			else uv = ''

		;assign x titles for both flares
		if (ctr eq ncircles-1) then xt='Roll Angle (Degrees)' else xt=''
		
		;assign main titles for flares
		if (ctr eq 0) then title=dates[fctr] else title=''
		
		;make the friggin plots
		if (fctr eq 0) and (ctr eq 0) then multiplot,[0,nflares,ncircles,0,1]
		if keyword_set(phase) then begin
			
			;phase plots
			plot,theta,arg,title=title,ytitle=uv,xtitle=xt,/xst,charsize=.75,psym=6,$
				yr=yr,ytickformat='(A1)',xtickformat='(A1)',/yst
			oplot,theta,marg,psym=5
			legend,['RHESSI','Model'],psym=[6,5],box=0,charsize=.75,spacing=.75
			
			;right hand ytitle
			if (fctr eq 1) and (ctr eq 2) then begin
				axis,yaxis=1,ytitle='Phase (Degrees)',ytickformat='(A1)',xtickformat='(A1)'
			endif
			
			if keyword_set(err_bars) then begin
				; from Ed, "recipe" for error bars in phases
				x = float(vis[min:max])
				y = imaginary(vis[min:max])
				amp = sqrt(x^2 + y^2)
				sarg = svis[min:max]
				err1 = atan(y+sarg*(x/amp),x+sarg*(-y/amp))*!radeg
				err2 = atan(y+sarg*(-x/amp),x+sarg*(y/amp))*!radeg
				low_err = dblarr(vis_per)
				hi_err = dblarr(vis_per)
				for j=0,vis_per-1 do begin ;low and high limits for error bars
					low_err[j] = min([err1[j],err2[j]])
					hi_err[j] = max([err1[j],err2[j]])
					;change error bars so they correspond with 360 degree phase shift
					if (hi_err[j] - low_err[j] gt 180) and (marg[j] gt 0) then begin
						temp = hi_err[j]
						hi_err[j] = low_err[j]+360.
						low_err[j] = temp
					endif else if (hi_err[j]-low_err[j] gt 180) and (marg[j] lt 0) $
						then begin
						temp = low_err[j]
						low_err[j] = hi_err[j]-360.
						hi_err[j] = temp
					endif else if $
						((marg[j] gt 0) and (arg[j] gt 0) and (hi_err[j] lt 0) and (low_err[j] lt 0)) $
						or ((marg[j] lt 0) and (arg[j] lt 0) and (hi_err[j] gt 0) and (low_err[j] gt 0)) $
						then begin
						hi_err[j] = -hi_err[j]
						low_err[j] = -low_err[j]
					endif
				endfor
			errplot,theta,low_err,hi_err
			
			;calculate sigma (svis) for phases and phase chi^2 for each circle
			svis_phase = hi_err + abs(low_err)
			chi2,arg,marg,svis_phase,circles,chi2,/phase
			chi2 = string(chi2[1],format='(d7.2)')
			legend,[chi2],/right,/bottom,box=0,charsize=.75
			endif
			multiplot
		
		;amplitude plots
		endif else begin
			;calculate amplitude chi^2 for each uv circle
			chi2,arg,marg,svis[min:max],circles,chi2
			chi2=string(chi2[1],format='(d7.2)')
			plot,theta,arg,title=title,ytitle=uv,xtitle=xt,/xst,charsize=.75,$
				ytickformat='(A1)',xtickformat='(A1)'
			oplot,theta,marg,linestyle=2
			legend,['RHESSI','Model'],linestyle=[0,2],box=0,charsize=.75,spacing=.75
			legend,[chi2],/right,/bottom,box=0,charsize=.75

			;right hand ytitle
			if (fctr eq 1) and (ctr eq 2) then begin
				axis,yaxis=1,ytitle='Visibility Amplitude',ytickformat='(A1)',xtickformat='(A1)'
			endif
			
			if keyword_set(err_bars) then begin
				oploterr,theta,arg,svis[min:max]
				oplot,theta,marg,linestyle=2
			endif
			multiplot
		endelse
	endfor
endfor

multiplot,/reset
multiplot,/def
device,/close
set_plot,'x'

return
end