Add marion_diffuse_tracking, to be used instead of marion_diffuse for vector surface_tilt inputs

[cbcrespo]

• Closes pvlib.iam.marion_integrate uses too much memory for vector inputs #1402
• I am familiar with the contributing guidelines
• I attest that all AI-generated material has been vetted for accuracy and is in compliance with the pvlib license
• 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.

As noted in #1402, marion_diffuse uses too much memory (and takes a long time) when the surface_tilt input is a vector, which occurs for trackers. marion_diffuse calls marion_integrate for each surface_tilt and each region (sky, horizon or ground). marion_integrate integrates over a solid angle for each of those. @cwhanse, @markcampanelli, @adriesse, @kdebrab and others discussed how this issue could be circumvented via interpolation.

This PR proposes a solution via a new function called marion_diffuse_tracking, to be used instead of marion_diffuse for trackers.

marion_diffuse_tracking creates vectors of 0-90 or 0-180º (for bifacial arrays) with a 0.5º step, marion_integrate is called for each angle in that vector and each region only once, and the results are used to create an interpolator function (using scipy.interpolate.PchipInterpolator). This interpolator is then used to obtain an IAM value for each value in surface_tilt.

Regardless of the length of surface_tilt, marion_integrate is called a fixed number of times, helping to contain the computational burden.

In addition, the interpolation is cached, meaning that e.g. for a PVSystem object with multiple tracker Arrays which have identical tilt vectors, the interpolator will only have to be created once for each region.

More details on time and memory savings, as well as deviation, here. For one example for a 8760-length vector, marion_diffuse_tracking uses 4% of the max memory and 2% of the time of marion_diffuse.

[ramaroesilva]

kudos to @kandersolar for raising the issue and everyone else who contributed to the conceptualization of this.

was now thinking that neither me and @cbcrespo have never modelled bifacial, and we have no idea how IAM is calculated for that (not even sure if marion_diffuse is applicable). would be important that someone who's experienced in this to take a look.

marion_diffuse_tracking creates vectors of 0-90 or 0-180º with a 0.5º step

likely it's not as simple as this....