function bn_vs_n, x

  common slope, a, c

  if x gt 0 then begin
    ga = gamma(a)
    dif = ga - 2. * igamma(a, x) * ga
    c = dif
    return, dif
  endif else begin
    return, c*1.2
  endelse

end




function sersic_bn, nmin, nmax, dn

  common slope, a, c

  c = 0d
  nrange = nmax - nmin
  nn = nrange / dn
  n = findgen(nn) / float(nn) * nrange + nmin + dn

  bn = dblarr(nn)
  for i = 0, nn-1 do begin
    a = 2d * n[i]
    if i eq 0 then guess=1d else guess = bn[i-1]
    bn[i] = newton(guess, 'bn_vs_n', /double)
  endfor
  n2 = 2d * n
  kn = 2d * !dpi * n * exp(bn) / bn^n2 * gamma(n2)

  return, {n:n, bn:bn, kn:kn}

  

end