close pvlib#1041 prune old comments

* remove commented code, especially any code relating to use of arctan2
  adopted in 1fb82cc
* remove comments that should be sourced directly from refernces, and
  instead list reference and Eqs. numbers
* tidy up, combine lines, etc

Author: mikofski

pvlib/tracking.py

@@ -21,7 +21,7 @@ class SingleAxisTracker(PVSystem):
 
     axis_azimuth : float, default 0
         A value denoting the compass direction along which the axis of
-        rotation lies. Measured in decimal degrees East of North.
+        rotation lies. Measured in decimal degrees east of north.
 
     max_angle : float, default 90
         A value denoting the maximum rotation angle, in decimal degrees,
@@ -291,7 +291,7 @@ def singleaxis(apparent_zenith, apparent_azimuth,
 
     axis_azimuth : float, default 0
         A value denoting the compass direction along which the axis of
-        rotation lies. Measured in decimal degrees East of North.
+        rotation lies. Measured in decimal degrees east of north.
 
     max_angle : float, default 90
         A value denoting the maximum rotation angle, in decimal degrees,
@@ -364,58 +364,24 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     if apparent_azimuth.ndim > 1 or apparent_zenith.ndim > 1:
         raise ValueError('Input dimensions must not exceed 1')
 
-    # Calculate sun position x, y, z using coordinate system as in [1], Eq 2.
+    # Calculate sun position x, y, z using coordinate system as in [2], Eq 1.
 
-    # Positive y axis is oriented parallel to earth surface along tracking axis
-    # (for the purpose of illustration, assume y is oriented to the south);
-    # positive x axis is orthogonal, 90 deg clockwise from y-axis, and parallel
-    # to the earth's surface (if y axis is south, x axis is west);
-    # positive z axis is normal to x, y axes, pointed upward.
-
-    # Equations in [1] assume solar azimuth is relative to reference vector
-    # pointed south, with clockwise positive.
-    # Here, the input solar azimuth is degrees East of North,
-    # i.e., relative to a reference vector pointed
-    # north with clockwise positive.
-
-    # NOTE: Equations in [2] agree with the reference frame used here, so
-    # adjustments are not required
-
-    # Rotate sun azimuth to coordinate system as in [1, 2]
-    # to calculate sun position.
-
-    # NOTE: sin(90-x) = cos(x) & cos(90-x) = sin(x)
+    # NOTE: solar elevation = 90 - solar zenith, then use trig identities:
+    # sin(90-x) = cos(x) & cos(90-x) = sin(x)
     sin_zenith = sind(apparent_zenith)
     x = sin_zenith * sind(apparent_azimuth)
     y = sin_zenith * cosd(apparent_azimuth)
     z = cosd(apparent_zenith)
 
-    # translate array azimuth from compass bearing to [1] coord system
-    # wholmgren: strange to see axis_azimuth calculated differently from az,
-    # (not that it matters, or at least it shouldn't...).
-
-    # NOTE: Coordinate system in [2] agrees with refernece frame used here, so
-    # adjustments are not required
-
-    # translate input array tilt angle axis_tilt to [1] coordinate system.
-
-    # In [1] coordinates, axis_tilt is a rotation about the x-axis.
-    # For a system with array azimuth (y-axis) oriented south,
-    # the x-axis is oriented west, and a positive axis_tilt is a
-    # counterclockwise rotation, i.e, lifting the north edge of the panel.
-    # Thus, in [1] coordinate system, in the northern hemisphere a positive
-    # axis_tilt indicates a rotation toward the equator,
-    # whereas in the southern hemisphere rotation toward the equator is
-    # indicated by axis_tilt<0.  Here, the input axis_tilt is
-    # always positive and is a rotation toward the equator.
-
-    # Calculate sun position (xp, yp, zp) in panel-oriented coordinate system:
-    # positive y-axis is oriented along tracking axis at panel tilt;
-    # positive x-axis is orthogonal, clockwise, parallel to earth surface;
-    # positive z-axis is normal to x-y axes, pointed upward.
-    # Calculate sun position (xp,yp,zp) in panel coordinates using [1] Eq 11
-    # note that equation for yp (y' in Eq. 11 of Lorenzo et al 2011) is
-    # corrected, after conversation with paper's authors.
+    # Assume the tracker reference frame is right-handed. Positive y-axis is
+    # oriented along tracking axis tilted by the axis tilt and rotated
+    # clockwise from north by the axis azimuth; positive x-axis is orthogonal,
+    # 90 deg clockwise from y-axis, and parallel to the earth's surface (if y-
+    # axis is south, x-axis is west); positive z-axis is normal to x, y axes,
+    # pointed upward.
+
+    # Calculate sun position (xp, yp, zp) in tracker coordinate system using
+    # [2] Eq 4.
 
     cos_axis_azimuth = cosd(axis_azimuth)
     sin_axis_azimuth = sind(axis_azimuth)
@@ -432,63 +398,36 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     # The ideal tracking angle wid is the rotation to place the sun position
     # vector (xp, yp, zp) in the (y, z) plane; i.e., normal to the panel and
     # containing the axis of rotation.  wid = 0 indicates that the panel is
-    # horizontal.  Here, our convention is that a clockwise rotation is
+    # horizontal. Here, our convention is that a clockwise rotation is
     # positive, to view rotation angles in the same frame of reference as
-    # azimuth.  For example, for a system with tracking axis oriented south,
-    # a rotation toward the east is negative, and a rotation to the west is
-    # positive.
-
-    # angle from x-y plane to projection of sun vector onto x-z plane
-    #     tmp = np.degrees(np.arctan(zp/xp))
-
-    # Obtain wid by translating tmp to convention for rotation angles.
-    # Have to account for which quadrant of the x-z plane in which the sun
-    # vector lies.  Complete solution here but probably not necessary to
-    # consider QIII and QIV.
-    #     wid = pd.Series(index=times)
-    #     wid[(xp>=0) & (zp>=0)] =  90 - tmp[(xp>=0) & (zp>=0)]  # QI
-    #     wid[(xp<0)  & (zp>=0)] = -90 - tmp[(xp<0)  & (zp>=0)]  # QII
-    #     wid[(xp<0)  & (zp<0)]  = -90 - tmp[(xp<0)  & (zp<0)]   # QIII
-    #     wid[(xp>=0) & (zp<0)]  =  90 - tmp[(xp>=0) & (zp<0)]   # QIV
-
-    # NOTE: Use arctan2 and avoid the tmp corrections.
+    # azimuth. For example, for a system with tracking axis oriented south, a
+    # rotation toward the east is negative, and a rotation to the west is
+    # positive. This is a right-handed rotation around the tracker y-axis.
 
     # Calculate angle from x-y plane to projection of sun vector onto x-z plane
-    # and then obtain wid by translating tmp to convention for rotation angles.
+    # using [2] Eq. 5.
 
-    # NOTE: if zp/xp = tan(90 - wid) = cot(wid) then tan(wid) = xp/zp
     wid = np.degrees(np.arctan2(xp, zp))
 
     # filter for sun above panel horizon
     zen_gt_90 = apparent_zenith > 90
     wid[zen_gt_90] = np.nan
 
-    # Account for backtracking; modified from [1] to account for rotation
-    # angle convention being used here.
+    # Account for backtracking
     if backtrack:
         # distance between rows in terms of rack lengths relative to side slope
         axes_distance = 1/gcr/cosd(side_slope)
-        # clip needed for low angles. GH 656
-        # temp = np.clip(axes_distance*cosd(wid - side_slope), -1, 1)
 
         # NOTE: account for rare angles below array, see GH 824
         temp = np.abs(axes_distance * cosd(wid - side_slope))
 
-        # backtrack angle
-        # (always positive b/c acosd returns values between 0 and 180)
-        # wc = np.degrees(np.arccos(temp))
-        # equation 14, ref [2]
+        # backtrack angle using [2], Eq. 14
         with np.errstate(invalid='ignore'):
             wc = np.degrees(-np.sign(wid)*np.arccos(temp))
 
-        # Eq 4 applied when wid in QIV (wid < 0 evalulates True), QI
-        # with np.errstate(invalid='ignore'):
-        #     errstate for GH 622
-        #     tracker_theta = np.where(wid < 0, wid + wc, wid - wc)
-
         # NOTE: in the middle of the day, arccos(temp) is out of range because
         # there's no row-to-row shade to avoid, & backtracking is unnecessary
-        # Equations 15-16, ref [2]
+        # [2], Eqs. 15-16
         with np.errstate(invalid='ignore'):
             tracker_theta = wid + np.where(temp < 1, wc, 0)
     else:
@@ -498,10 +437,10 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     # system-plane normal
     tracker_theta = np.clip(tracker_theta, -max_angle, max_angle)
 
-    # calculate panel normal vector in panel-oriented x, y, z coordinates.
-    # y-axis is axis of tracker rotation.  tracker_theta is a compass angle
+    # Calculate panel normal vector in panel-oriented x, y, z coordinates.
+    # y-axis is axis of tracker rotation. tracker_theta is a compass angle
     # (clockwise is positive) rather than a trigonometric angle.
-    # the *0 is a trick to preserve NaN values.
+    # NOTE: the *0 is a trick to preserve NaN values.
     panel_norm = np.array([sind(tracker_theta),
                            tracker_theta*0,
                            cosd(tracker_theta)])
@@ -512,99 +451,49 @@ def singleaxis(apparent_zenith, apparent_azimuth,
     # calculate angle-of-incidence on panel
     aoi = np.degrees(np.arccos(np.abs(np.sum(sun_vec*panel_norm, axis=0))))
 
-    # calculate panel tilt and azimuth
-    # in a coordinate system where the panel tilt is the
-    # angle from horizontal, and the panel azimuth is
-    # the compass angle (clockwise from north) to the projection
-    # of the panel's normal to the earth's surface.
-    # These outputs are provided for convenience and comparison
-    # with other PV software which use these angle conventions.
+    # Calculate panel tilt and azimuth in a coordinate system where the panel
+    # tilt is the angle from horizontal, and the panel azimuth is the compass
+    # angle (clockwise from north) to the projection of the panel's normal to
+    # the earth's surface. These outputs are provided for convenience and
+    # comparison with other PV software which use these angle conventions.
 
-    # project normal vector to earth surface.
-    # First rotate about x-axis by angle -axis_tilt so that y-axis is
-    # also parallel to earth surface, then project.
+    # Project normal vector to earth surface. First rotate about x-axis by
+    # angle -axis_tilt so that y-axis is also parallel to earth surface, then
+    # project.
 
     # Calculate standard rotation matrix
     rot_x = np.array([[1, 0, 0],
                       [0, cosd(-axis_tilt), -sind(-axis_tilt)],
                       [0, sind(-axis_tilt), cosd(-axis_tilt)]])
 
-    # panel_norm_earth contains the normal vector
-    # expressed in earth-surface coordinates
-    # (z normal to surface, y aligned with tracker axis parallel to earth)
+    # panel_norm_earth contains the normal vector expressed in earth-surface
+    # coordinates (z normal to surface, y aligned with tracker axis parallel to
+    # earth)
     panel_norm_earth = np.dot(rot_x, panel_norm).T
 
-    # projection to plane tangent to earth surface,
-    # in earth surface coordinates
+    # projection to plane tangent to earth surface, in earth surface
+    # coordinates
     projected_normal = np.array([panel_norm_earth[:, 0],
                                  panel_norm_earth[:, 1],
                                  panel_norm_earth[:, 2]*0]).T
 
     # calculate vector magnitudes
     projected_normal_mag = np.sqrt(np.nansum(projected_normal**2, axis=1))
 
-    # renormalize the projected vector
-    # avoid creating nan values.
+    # renormalize the projected vector, avoid creating nan values.
     non_zeros = projected_normal_mag != 0
     projected_normal[non_zeros] = (projected_normal[non_zeros].T /
                                    projected_normal_mag[non_zeros]).T
 
     # calculation of surface_azimuth
-    # 1. Find the angle.
-    #     surface_azimuth = pd.Series(
-    #         np.degrees(np.arctan(projected_normal[:,1]/projected_normal[:,0])),
-    #         index=times)
     surface_azimuth = \
         np.degrees(np.arctan2(projected_normal[:, 1], projected_normal[:, 0]))
 
-    # 2. Clean up atan when x-coord or y-coord is zero
-    #     surface_azimuth[(projected_normal[:,0]==0) & (projected_normal[:,1]>0)] =  90
-    #     surface_azimuth[(projected_normal[:,0]==0) & (projected_normal[:,1]<0)] =  -90
-    #     surface_azimuth[(projected_normal[:,1]==0) & (projected_normal[:,0]>0)] =  0
-    #     surface_azimuth[(projected_normal[:,1]==0) & (projected_normal[:,0]<0)] = 180
-
-    # 3. Correct atan for QII and QIII
-    #     surface_azimuth[(projected_normal[:,0]<0) & (projected_normal[:,1]>0)] += 180 # QII
-    #     surface_azimuth[(projected_normal[:,0]<0) & (projected_normal[:,1]<0)] += 180 # QIII
-
-    # 4. Skip to below
-
-    # at this point surface_azimuth contains angles between -90 and +270,
-    # where 0 is along the positive x-axis,
-    # the y-axis is in the direction of the tracker azimuth,
-    # and positive angles are rotations from the positive x axis towards
-    # the positive y-axis.
-    # Adjust to compass angles
-    # (clockwise rotation from 0 along the positive y-axis)
-    #    surface_azimuth[surface_azimuth<=90] = 90 - surface_azimuth[surface_azimuth<=90]
-    #    surface_azimuth[surface_azimuth>90] = 450 - surface_azimuth[surface_azimuth>90]
-
-    # finally rotate to align y-axis with true north
-    # PVLIB_MATLAB has this latitude correction,
-    # but I don't think it's latitude dependent if you always
-    # specify axis_azimuth with respect to North.
-    #     if latitude > 0 or True:
-    #         surface_azimuth = surface_azimuth - axis_azimuth
-    #     else:
-    #         surface_azimuth = surface_azimuth - axis_azimuth - 180
-    #     surface_azimuth[surface_azimuth<0] = 360 + surface_azimuth[surface_azimuth<0]
-
-    # the commented code above is mostly part of PVLIB_MATLAB.
-    # My (wholmgren) take is that it can be done more simply.
-    # Say that we're pointing along the postive x axis (likely west).
-    # We just need to rotate 90 degrees to get from the x axis
-    # to the y axis (likely south),
-    # and then add the axis_azimuth to get back to North.
-    # Anything left over is the azimuth that we want,
-    # and we can map it into the [0,360) domain.
-
-    # 4. Rotate 0 reference from panel's x axis to it's y axis and
-    #    then back to North.
+    # Rotate 0 reference from panel's x-axis to its y-axis and then back to
+    # north.
     surface_azimuth = 90 - surface_azimuth + axis_azimuth
 
-    # 5. Map azimuth into [0,360) domain.
-    # surface_azimuth[surface_azimuth < 0] += 360
-    # surface_azimuth[surface_azimuth >= 360] -= 360
+    # Map azimuth into [0,360) domain.
     with np.errstate(invalid='ignore'):
         surface_azimuth = surface_azimuth % 360