delta := 2.5971307458876746E-6; deltav := 0.004772707592564195; A := 102.62012490987308; # Conditions pi_r := 348.0; delta_r := 68.0; i := 1.0; k0 := 5447; deltap := deltav+delta; # Following are from http://terrytao.wordpress.com varpi := (1-delta_r*delta)/pi_r; theta := deltap / (1/4 + varpi); thetat := min( ((i*deltap - delta)/2.0 + varpi) / (1/4 + varpi), 1); 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); alpha := j^2 / (4 * (k0-1)); e := exp( A + (k0-1) * int( exp(-(A+2*alpha)*t)/t, t=deltat..theta, numeric ) ); gd := (j^2/2) * BesselJ(k0-3,j)^2; 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";