"""Tutorial 23 - Irradiance map over flat ground (analytic check) and over a DEM."""
import math
import numpy as np
from scipy.ndimage import gaussian_filter
import less

horizontal_sun = np.array([1.6203, 0.7190])     # expected horizontal irradiance, W m-2 nm-1
flat = less.Scene()
flat.size = 40.0
flat.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
flat.illumination = less.Illumination(
    source=less.Sun(zenith=30, azimuth=180, wavelengths=[550, 850],
                    irradiance=horizontal_sun / math.cos(math.radians(30))),
    atmosphere=less.NoAtmosphere())
check = flat.simulate(less.IrradianceMapSensor(bands=[550, 850], resolution=1.0, photons_per_m2=400))
print("down:", check.downwelling_irradiance().mean(axis=(0, 1)), "expected", horizontal_sun)
print("up:  ", check.upwelling_irradiance().mean(axis=(0, 1)), "expected", 0.15 * horizontal_sun)

rng = np.random.default_rng(42)
dem = np.zeros((128, 128))
for sigma, amplitude in [(18, 6), (9, 3), (4, 1.5), (2, 0.5)]:
    dem += amplitude * gaussian_filter(rng.standard_normal((128, 128)), sigma=sigma)
dem = (dem - dem.min()) * 12.0 / (dem.max() - dem.min())
hills = less.Scene()
hills.size = 100.0
hills.terrain = less.Terrain(property=less.Lambertian(reflectance=0.30), dem=dem.astype(np.float32))
hills.illumination = less.Illumination(source=less.Sun(zenith=45, azimuth=180), atmosphere=less.NoAtmosphere())
result = hills.simulate(less.IrradianceMapSensor(bands=[550], resolution=0.78, photons_per_m2=2000))
np.save("23_downwelling_horizontal.npy", result.downwelling_irradiance("horizontal")[:, :, 0])
np.save("23_downwelling_slope.npy", result.downwelling_irradiance("slope")[:, :, 0])