Fix single axis tracker algorithm

doc/src/records/pvarray.md

@@ -102,6 +102,16 @@ The tilt of the photovoltaic array from horizontal.  Values outside the range 0
                             else 0
   ------------------------------------------------------------------
 
+The following figures illustrate the use of both pvTilt and pvAzm for various configurations:
+
+![Fixed, south facing, tilted at 40^o^](media/pv_fixed.png)
+
+![One-axis tracker, south facing, tilted at 20^o^](media/pv_tilted_tracker_south.png)
+
+![One-axis tracker, horizontal aligned North/South (more common)](media/pv_horiz_tracker_south.png)
+
+![One-axis tracker, horizontal aligned East/West (less common)](media/pv_horiz_tracker_east.png)
+
 **pvAzm=*float***
 
 Photovoltaic array azimuth (0 = north, 90 = east, etc.). If a value outside the range 0^o^ $\leq$ *x* $<$ 360^o^ is given, it is normalized to that range. For one-axis tracking, defines the azimuth of the rotation axis. Not used for two-axis tracking arrays. Should be omitted if pvVertices is given.

src/CNRECS.DEF

@@ -4454,6 +4454,7 @@ RECORD PVARRAY "PVArray" *RAT		// input / runtime photovoltaics array
   *h *e FLOAT pv_aoi		// angle of incidence (radians)
   *h *e ANGLE pv_panelTilt	// tilt of pv panel (different from array tilt for tracking systems), radians
   *h *e ANGLE pv_panelAzm	// azimuth of pv panel (different from array tilt for tracking systems), radians
+  *h *e ANGLE pv_panelRot	// rotation of pv panel for 1-axis tracking systems, radians clockwise from vertical
   *h *e FLOAT pv_poa		// plane of array incidence (before shading), Btu/h-ft2
   *h *e FLOAT pv_poaBeam	// plane of array beam incidence (before shading), Btu/h-ft2
   *h *e FLOAT pv_radIBeam	// beam radiation incident on array, Btu/h-ft2

src/PVCalc.cpp

@@ -303,35 +303,44 @@ RC PVARRAY::pv_CalcPOA()
 	{
 	case C_PVARRCH_1AXT:
 	{
-		float sinz = sin(acos(cosz));
-		float x = (sinz*sin(azm - pv_azm)) / (sinz*cos(azm - pv_azm)*sin(pv_tilt) + cosz*cos(pv_tilt));
-		float psi;
-		if (x < 0.f && (azm - pv_azm) > 0.f) {
+		// Based on Rotation Angle for the Optimum Tracking of One-Axis Trackers by William F.Marion and Aron P.Dobos
+		const float sinz = sin(acos(cosz));  // sin of zenith angle
+		float azm_delta = azm - pv_azm;
+
+		// normalize azm_delta between -Pi and Pi
+		if (azm_delta >= kPi)
+			azm_delta = azm_delta - k2Pi;
+
+		if (azm_delta <= -kPi)
+			azm_delta = azm_delta + k2Pi;
+
+		const float x = (sinz*sin(azm_delta)) / (sinz*cos(azm_delta)*sin(pv_tilt) + cosz*cos(pv_tilt));  // from Equation #7
+		float psi = 0.f;
+
+		if (x < 0.f && azm_delta > 0.f) {
 			psi = kPi;
 		}
-		else if (x > 0.f && (azm - pv_azm) < 0.f){
+		else if (x > 0.f && azm_delta < 0.f){
 			psi = -kPi;
 		}
-		else {
-			psi = 0.f;
-		}
-		float r = atan(x) + psi;
+
+		pv_panelRot = atan(x) + psi;  // Equation #7
+
 		const float rlim = 0.25f*kPi;
-		r = max(-rlim, min(rlim, r));
-		pv_panelTilt = acos(cos(r)*cos(pv_tilt));
+		pv_panelRot = max(-rlim, min(rlim, pv_panelRot));
+		pv_panelTilt = acos(cos(pv_panelRot)*cos(pv_tilt));  // Equation #1
 		if (pv_panelTilt == 0.f) {
 			pv_panelAzm = pv_azm;
 		}
 		else
 		{
-			float rx = bracket(-1.f, sin(r) / sin(pv_panelTilt), 1.f);
-			float asrx = asin(rx);
-			if (r >= -kPiOver2 && r <= kPiOver2)
-				pv_panelAzm = pv_azm + asrx;
-			else if (r >= -kPi && r < -kPiOver2)
-				pv_panelAzm = pv_azm - asrx - kPi;
-			else // if (r > kPiOver2 && r <= kPi)
-				pv_panelAzm = pv_azm - asrx + kPi;
+			float const asrx = asin(bracket(-1.f, sin(pv_panelRot) / sin(pv_panelTilt), 1.f));
+			if (pv_panelRot >= -kPiOver2 && pv_panelRot <= kPiOver2)
+				pv_panelAzm = pv_azm + asrx;        // Equation #2
+			else if (pv_panelRot >= -kPi && pv_panelRot < -kPiOver2)
+				pv_panelAzm = pv_azm - asrx - kPi;  // Equation #3
+			else // if (pv_panelRot > kPiOver2 && pv_panelRot <= kPi)
+				pv_panelAzm = pv_azm - asrx + kPi;  // Equation #4
 		}
 	}
 		break;