"""
Reverse transposition limitations
====================================

Unfortunately, sometimes there is not a unique solution.

Author: Anton Driesse

"""

# %%
#
# Introduction
# ------------
# When irradiance is measured on a tilted plane, it is useful to be able to
# estimate the GHI that produces the POA irradiance.
# The estimation requires inverting a GHI-to-POA irradiance model,
# which involves two parts:
# a decomposition of GHI into direct and diffuse components,
# and a transposition model that calculates the direct and diffuse irradiance
# on the tilted plane.
# Recovering GHI from POA irradiance is termed "reverse transposition."
#
# Unfortunately, for a given POA irradiance value, sometimes there is not a
# unique solution for GHI.
# Different GHI values can produce different combinations of direct and
# diffuse irradiance that sum to the same POA irradiance value.
#
# In this example we look at a single point in time and consider a full range
# of possible GHI and POA global values as shown in figures 3 and 4 of [1]_.
# Then we use :py:meth:`pvlib.irradiance.ghi_from_poa_driesse_2023` to estimate
# the original GHI from POA global.
#
# References
# ----------
# .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
#     Perez diffuse sky model for forward and reverse transposition.
#     Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
#

import numpy as np

import matplotlib
import matplotlib.pyplot as plt

from pvlib.irradiance import (erbs_driesse,
                              get_total_irradiance,
                              ghi_from_poa_driesse_2023,
                              )

matplotlib.rcParams['axes.grid'] = True

# %%
#
# Define the conditions that were used for figure 3 in [1]_.
#

dni_extra = 1366.1
albedo = 0.25
surface_tilt = 40
surface_azimuth = 180

solar_azimuth = 82
solar_zenith = 75

# %%
#
# Define a range of possible GHI values and calculate the corresponding
# POA global.  First estimate DNI and DHI using the Erbs-Driesse model, then
# transpose using the Perez-Driesse model.
#

ghi = np.linspace(0, 500, 100+1)

erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra)

dni = erbsout['dni']
dhi = erbsout['dhi']

irrads = get_total_irradiance(surface_tilt, surface_azimuth,
                              solar_zenith, solar_azimuth,
                              dni, ghi, dhi,
                              dni_extra,
                              model='perez-driesse')

poa_global = irrads['poa_global']

# %%
#
# Suppose you measure that POA global is 200 W/m2. What would GHI be?
#

poa_test = 200

ghi_hat = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth,
                                    solar_zenith, solar_azimuth,
                                    poa_test,
                                    dni_extra,
                                    full_output=False)

print('Estimated GHI: %.2f W/m².' % ghi_hat)

# %%
#
# Show this result on the graph of all possible combinations of GHI and POA.
#

plt.figure()
plt.plot(ghi, poa_global, 'k-')
plt.axvline(ghi_hat, color='g', lw=1)
plt.axhline(poa_test, color='g', lw=1)
plt.plot(ghi_hat, poa_test, 'gs')
plt.annotate('GHI=%.2f' % (ghi_hat),
             xy=(ghi_hat-2, 200+2),
             xytext=(ghi_hat-20, 200+20),
             ha='right',
             arrowprops={'arrowstyle': 'simple'})
plt.xlim(0, 500)
plt.ylim(0, 250)
plt.xlabel('GHI [W/m²]')
plt.ylabel('POA [W/m²]')
plt.show()

# %%
#
# Now change the solar azimuth to match the conditions for figure 4 in [1]_.
#

solar_azimuth = 76

# %%
#
# Again, estimate DNI and DHI using the Erbs-Driesse model, then
# transpose using the Perez-Driesse model.
#

erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra)

dni = erbsout['dni']
dhi = erbsout['dhi']

irrads = get_total_irradiance(surface_tilt, surface_azimuth,
                              solar_zenith, solar_azimuth,
                              dni, ghi, dhi,
                              dni_extra,
                              model='perez-driesse')

poa_global = irrads['poa_global']

# %%
#
# Now reverse transpose all the POA values and observe that the original
# GHI cannot always be found.  There is a range of POA values that
# maps to three possible GHI values, and there is not enough information
# to choose one of them.  Sometimes we get lucky and the right one comes
# out, other times not.
#

result = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth,
                                   solar_zenith, solar_azimuth,
                                   poa_global,
                                   dni_extra,
                                   full_output=True,
                                   )

ghi_hat, conv, niter = result
correct = np.isclose(ghi, ghi_hat, atol=0.01)

plt.figure()
plt.plot(np.where(correct, ghi, np.nan), np.where(correct, poa_global, np.nan),
         'g.', label='correct GHI found')
plt.plot(ghi[~correct], poa_global[~correct], 'r.', label='unreachable GHI')
plt.plot(ghi[~conv], poa_global[~conv], 'm.', label='out of range (kt > 1.25)')
plt.axhspan(88, 103, color='y', alpha=0.25, label='problem region')

plt.xlim(0, 500)
plt.ylim(0, 250)
plt.xlabel('GHI [W/m²]')
plt.ylabel('POA [W/m²]')
plt.legend()
plt.show()