pro gamow,t_plasma,z1=z1,z2=z2,a1=a1,a2=a2

; displays gamow peak for input value of t (in K)  
; default values for protons
  if not keyword_set(z1) then z1=1. 
  if not keyword_set(z2) then z2=1.
  if not keyword_set(a1) then a1=1.
  if not keyword_set(a2) then a2=1.

  kt=findgen(200)/10.
  kt1=kt(1)
  kt=kt+kt1 ; shift right to avoid zero at first energy point
  
  kb=1.3806e-16
  kt_plasma=kb*t_plasma/1.602e-9 ; in kev

  ekin=kt*kt_plasma
  p1=exp(-ekin/kt_plasma)
  
  b=31.28*z1*z2*sqrt(a1*a2/(a1+a2)) ; Gamow energy in keV
  p2=exp(-b/sqrt(ekin))
  ptot=p1*p2
  
  window,0
  ymin=max(ptot)/100.
  plot_io,kt,p1,yrange=[ymin,1.],line=1, $
	   title='Gamow peak',xtitle='Ekin/kT_plasma', $
	   ytitle='probability',charsize=1.5; logarithmic display
  oplot,kt,p2,line=2
  oplot,kt,ptot
  window,1
  plot,kt,ptot,title='Gamow peak',xtitle='Ekin/kT_plasma', $
	ytitle='probability', charsize=1.5 ; linear display
  return
end