Convenience function for combining arbitrary time-length exposures

I had at one time combined various length exposures via a scaling factor inversely proportional to exposure time, using masked arrays where each contributing exposure length had clip/mask points (primarily white clip in long exposures and dark in short). Code is long gone... This was simply via an image pre-scale before ma summing. Other scalings are probably useful as well.

The result is high dynamic range, with variable noise in 2D space. I'm sure a 2D STD could be calculated per pixel, but I did not. It would seem that the variance would scale with the multiplier and by n**0.5 for n non-masked contributions.

This is computationally similar to these topics https://github.com/astropy/ccdproc/issues/301 and https://github.com/astropy/ccdproc/issues/360

However since image registration is desirable as an interim step it may be better as a recipe than a method on an image cube.