[FIX] iam.physical returns nan for aoi > 90° when n = 1 (#1706)

[kdebrab]

• Closes regression: iam.physical returns nan for aoi > 90° when n = 1 #1706
• I am familiar with the contributing guidelines
• Tests added
• Updates entries in docs/sphinx/source/reference for API changes.
• Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• Pull request is nearly complete and ready for detailed review.
• Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

[cwhanse]

The nan appears to enter here and on the following line, because n1costheta1 and n2costheta1 are both zero when AOI>90, n=1 and n_ar is None. I wonder if there is a condensed mathematical expression for the IAM in the limit case of n1=1, n2=1, and n3=1.

[kdebrab]

The nan appears to enter here and on the following line, because n1costheta1 and n2costheta1 are both zero when AOI>90, n=1 and n_ar is None. I wonder if there is a condensed mathematical expression for the IAM in the limit case of n1=1, n2=1, and n3=1.

That condensed mathematical expression is np.exp(-K * L / costheta) / np.exp(-K * L) (I think) but it has the same problem of a division by zero for AOI>90. Moreover, we also need to support n2 as a vector (see #1616 (comment)), where potentially some elements are 1 and others not. For both cases we need the lines included in the current pull request.
In the end, I'm not sure whether it's worth to have a special path for a minor speedup of a rather uncommon use case.

[mikofski]

Thanks

[kandersolar]

Thanks @kdebrab! Can you add a bug fix entry in the 0.9.6 what's new file?

One thing is that the new test produces some runtime warning clutter: https://github.com/pvlib/pvlib-python/actions/runs/4510438842/jobs/7941391259?pr=1707#step:9:58

It would be nice to prevent iam.physical from emitting those, but n=1 seems uncommon enough that I'd be okay with just suppressing them in the test function if there's not a clean way of preventing them at the source.

[adriesse]

There is an interesting thing happening with this sequence of inexact trig calculations:

costheta = np.maximum(0, cosd(aoi))  # --> 6.123233995736766e-17
sintheta = np.sqrt(1 - costheta**2)  # --> 1.0
costheta = np.sqrt(1 - sintheta**2) # --> 0.0

If the first line were producing a true zero this function would blow up sooner and more often; therefore, the previous use of zeroang = 1e-06 might be the preferred solution.

We shouldn't really rely on cosd(90) != 0 anywhere, although we probably unwittingly do.

[adriesse]

I would follow this up with a suggestion (now retracted) to rewrite tools.cosd as

np.sqrt(1 - np.square(np.sin(angle * (np.pi / 180))))

[cwhanse]

I would follow this up with a suggestion to rewrite tools.cosd as

np.sqrt(1 - np.square(np.sin(90 * (np.pi / 180))))

I'm puzzled by this. I assume you meant theta * (np.pi/180)). How would it properly evaluate cosd(135)? It can't return a negative value, since np.sqrt only returns the principal root.

[adriesse]

I'm puzzled by this

Sorry, didn't think that through all the way.

Still, I hope this does not detract from the preceding (hopefully) valid observation about what's going wrong in this function.

[kandersolar]

I'm going to merge this since it's ready to go and an improvement. Additional improvements can always be future PRs :)

Thanks again @kdebrab for noticing and fixing this bug!