The C routine dvolume

The C routine dvolume first checks the boundary conditions. Then, depending on the value of $\Omega_K$ the proper motion distance (eq.  4) is calculated. On doing this, the same problem with the IDL environment as described in VIE occurs. It is solved in the same way as described there.
Actually the routine dvolume calculates three values for d_M which are called dm, dm_p, dm_m, p and m standing for plus and minus. dm_p is calculated at z+dz and dm_m at z-dz where dz is 10^-4. These values are then used for the approximation of d(d_m)/dz. The first term on the right handside of eq.  9 is calculated straightforward. For the second term of eq.  9, the following approximation is used:

$$\frac{d(d_m)}{dz} \simeq \frac{dm\_p - dm\_m}{dz}$$ (12)

which is a pretty rough approximation but it is good if dz is small enough (in calc_lib.so 10^-4 is used which is sufficient).
After eq.  9 has been calculated as described above, the result is returned.