Fluxes of coincident stars are not superimposing properly

I noticed that
test double stars.pdf
in the case multiple stars are projected on the same pixel of detector the resulting image is not based on the sum of the different fluxes. Instead, the information of the (last?) star erase the previous information.

Attached an example where in a corner I put (on a pixel) 4 stars of mag 10 and on the opposite a star of mag ~8.49485 (equivalent to the sum of the four fluxes of mag10).

Probably, making a sum instead than a substitution in the making of the object image will solve the problem.
This is a serious limitation in the case many stars (>> thousands) are used on random positions.

I checked on the generated fits file in the example that the peaks are respectively

3.12 10**6
1.24 10**7
Exactly a factor four

Thanks @carmeloarci !
These lines:

ScopeSim/scopesim/optics/fov.py

Lines 397 to 398 in d762bcf

    
           canvas_image_hdu.data[y[mask], x[mask] 
        
                                 ] += flux[mask] * weight[mask]

Should be replaced with something like

                for yi, xi, fluxi, weighti in zip(
                    y[mask], x[mask], flux[mask], weight[mask]
                ):
                    canvas_image_hdu.data[yi, xi] += fluxi * weighti

where x and y are the pixel positions of the 5 sources, and flux and weight correspond to the 5 magnitude values.

If there is overlap in the x and y values, then the += will only use the last value.

See this 1d example:

>>> import numpy as np
>>> data = np.zeros(8)
>>> x = [1, 1, 1, 1, 5]
>>> flux = [2, 2, 2, 2, 8]
>>> data[x] += flux
>>> data
array([0., 2., 0., 0., 0., 8., 0., 0.])

We should indeed fix this, and also write a test for this to prevent regressions.

Stacking stars should work properly now on dev_master, thanks for the issue and the easily reproducible example.

Thanks! I also tested the new version from github. Test passed

	canvas_image_hdu.data[y[mask], x[mask]
	] += flux[mask] * weight[mask]