I have finally made several tests and came to the following code corrections:
the fill fill_where_nan works ok, so I made no changes regarding it
added resolution = 1 for the slope/aspect estimation
changed standard resampling in sky_illumination to pyramid_scale = 2 and max_pyramid_radius = 20
added a few checks to see if the computed values are at least zero (if not they are in shade and should be zero)
I have finally made several tests and came to the following code corrections: