oldtests/gsl-intgamma/quad.mathematica

   1 SigmaIntegrand[n_, k_, R_, x_] := Exp[-(R x)^2 + (k/(2 x))^2] x^(2 n) ;
   2 l:= ReadList["values.in", Table[Number, {4}]];
   3 For[i = 1, i <= Length[l], i++,
   4   p := l[[i]];
   5   n := p[[1]];
   6   k := p[[2]];
   7   r := p[[3]];
   8   eta := p[[4]];
   9   res := NIntegrate[SigmaIntegrand[n, k, r, x], {x,eta,Infinity}];
  10   Print[InputForm[n], " ", CForm[k], " ", CForm[r], " ", CForm[res]]
  11 ];
  12