OpenMP support for gravmag.prism

Had to change the coordinate system shift (x[0] = x2 - xp[l]) to use
dx = x[i] - xp[l] inside the integration loops because I can't alter the
x array in parallel.

Benchmarks show 2x speedup in an i5 processor with 4 threads.
harvester doesn't seem to be much faster though.