Fix calc_surface_orientation broadcasting for >1-D tracker_theta

[ArthurOnnoTerabase]

The unit-normal-based implementation introduced in v0.15.1 (#2702) used np.column_stack inside _unit_normal, which only handles 1-D inputs correctly. For 2-D (or higher rank) tracker_theta the column-stacked array had the wrong shape, and the subsequent np.where in calc_surface_orientation raised a ValueError.

Switch to np.stack(..., axis=-1) and use ellipsis indexing in calc_surface_orientation so the function preserves the input rank for arbitrary-shape tracker_theta. np.atleast_1d preserves the historical (1, 3) return shape of _unit_normal for scalar input.

Also adds a parametrised regression test that exercises 2-D and 3-D tracker_theta inputs.

Closes #2747

• Closes calc_surface_orientation broadcasting error for 2-D inputs in 0.15.1 (regression vs 0.15.0) #2747
• 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.

[cwhanse]

Thanks @ArthurOnnoTerabase