function mom=momlambert(dig,N,a,eta)
digits(dig);
for k=1:2*N
  g=feval(symengine,'igamma',a+1,(k-a)*log(vpa(1/eta)));
  igam=vpa(g);
  mom(k)=igam/vpa((k-vpa(a))^vpa(a+1));
end
