PRO axaf1_observ, resolving

@effic_com
@axaf1_com

@b_on_w
if(n_elements(resolving) EQ 0) then resolving = 900.
if(n_elements(ax1_file) EQ 0) then ax1_file = "s_Capella_hires.ax1"
if(n_elements(ax_time) EQ 0) then ax_time = 10.E3

; Read the input file:
all_data = fltarr(1000,4)
OPENR, 37, ax1_file
Ndata = 0
tmp0 = 0.
tmp1 = 0.
tmp2 = 0.
tmp3 = 0.
WHILE (not EOF(37)) DO BEGIN
  READF, 37,FORMAT='(4F15.8)', tmp0, tmp1, tmp2, tmp3
  all_data(Ndata,0) = tmp0
  all_data(Ndata,1) = tmp1
  all_data(Ndata,2) = tmp2
  all_data(Ndata,3) = tmp3
  Ndata = Ndata+1
  ENDWHILE
CLOSE, 37

; Fill the arrays:
es = fltarr(Ndata)
spec_heg = es
spec_meg = es
spec_total = es
es = all_data(0:Ndata-1,0)
spec_heg = all_data(0:Ndata-1,1)
spec_meg = all_data(0:Ndata-1,2)
spec_total = all_data(0:Ndata-1,3)


; Convolve with a Gaussian
; resolving = 900.
g_width = es(0)/(resolving*2.35)
dx = es(1)-es(0)
; Generate the Gaussian kernal: (to +/- 3 sigma)
n_kernal = 2*FIX(3.5*g_width/dx) + 1
print, ' Kernal size = ', n_kernal
kernal = fltarr(n_kernal)
; Fill the kernal
for ik=0,n_kernal-1 do begin
  mid_val = FLOAT(ik)-FLOAT((n_kernal-1)/2)
  norm_min = dx*(mid_val-0.5)/g_width
  norm_max = dx*(mid_val+0.5)/g_width
;  The IDL errorf is integral from 0 to x of gauss with sigma = 0.707
;   and normalization (-inf to + inf) = 2.0 ... arg!
  kernal(ik) = 0.5*(errorf(norm_max*0.707)-errorf(norm_min*0.707))
;  print, norm_min, norm_max, kernal(ik)
end
print, '  Kernal:'
print, kernal
print, ' '
ax_kernal = kernal

blur_counts = CONVOL(spec_total,ax_kernal)*ax_time

; Add Poisson noise: convert counts to observed counts
blur_counts = FLOAT(g_poiss(blur_counts))

; Output to a file:
; and output to a file = ax1_file.o
IF(1 EQ 1) then begin
ax_outfile = ax1_file + ".o"
OPENW, 37, ax_outfile
FOR id = 0,n_elements(blur_counts)-1 DO BEGIN
  PRINTF, FORMAT='(2F15.8)', 37, es(id), blur_counts(id)
END
CLOSE, 37
END ; end of if(1 EQ 1)

; Calculate some statistics
t_counts = TOTAL(blur_counts)

print, " Total counts in "+STRING(ax_time)+" seconds: "+STRING(t_counts)

;max = resolving
max = 1000.
if(max LT 100.) then max = 100.
!MTITLE = "Simulated HETG Observation of Capella"
plot_yt = "Counts in bin (bin width is  E/"+$
	STRING(es(0)/(es(1)-es(0)),FORMAT='(F5.0)')+")"
plot, es, blur_counts, PSYM=10, $
	XTITLE = "Energy (keV)", YTITLE = plot_yt, $
	XRANGE = [0.9, 1.5], YRANGE = [-0.1*max,max]
t_string = " Exposure: "+STRING(ax_time,FORMAT='(F7.0)')+" seconds"
c_string = " Total counts:"+STRING(t_counts,FORMAT='(I6)')+$
	STRING(FORMAT='("  (",F4.1," -",F4.1,"keV)")',MIN(es),MAX(es))
r_string = " Resolving power is: "+STRING(resolving,FORMAT='(F5.0)')
x_loc = 0.60
xyouts, x_loc, 0.90, t_string, /NORMAL
xyouts, x_loc, 0.87, c_string, /NORMAL
xyouts, x_loc, 0.84, r_string, /NORMAL

oplot, [0.8, 1.6], [0.0,0.0], LINESTYLE=1
; add atro lines!
ral = max/1000.
astro_lines
al_all_lines, -45.*ral, -10.0*ral
for il= 49, 55 do al_line_name, il, -155.*ral
al_line_name, 57, -155.*ral
for il= 59, 65 do al_line_name, il, -155.*ral
al_line_name, 67, -155.*ral
al_line_name, 69, -155.*ral
RETURN
END