Avoid division by zero error if integral flux is 0 in eval_gradients().

This fixes issue #1141.
Add some debug outputs.
Update after bug fix
Add information about energy dispersion implementation.