varpi_v := 490.9173761236293; deltap := 0.00987167117510929; A := 965.0207482345567; pi_r := 50; delta_r := 5; k0 := 272; # Following are from http://terrytao.wordpress.com varpi := 1/pi_r - varpi_v / 10^5; delta := (1 - pi_r * varpi) / delta_r; theta := deltap / (1/4 + varpi); # Gergely's improved value for thetat thetat := ((deltap - delta)/2 + varpi) / (1/4 + varpi); deltat := delta / (1/4 + varpi); j := BesselJZeros(k0-2,1); eps := 1 - j^2 / (k0 * (k0-1) * (1+4*varpi)); kappa1 := int( (1-t)^((k0-1)/2)/t, t = theta..1, numeric); kappa2 := (k0-1) * int( (1-t)^(k0-1)/t, t=theta..1, numeric); # using Gergely and Eytan's improved kappa_3 alpha := j^2 / (4 * (k0-1)); e := exp( A + (k0-1) * int( exp(-(A+2*alpha)*t)/t, t=deltat..theta, numeric ) ); # using Gergely's exact expression for denominator gd := (j^2/2) * BesselJ(k0-3,j)^2; # using Eytan's exact expression for numerator tn := sqrt(thetat)*j; gn := (tn^2/2) * (BesselJ(k0-2,tn)^2 - BesselJ(k0-3,tn)*BesselJ(k0-1,tn)); kappa3 := (gn/gd) * e; eps2 := 2*(kappa1+kappa2+kappa3); # we win if eps2 < eps epsdiff := evalf(eps2-eps); save epsdiff, "output.txt";