;----define the gamma distribution
;function gamma_fit, shape, slope, N0, r
;  return, N0*r^shape*exp(-slope*r)/gamma(shape+1)
;end
;
pro gamma_fit, r, A, f
  shape=A[0]
  slope=A[1]
  N0=A[2]
  f=N0*r^shape*exp(-slope*r)/gamma(shape+1)
end