HP42s first major program (Double Integral) Best way to approach?

[Albert Chan]

Borrowing my AGM2 method of calculating ellipse perimeter, but factor out AGM.

Code:

local sqrt, pi = math.sqrt, math.pi

function agm(a, b)
    local t = 0.5
    local s = a*a
    while true do
        local k = b-a
        local s1 = s - t*k*k
        if s == s1 then return a+0.5*k, s end
        b = sqrt(a*b)
        a = a + 0.5*k
        t = t + t
        s = s1
    end
end

function perimeter(a, b)    -- ellipse perimeter
    a, b = agm(a+b, 2*sqrt(a*b))
    return pi*b/a           -- b/a = "diameter"
end

We can get *both* K(m), E(m) from 1 call of agm(), thus bore-hole removed volume.
see equation 12: https://mathworld.wolfram.com/SteinmetzSolid.html

Code:

function KE(m)    -- return K(m), E(m), m < 1
    m = sqrt(1-m*m)
    local K, E = agm(1+m, 2*sqrt(m))
    K = pi/K
    return K, K*E/4
end

function hv(d,D)  -- bore-hole removed volume, thru solid pipe
    local K, E = KE(d/D)
    return D/3 * (d*d*(K+E) - D*D*(K-E))
end

lua> a, b, c = 24*2, 29*2, 12*2    -- all diameters
lua> hv(c,b) - hv(c,a)                   -- assumed penetrating pipe both side
4651.24643661709