pro gser,gamser,a,x,gln
  ITMAX=100
  EPS=3.0e-7
  ;Returns the incomplete gamma function P(a, x) evaluated by its series representation as gamser.
  ;Also returns ln Gamma(a) as gln.
  ;int n;
  ;float sum,del,ap;
  gln=alog(gamma(a));
  if (x <= 0.0) then begin
    if (x < 0.0) then message,"x less than 0 in routine gser";
    gamser=0.0;
    return;
  end else begin
    ap=a;
    sum=1.0/a;
    del=sum;
    for n=1,ITMAX do begin
      ap++;
      del *= x/ap;
      sum += del;
      if (abs(del) lt abs(sum)*EPS) then begin
        gamser=sum*exp(-x+a*alog(x)-(gln));
        return;
      end
    end
    message,nrerror("a too large, ITMAX too small in routine gser");
    return;
  end
end

pro gcf,gammcf,a,x,gln
  ITMAX=100; Maximum allowed number of iterations.
  EPS=3.0e-7; Relative accuracy.
  FPMIN=1.0e-30; Number near the smallest representable floating-point number.
  ;Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation
  ;as gammcf. Also returns ln Gamma(a) as gln.
  ;int i;
  ;float an,b,c,d,del,h;
  gln=alog(gamma(a));
  b=x+1.0-a; Set up for evaluating continued fraction
           ; by modified Lentz’s method (§5.2)
           ; with b0 = 0.
  c=1.0/FPMIN;
  d=1.0/b;
  h=d;
  i=1
  del=10
  while i le ITMAX and abs(del-1) ge EPS do begin
    an = -i*(i-a);
    b += 2.0;
    d=an*d+b;
    if (abs(d) lt FPMIN) then d=FPMIN;
    c=b+an/c;
    if (abs(c) lt FPMIN) then c=FPMIN;
    d=1.0/d;
    del=d*c;
    h *= del;
  end
  if (i > ITMAX) then message,"a too large, ITMAX too small in gcf";
  gammcf=exp(-x+a*alog(x)-(gln))*h; Put factors in front.
end
function gammq,a,x
; Returns the incomplete gamma function Q(a, x) = 1 - P(a, x).
;void gcf(float *gammcf, float a, float x, float *gln);
;void gser(float *gamser, float a, float x, float *gln);
  if (x lt 0.0 or a le 0.0) then message, "Invalid arguments in routine gammq";
  if (x lt (a+1.0)) then begin ;Use the series representation
    gser,gamser,a,x,gln;
    return,1.0-gamser; and take its complement.
  end else begin
    ; Use the continued fraction representation.
    gcf,gammcf,a,x,gln;
    return,gammcf;
  end
end