(old) Bugfix in Igeom=1 Lcons=3 and output of stability variables

[ErolBa]

This pull request involves the following modifications:

1. Primary modification: fixes a bug when running with Igeometry=1 and Lconstraint=3, adds the total pflux as a conserved quantity (otherwise the Beltrami solver is missing a conserved quantity).
2. Outputs .hessian matrix for stability analysis
3. Calculates and outputs plasma beta
4. Outputs_ quantities related to stability in the .h5 file (eigvals, eigvecs, and the hessian matrix)

[ErolBa]

The tests pass (if trivially Lconstraint and LBeltrami are not checked, we can change those). However, the test for slab (G=1) seems lacking, and there should be a test involving different Lconstraints. Particularly, with Lcons=3 where the bug was present.