Quantifying the dynamic wing morphing of hovering hummingbird

Animal wings are lightweight and flexible; hence, during flapping flight their shapes change. It has been known that such dynamic wing morphing reduces aerodynamic cost in insects, but the consequences in vertebrate flyers, particularly birds, are not well understood. We have developed a method to reconstruct a three-dimensional wing model of a bird from the wing outline and the feather shafts (rachides). The morphological and kinematic parameters can be obtained using the wing model, and the numerical or mechanical simulations may also be carried out. To test the effectiveness of the method, we recorded the hovering flight of a hummingbird (Amazilia amazilia) using high-speed cameras and reconstructed the right wing. The wing shape varied substantially within a stroke cycle. Specifically, the maximum and minimum wing areas differed by 18%, presumably due to feather sliding; the wing was bent near the wrist joint, towards the upward direction and opposite to the stroke direction; positive upward camber and the ‘washout’ twist (monotonic decrease in the angle of incidence from the proximal to distal wing) were observed during both half-strokes; the spanwise distribution of the twist was uniform during downstroke, but an abrupt increase near the wrist joint was found during upstroke.

on the wing surface , whereas in some cases natural features or projected random dots are used. In the present study, the direct contact with the bird was prohibited by the zoo. We therefore decided to use the feather shafts (rachides or rachises) as the major feature on the wing surface. They are supplemented by the boundary between primary and secondary feathers, and the boundary between covert feathers and the other parts. These features are not sharp enough for the automatic reconstruction with the adequate accuracy, partly because of the motion blur due to the slow shutter speed (1/3000 s), which was inevitable because we dependent upon the natural lighting (the addition of any artificial lighting was prohibited). Therefore we needed to track the natural features manually, which required us to select limited number of image frames to be analysed.
Consequently, the wingbeat cycle for the detailed analysis was selected based on the preliminary analysis checking the wingbeat frequency for each stroke cycle. Four separate stable hover-drink sequences were analysed (figure S2) each consisting of 45, 26, 12, and 29 wingbeat cycles. Figure S2A shows the box plots and swarm plots for the image frames per cycle (dividing 2000 frames per seconds by this number yields the wingbeat frequency). Figure S2B-E) illustrate the time history of wingbeat frequency. It seems the bird tends to increase the wingbeat frequency when it ends drinking (except for the 4th sequence) but overall the frequency was within the range of 28-32 Hz. The first sequence (figure S2B) was the longest of them, and there was a five continuous wingbeats with the same frequency at 29.4 Hz before the end of drinking. We selected the middle of them for the further analysis for wing kinematics and deformation (annotated by orange texts in figure S2A and B; also highlighted by the green stripe in figure S2B). This is the median frequency of the 1st sequence and very close to the median frequency for the pooled data (29.9 Hz, figure S2A). We also tracked the wingbase and wingtip for the stable cycles around the selected cycle (near the green stripe in S2A). The calculated wingbeat amplitudes were within 105 and 115 degrees. For the selected cycle, we did much more detailed analysis, which finally resulted in slightly different values mentioned above (28.8 Hz and 103 deg, Table 1).

Determination of the global coordinates
To calculate the stroke plane and stroke plane coordinates, the global coordinates need to be obtained beforehand.
The global Z axis (vertical upward) was obtained from the metallic ruler, therefore we have the horizontal plane. To obtain the other two axes, namely, backward (global X axis) and lateral right (global Y axis), we did the following. The left wingtip was tracked for four time frames (t = 0.077T, 0.19T, 0.36T, 0.54T ). The vector connecting the left and right wingtips were averaged. Since the left wingtip was slightly lower than the right wingtip (the left-right wingtip vector was tilted from the horizontal plane by approximately 1.5 deg), we did not use this directly as the global Y axis. Instead, we called this left-right wingtip vector a Y' axis and assumed the yaw component of the Y' axis is correct but the roll component is not. Therefore, the global X axis was determined as the perpendicular direction to the Z and Y' axes. Next, the Y axis was determined as the perpendicular direction to the Z and X axes.
The cause of the slightly lower left wingtip compared to the right wingtip can be either one or combination of the followings: the error in the metallic ruler edge direction from the true gravitational acceleration vector, the calibration error of the left wingtip (since the left wingtip is nearly at the boundary of the calibration volume), or the slight asymmetry of the wing kinematics.

Line segments and triangle mesh
The line segments and the triangular elements (triangle mesh) are shown in figure S3 for the illustration purpose.

Shortest path length and wing length
The difference between shortest path (and its length) and the wing length are illustrated in figure S4. The theoretical minimum shortest-path length L SP is achieved when the wing is totally flat, when the shortest-path length is equal to the wing length L W . Otherwise, where there is any wing bending, the shortest-path length is greater than the wing length.

Selection of the stroke plane and stroke plane angle
The stroke plane (and therefore stroke plane angle) depends on the choice of the characteristic point used for the liner regression. If one choose another point than the wingtip, stroke plane would have been different. We checked this by using the points along the shortest-path in 5% L SP spacings (figure S5). The results shows the stroke plane angles are within 2 deg variation (minimum 10.9 deg to maximum 12.9 deg). Therefore we chose wingtip for defining the stroke plane angle, which has been the standard for the previous literatures for flapping flight. Note here the stroke plane concept is connected to the flapping wing aerodynamic model, where wing  Figure S3: Example of the line segments (A) and triangle mesh (B). In A, wing outline is coloured by red and orange; rachides are coloured by light blue; the covert edge is coloured by dark green and light green; and the primary-secondary boundary is coloured by black (dual colours are for the illustration purpose to distinguish individual line segment). The wing area for each region and for the overall wing was calculated using the triangle mesh in B.

Camber at mid-chord
The camber in the main text was defined as the maximum height h max for each wing cross-section divided by the local chord length c. Here we also calculated the camber using the height at mid-chord (figure S6, table 0.59T

Shortest path
Wing length Figure S4: Shortest path length is the length along the shortest path (magenta curve), which is always greater than or equal to the and the wing length (black line with arrowheads). In figure S7 a detailed schematic for how to obtain AoA. A 3D relative wind vector (not shown) firstly obtained from the central difference from the P n at last and next time steps (blue and green wings, respectively).
Then the component parallel to the wing cross-section at current time step is extracted, and finally the AoA is calculated. Similarly, from the component perpendicular to the cross-section the SSA is calculated (not shown).
The time history of SSA is shown in figure S8. The SSA for the wing distal sections were kept either 0 deg  or ± 180 deg for all the time except during stroke reversals (especially supination). Wing proximal sections, especially 10% section experienced the nonzero values at early half-strokes, but as we saw in S9 the relative wind speed (magnitude of the relative wind vector) itself is small for the current flight (hovering) and hence side force or the consequent aerodynamic power would not be large.

Relative wind speed and local chord length
Relative wind speed (magnitude of the relative wind vector for each wing cross-section) and the local chord length at wing sections are shown in figures S9 and S10. It is evident the wing distal sections experience much Figure S10: Local chord length. Gray shaded region indicates downstroke period.
higher relative wind speed than the proximal sections while having slightly shorter chord lengths. For example, the 60% section has roughly 80% chord length of the 20% section but the peak relative wind speed reached around four times larger. Since the amount of aerodynamic force is proportional to the square of the wind speed (and linearly proportional to the chord length), the wing distal sections would be the major contributor for the aerodynamic force. Note the difference in the angle of attack (AoA) and wing shape (camber) at different sections would surely affect the aerodynamic force coefficient and hence aerodynamic force, but that won't be enough to compensate the impact from the speed.