e = Sqr(2 * f - f ^ 2)
C = a / Sqr(1 - e ^ 2)
e2 = e / Sqr(1 - e ^ 2)
beita0 = 1 - (3 / 4) * e2 ^ 2 + (45 / 64) * e2 ^ 4 - (175 / 256) * e2 ^ 6 +
(11025 / 16384) * e2 ^ 8
beita2 = beita0 - 1
beita4 = (15 / 32) * e2 ^ 4 - (175 / 384) * e2 ^ 6 + (3675 / 8192) * e2 ^ 8
beita6 = (-35 / 96) * e2 ^ 6 + (735 / 2048) * e2 ^ 8
beita8 = (315 / 1024) * e2 ^ 8
For i = 1 To k
B0 = x(i) / (C * beita0)
a1 = a * Cos(B0) / Sqr(1 - e ^ 2 * (Sin(B0)) ^ 2)
l0 = y(i) / a1
Do
fxb = 0
fxbl = 0
fybl = 0
t = Tan(B0)
yita = e2 * Cos(B0)
n = a / Sqr(1 - e ^ 2 * (Sin(B0)) ^ 2)
a2 = (1 / 2) * n * Sin(B0) * Cos(B0)
a3 = (1 / 6) * n * (Cos(B0)) ^ 3 * (1 - t ^ 2 + yita ^ 2)
a4 = (1 / 24) * n * Sin(B0) * (Cos(B0)) ^ 3 * (5 - t ^ 2 + 9 * yita ^ 2 + 4 *
yita ^ 4)
a5 = (1 / 120) * n * (Cos(B0)) ^ 5 * (5 - 18 * t ^ 2 + t ^ 4 + 14 * yita ^ 2 -
58 * yita ^ 2 * t ^ 2)
a6 = (1 / 720) * n * Sin(B0) * (Cos(B0)) ^ 5 * (61 - 58 * t ^ 2 + t ^ 4)
fxb = fxb + (C * beita6 + C * beita8 * (Cos(B0)) ^ 2) * (Cos(B0)) ^ 2
fxb = (fxb + C * beita4) * (Cos(B0)) ^ 2