Now reads output wordsize from header; supports double or single.
Switched to a less efficient but more robust root finder.

When there are multiple roots, this will find the first or
the last root, depending on the threshold.  This is important
when the vertical density profile, f_xy, is not monotonic.
Added calculation of <f (1-f)>_{x,y} and the \Theta mixing function.

The array <f (1-f)>_xy is saved to the npz file and \Theta is saved
to the .dat file.

Also, a bit of cleaning.