## ABSTRACT

Animal flight requires aerodynamic power, which is challenging to determine accurately *in vivo*. Existing methods rely on approximate calculations based on wake flow field measurements, inverse dynamics approaches, or invasive muscle physiological recordings. In contrast, the external mechanical work required for terrestrial locomotion can be determined more directly by using a force platform as an ergometer. Based on an extension of the recent invention of the aerodynamic force platform, we now present a more direct method to determine the *in vivo* aerodynamic power by taking the dot product of the aerodynamic force vector on the wing with the representative wing velocity vector based on kinematics and morphology. We demonstrate this new method by studying a slowly flying dove, but it can be applied more generally across flying and swimming animals as well as animals that locomote over water surfaces. Finally, our mathematical framework also works for power analyses based on flow field measurements.

## INTRODUCTION

The energy required to locomote on land, in water and in the air is a critical parameter for understanding the muscle physiology, comparative biomechanics and movement ecology of animals (Alexander, 1999; Alexander, 2003). The external work required to locomote has been analyzed extensively in terrestrial locomotion (Biewener, 1992; Biewener, 2003; Biewener et al., 1988; Gordon et al., 2017), in particular since force platforms started to be used as ergometers (Cavagna, 1975). The external aerodynamic or hydrodynamic power required to locomote in fluids is less straightforward to determine (Alexander, 2003; Biewener and Patek, 2018; Dudley, 2002; Pennycuick, 2008; Videler, 2012; Vogel, 1994). Flapping flight is a particularly energetically demanding form of locomotion, and accurately determining the aerodynamic power required to fly has proven challenging (Alexander, 2003; Biewener and Patek, 2018; Dudley, 2002; Pennycuick, 2008).

The most widely used method (Chin et al., 2017) to determine the aerodynamic power required for flight is probably aerodynamic modelling of flapping wings based on quasi-steady blade element models (Ellington, 1984a; Pennycuick, 2008; Weis-Fogh, 1972, 1973). These models have been improved using dynamically scaled robotic models, which rely on simplifying the body kinematics, morphology, and fluid-structure interaction of animals (Dickinson et al., 1999; Dickson et al., 2008). Another valuable modelling approach is using computational fluid dynamics models based on kinematics models or measurements (Liu, 2002). All these approaches involve computing the aerodynamic force based on a significant number of assumptions that have yet to be better validated with *in vivo* experiments. Unfortunately, current *in vivo* methods for determining aerodynamic power also have significant limitations that constrain their ability to rigorously validate models. Aerodynamic power calculations based on wake flow measurements (Muijres et al., 2011; Von Busse et al., 2014) assume the power can be calculated based on measuring the wake in a transverse plane behind or below the animal. By making this assumption, they only consider the “outflow” surface of the control volume and neglect contributions from the rest of the volume, which could be significant (Lentink et al., 2015). Another approach is inverse dynamics, in which measured body kinematics and morphology are used to determine both the velocity and acceleration distribution of the body. This requires taking the second time-derivative of kinematics, which is known to be significantly noisier than the first derivative needed to determine the velocity distribution (Hedrick et al., 2003; Ros et al., 2011, 2015). Further, the inertia tensor distribution needs to be determined, which currently requires sacrificing the animal after the measurements. Flight muscle work loop measurements, on the other hand, are more informative in terms of determining the total power required to move the wing (Biewener et al., 1992; Tobalske et al., 2003). Accurate calibration of the muscle force measurement based on surgically implanted strain gauges remains challenging however (Biewener, 2011; Biewener and Patek, 2018; Biewener et al., 1992; Jackson et al., 2011; Tobalske and Biewener, 2008), and the invasive nature of the experiment can change the animal's flight behavior (Tobalske et al., 2005). Finally, muscle recordings alone cannot be used to establish the external aerodynamic power required to fly, because they include the additional power to overcome friction and inertia in the musculoskeletal system. Despite limitations, these methods have all been instrumental for developing our present understanding of the energetics of animal flight. This understanding can be advanced further if there is a better and more convenient method for determining the aerodynamic power *in vivo* that relies on smaller assumptions. Here we present a new method to determine the aerodynamic power required to fly based on taking the vector dot product of the aerodynamic force vector on the wing with the representative velocity vector of the wing. We directly measure the wing velocity and the vertical and horizontal aerodynamic forces, and calculate the lateral aerodynamic force based on a first-principles derivation. We demonstrate this simple and accurate method by analyzing the aerodynamic power that a dove exerts to fly slowly.

## MATERIALS AND METHODS

Here, we present an accurate method for measuring the time-resolved aerodynamic power a bird generates in flight as demonstrated for a single perch to perch flight of a dove (Fig. 1A). The basic idea of this method relies on the fact that the total aerodynamic power for both wings *P*_{aero} equals two times the force vector for each wing dotted with the representative velocity vector for each wing :
(1)

As we describe in the following sections, we directly measure the representative velocity vector using wing kinematics (Fig. 1C,D) and directly measure horizontal and vertical forces to recover two of the three components of the force vector (Fig. 2A). By combining these measured quantities, we can calculate lift and drag (Fig. 2B), which allows us to find the third component of wing force (Fig. 2D). Together, these measurements combine to produce a direct method for computing aerodynamic power (Fig. 3).

### Animal experiment

The experimental setup consisted of a structured-light system (Fig. 1A) to image the 3D surface of the dove (Fig. 1B), synchronized with an aerodynamic force platform (AFP; Fig. 1A) to measure the vertical and horizontal aerodynamic forces produced by the dove (Fig. 2A). We analyzed the second wingbeat after takeoff of a 2-year-old male near-white Ringneck dove [*Streptopelia risoria* (Linnaeus 1758) 140–142 g, 0.48 m wingspan], which was trained to fly between perches 0.65 m apart inside the AFP. An open window in the acrylic side wall of the AFP was used to bring the dove in and out of the AFP. The window was closed with an additional acrylic panel prior during each recorded flight. The perches were mounted 0.36 m above the bottom plate of the AFP and the residual descent angle between the takeoff and landing perch was 2 deg. Training involved a combination of positive reinforcement and light tapping on the tail to cue a flight. All experiments were in accordance with Stanford University's Institutional Animal Care and Use Committee.

### Wing velocity calculation based on kinematics measurements

We calculated the time-resolved representative velocity vector of the dove's wing (Fig. 1C) based on 3D kinematics measurements (Fig. 1B). The kinematics were measured by tracking marker points on the dove as well as by measuring the 3D surface of the dove's body, tail and left wing using our 3D calibrated (Deetjen and Lentink, 2018) structured light method (Deetjen et al., 2017) at 1000 frames s^{−1} (Fig. 1A,B). Additionally, we manually tracked the following feature points using triangulation: ninth primary wingtip, shoulder, wrist, middle of the back, left and right feet, left and right eyes, tip and base of the beak, and top of the head (Fig. 1B). For the wingtips, shoulder, wrist, back and the top of the head, we attached square retro-reflective marker tape and identified their centers when they were visible. The remaining frames were estimated with a combination of manual estimation and interpolation. We combined these collected kinematics and surface data to fit a smooth morphing surface to the body, tail, and wings of the dove (see Appendix 1 for more details). Finally, we assumed bilateral symmetry for the right wing, because we focused our cameras and projectors on the left wing to maximize resolution and because the flight path was approximately straight (6.8 deg off the centerline based on tracking kinematics points on the body and head). In Appendix 3, the yaw angle of 6.8 deg is taken into account to increase the accuracy of the resulting aerodynamic power. The combined wing surface reconstruction and point tracking enable a high-fidelity computation of the representative wing velocity vector based on integrating the distributed velocity along the entire wing (Fig. 1C).

The primary method we used for calculating the representative wing velocity uses the velocity vector at the quarter-chord location of blade elements (25% chord length behind the leading edge). We used this location because we estimated that the center of pressure is located at the quarter-chord and thus, the aerodynamic force acts at this location. For thin-airfoils at low angles of attack, the center of pressure is assumed to be located at the quarter-chord (Pennycuick, 2008). For flapping wings, the angle of attack varies and attains high values. Hence, the center of pressure, which is a function of the angle of attack, shifts in the direction of 50% chord length (Dickson et al., 2008; Himmelskamp, 1947). The precise location of the center of pressure influences the aerodynamic power mid-stroke as shown in Fig. S1. To evaluate this, we calculated the average angle of attack during the downstroke of the dove in our study, 49±13 deg, and found the estimated center of pressure on the dove's wings shifted to roughly 27±11% chord. This estimate is based on center of pressure measurements made with a robot flapping wing, which underpins the established quasi-steady model for flapping wings (Dickson et al., 2008) we implemented. It was confirmed by the small difference between the aerodynamic power determined using the quarter-chord method versus the quasi-steady method based on the representative wing velocity vector (Fig. S1). Hence, we used the quarter-chord as a reasonable approximation for the location of the representative velocity vector of each chord following aircraft aerodynamics convention (Drela, 2014; Leishman, 2006).

In total, we compared three methods for computing the representative wing velocity vector (Fig. 1D): (1) Using the velocity vector at the quarter-chord from the surface measurement; (2) using a chord-averaged velocity vector from the surface measurement; (3) only using the velocity of the shoulder and ninth primary from manually tracked kinematics. For all three methods, the representative wing velocity vector used to compute total aerodynamic power (Eqn 1) is given by:
(2)where *S*_{j} is the surface area of the *j*^{th} of *n* blade elements and is the associated velocity vector (derivation in Appendix 2). We use blade elements with a spanwise width of 1 mm to compute the representative wing velocity and to fit the measured surface of the wing. For the three methods used to compute , we define the velocity of each blade element as: (1) the velocity at its quarter-chord, (2) the velocity averaged over the entire chord, or (3) the linearly interpolated velocity based on the measured velocities at the shoulder and ninth primary:
(3)where is the velocity of the ninth primary and is the velocity of the shoulder. The blade elements start at the shoulder (*j* =1) and end at the wingtip (*j* =*n*) and the 0.5 constants are used to estimate the velocity at the spanwise center of each element.

### Aerodynamic force measurement and reconstruction

We determined the time-resolved aerodynamic force vector generated by each dove's wing by measuring the net aerodynamic forces in the horizontal and vertical directions (Fig. 2A). We then combined these 2D forces with our 3D wing kinematics measurements to reconstruct the lateral force and thus 3D force vector (Fig. 2D). The vertical and horizontal aerodynamic forces of the dove were measured using a 2D AFP (Lentink, 2018; Lentink et al., 2015). Previous 1D versions of the AFP (Chin and Lentink, 2016; Ingersoll and Lentink, 2018) measured vertical forces by instrumenting the floor and ceiling of a flight chamber with carbon fiber composite panels. The new 2D AFP used in this study (Fig. 1A) adds two more panels on the front and back sides of the flight chamber (1 m length×1 m height×0.6 m width) to measure horizontal forces as well. Each of the four panels is connected in a statically determined manner to three Nano 43 sensors (six-axis, SI-9-9.125 calibration; ATI Industrial Automation) sampling at 2000 Hz with a resolution of 2 mN. These force plates mechanically integrate the unsteady pressure and shear distributions over its walls, which form a closed mechanical control surface integral of the Navier-Stokes equations over the volume in which the dove flies. Thus, they sum up to the time-resolved net aerodynamic force generated inside by the dove (Lentink, 2018; Lentink et al., 2015). The pressure and shear stress distributions on the walls result from the motion of the dove, which generates airflow inside that is updated with the speed of sound (Hightower et al., 2017; Lentink, 2018; Lentink et al., 2015). The associated phase delay due to the traveling time of the acoustics waves in the volume is on the order of 1.5 ms (for 0.5 m distance from the wall and a speed of sound of 340 m s^{−1}). This time delay is negligible compared with the period of the cut-off frequency of the force plates (80 Hz, 12.5 ms) and the wingbeat period (10 Hz, 100 ms). Consequently the AFP recovers the net aerodynamic force generation with almost no phase delay, similar to terrestrial force platforms, which are also constrained by their natural vibration frequency. We also measured takeoff and landing forces by mounting the perches on carbon fiber beams, each of which was connected in a statically determined manner to three Nano 43 sensors set on mechanically isolated support structures. The force measurements were filtered using an eighth-order Butterworth filter with a cutoff frequency of 80 Hz for the plates and 60 Hz for the perches, or about 8 and 6 times the flapping frequency of a dove, respectively. This enabled us to filter out noise from the setup, because the natural frequencies of the force plates were all above 90 Hz, and the perches had natural frequencies above 70 Hz.

The setup was validated by tethering a quadcopter to an instrumented beam and then comparing vertical and horizontal forces measured by the beam with forces measured by the AFP as reported previously (Lentink et al., 2015). Over a period of 10 s and a sampling rate of 1000 Hz, the impulse ratio and the mean force ratio of the 2D AFP to the beam were both 1.00±0.02 (*n*=10 trials) in the vertical direction and 1.00±0.01 in the horizontal direction (*n*=20 trials). Additionally, for a flight that starts and ends at rest, we expect that the total vertical impulse imparted by the legs and wings should equal full bodyweight (bw) support (Chin and Lentink, 2017), and that the net horizontal impulse should equal zero. Integrating the forces from takeoff to landing for the dove's flight, we measured a vertical impulse ratio (impulse from legs and wings divided by impulse due to bodyweight) of 1.03, and a horizontal impulse (impulse from legs and wings divided by bodyweight) of 0.02 bw·s.

We converted the direct measurement of the net 2D aerodynamic force produced by the dove into the 2D forces generated by each wing by assuming bilateral symmetry. The vertical and horizontal components of the aerodynamic force on the left and right wing, *F*_{wing,z} and *F*_{wing,x}, respectively, equal:
(4)
(5)where *F*_{AFP,z} and *F*_{AFP,x} are the net vertical and horizontal forces directly measured by the AFP. The body forces *F*_{body,z} and *F*_{body,x}, as well as the tail forces *F*_{tail,z} and *F*_{tail,x} are very small (∼0.5±0.4% bw) compared with the wing forces due to the low body velocity (1.2±0.3 m s^{−1}). We estimated these forces by modeling the body as an ellipsoid (Alonso et al., 2010) and the tail as the average between a delta wing and flat plate at low Reynolds numbers (Earnshaw and Lawford, 1964; Taira et al., 2007).

We used the 2D aerodynamic forces on each wing together with the direction of the wing's measured velocity to calculate the third, lateral, aerodynamic force component (Fig. 2D). The 3D aerodynamic force can be decomposed into the lift and drag vector components generated by the wing (Fig. 2B). We then leverage the measured wing velocity to find the directions of these lift and drag vectors:
(6)Here, *D* and *L* are the magnitudes of the drag and lift vectors and **D̂** and **L̂** are the representative unit vector directions for each wing (Eqns 10,11). We can also rewrite Eqn 6 to explicitly show the three equations that are formed:
(7)
(8)
(9)where x, y and z subscripts indicate the x, y or z components of the corresponding vectors in Eqn 6. Using the AFP, we directly measure *F*_{wing,x} and *F*_{wing,z}, and using the measured velocity, we can directly determine the direction of drag, which by definition acts in the opposite direction as velocity:
(10)

Additionally, we know that the direction of lift **L̂** is perpendicular to velocity and we assume that it also acts perpendicular to the radial vector of the wing, because the aerodynamic force is governed by pressure. Thus, the direction of overall lift produced by the wing equals:
(11)where the radial vector is the vector from the shoulder to the ninth primary. The negative sign is used for the left wing and the positive sign for the right wing. When we put this all together, we can solve for the lateral aerodynamic force *F*_{wing,y} in Eqn 8 by solving for *D* and *L* as follows:
(12)
(13)

For greater accuracy in cases where the flight path of the bird is misaligned with the world reference frame, as defined by the AFP, the following equations for yaw adjusted drag, *D*_{A}, and lift, *L*_{A}, should be used in place of *D* and *L*:
(14)
(15)
(16)where θ is the yaw angle of the bird (derivation in Appendix 3; effect of yaw adjustment in Fig. S4). Future studies could improve on this result by measuring the kinematics for both the left and right wings separately to relax the bilateral symmetry constraint. By plugging the magnitude of lift and drag into Eqn 8 (*L* and *D* for small yaw angles or *L*_{A} and *D*_{A} for large yaw angles), we then solved for the lateral aerodynamic force. For large yaw angles, the horizontal force in Eqn 7 also needs to be updated from the direct measurement by plugging in the yaw adjusted magnitude of lift and drag. Finally, we calculated aerodynamic power using the full 3D aerodynamic force vector dotted with the 3D representative velocity vector of the wing, as summarized in Fig. 3D.

### Sensitivity analysis of aerodynamic power

Near stroke reversal, the aerodynamic power is sensitive to the directions of the lift and drag vectors (Fig. 2B,C). This sensitivity arises from the methodology we use to compute the lateral aerodynamic force, and here we show why the results are sensitive at stroke reversal and how we mitigate this issue.

Before mitigating the effects of computed aerodynamic power sensitivity, we first need to identify why it occurs near stroke reversal and quantify it. Sensitivity occurs when the denominator of Eqn 14 (*S̅*): is close to zero, causing the magnitude of the drag vector to spike (Fig. 2C). The magnitude of drag contributes directly to the aerodynamic power, making this particular denominator critical in understanding the overall sensitivity of the aerodynamic power calculation to wing orientation. In the absence of yaw (θ=0), the denominator is made up of the horizontal (subscript x) and vertical (subscript z) components of the unit lift vector (**L̂**) and the unit drag vector (**D̂**), or in other words, the 3D unit lift and drag vectors projected into the 2D sagittal plane of the dove (2D vectors **L̂**_{xz} and **D̂**_{xz} respectively). To build intuition for why lift and drag are sensitive at stroke reversal, consider the following cases when the denominator approaches zero: (1) The drag vector projected into the 2D sagittal plane (**D̂**_{xz}) is small, and (2) the lift and drag vectors projected into the 2D sagittal plane (**L̂**_{xz} and **D̂**_{xz}) either point in the same direction or the opposite direction. Case 1 occurs when the wing velocity, and thus drag, is pointed primarily in the lateral (y) direction, while case 2 occurs when either lift or drag is pointed primarily in the lateral direction. This explains why the highest aerodynamic power sensitivity is at stroke reversal when the wing is moving primarily in the lateral direction.

Now that we have a way of quantifying the level of sensitivity, we can mitigate its effects. In order to mitigate spikes in computed lateral force in sensitive regions, we smoothed the computed aerodynamic force in the lateral direction *F*_{wing,y} with weights *W* determined by the following formula:
(17)

We designed this formula, based on classical fuzzy logic (Bansod et al., 2005; Zadeh, 1988), to decrease the weight on points for which *S̅* (Eqn 16) approaches zero. When the magnitude of the denominator drops below a value of *c*_{0}=0.15, the weight equals zero because the result is too sensitive to be used. Conversely, if the magnitude of the denominator is greater than *c*_{1}=0.35, all values are weighted equally with a value of one. The *c*_{0} and *c*_{1} parameters were tuned by hand, and the logarithm serves to smooth the transition at *c*_{1}. We plotted the computed lift and drag magnitudes in Fig. 2B and the original versus the smoothed aerodynamic forces in Fig. 2D to show this approach is practical and useful to resolve mathematical singularities.

## RESULTS AND DISCUSSION

Here, we present how the aerodynamic force platform (AFP) can be used as an ergometer by computing the aerodynamic power required for flapping flight directly as the vector dot product of force and velocity. By using this method, we are able to compute aerodynamic power (Fig. 3) based on first principles: we directly measured the velocity vector, we directly measured two of three components of force, and the third component of force is a logical outcome based on first principles. We demonstrate the utility of this new method for the second stroke after takeoff of a dove by computing the x, y and z components of aerodynamic power (Fig. 3A). We found that the majority of the aerodynamic power is produced during the downstroke (Fig. 3C) and that its magnitude is primarily determined by the product of force and velocity in the vertical (z) direction. Further, the largest contribution of power in the forward (x) direction occurs during the latter half of the downstroke, and the largest contribution of power in the lateral (y) direction occurs at the beginning of the downstroke.

To determine the influence of the methods used to calculate the representative wing velocity vector on aerodynamic power, we compared the quarter-chord based velocity calculation with the two other methods that require less information (Fig. 3B,C). While the quarter-chord based calculation requires the leading and trailing edges of the wing to be precisely identified, the chord-averaged velocity method can be used when the leading and trailing edges are difficult to determine (i.e. when the location of the quarter-chord point cannot be determined). Despite this simplification, both the quarter-chord and chord-averaged methods give similar time-resolved trends throughout most of the stroke. This makes sense because both methods leverage the full 3D reconstructed wing surface, suggesting that the chord-averaged calculation, which is more robust to error in the wing surface reconstruction, could be the best approach when surface data are noisy. The main difference occurs during mid-downstroke, when the chord-averaged method overestimates power (for this single flight). This can be attributed to the difference in the estimated location of the center of pressure (Fig. S1), which influences the aerodynamic power mid-downstroke. In contrast, the quarter-chord estimate closely matches the aerodynamic power for a center of pressure location computed based on angle of attack (Dickson et al., 2008). When we compare the calculated power for the well-resolved (quarter-chord method) versus the minimalistic wing (shoulder and ninth primary) kinematics, we again found that the methods give overall similar time-resolved results. The power calculated using the minimalistic method deviates during both the mid-downstroke and the down-to-upstroke transition, on average overestimating the power (for this single flight). This is due to differences in both the direction (Fig. 1E) and the magnitude (Fig. 1D) of the velocity vector compared with the quarter-chord method. Additionally, the minimalistic model gives less well-resolved results at the end of the downstroke, because the denominator is closer to zero over a longer time window (Fig. 2C). Still, we saw that even coarse wing kinematics are sufficient to provide informative power measurements using the AFP as an ergometer. Future computational fluid dynamics or robot flapping wing studies could provide further insight into how the velocity calculation can be improved beyond what we present here.

Finally, we compared our method to calculate the aerodynamic power based on *in vivo* kinematics and aerodynamic forces with two commonly used aerodynamic models: the quasi-steady model based on lift-drag polars of spinning animal wings (Dickinson et al., 1999; Dickson et al., 2008) and the actuator disk model (Ellington, 1984b; Pennycuick, 2008) for predicting the aerodynamic power of flapping animal wings (for details and results see Fig. 3B,C and Fig. S2). Both models are based on *in vivo* kinematics, but whereas the quasi-steady model predicts the stroke-resolved aerodynamic power, the actuator disk model only predicts the stroke-averaged power. Comparing the predictive performance of these models for this single flight, we found that the quasi-steady model identifies a peak in aerodynamic power during mid-downstroke, whereas our direct method identifies a dip. This dip in aerodynamic power can be attributed to a dip in the angle of attack during mid-downstroke (Fig. S3A), during which the lift peaks and the drag dips (Fig. 2B). Additionally, during the upstroke, the quasi-steady model overpredicts power by a significant margin (Fig. 3B,C). The stroke-averaged power prediction is, however, much closer and only 22% higher than our method. The actuator disk model, on the other hand, underestimated the stroke-averaged power prediction by 28% compared with our method. This shows that both stroke-averaged model predictions are reasonable, and of these two models, the actuator disk model is easiest to implement. Furthermore, because the time-resolved aerodynamic power prediction by the quasi-steady model (based on lift-drag polars for a fruit fly shaped wing) is erroneous during mid-downstroke and mid-upstroke, it is of limited use to analyze the time-resolved flight of doves. Even with more accurate lift-drag polars for dove wings, the aerodynamic power discrepancy mid-downstroke cannot be properly addressed. This is because the angle of attack of spinning wing polars does not have a mid-downstroke dip (Crandell and Tobalske, 2011; Dickinson et al., 1999; Dickson et al., 2008; Kruyt et al., 2014; Lentink and Dickinson, 2009a; Usherwood, 2005, 2010), as we measured *in vivo* (Fig. S3A). Similarly, the mid-upstroke peak predicted by the quasi-steady model cannot be prevented with more accurate polars (Crandell and Tobalske, 2011; Dickinson et al., 1999; Dickson et al., 2008; Kruyt et al., 2014; Lentink and Dickinson, 2009b; Usherwood, 2005; Usherwood, 2010). Including additional terms in our quasi-steady model, such as added mass effects and rotational lift (Dickinson et al., 1999; Dickson et al., 2008), did not resolve the discrepancy either. This is because added mass effects, which depend on wing acceleration, and rotational lift effects, which depend on angle of attack rotation, both peak at stroke reversal and thus have minimal effect mid-stroke (Figs S2, S3B). In contrast, the stroke-averaged model predictions are simpler and more reliable and thus preferable as estimates. Consequently, we found the following key advantages of our method: it requires the smallest number of assumptions, it provides stroke-resolved aerodynamic power, and it is based primarily on *in vivo* data and first-principle equations (except for the center of pressure model, which is corroborated from robot flapping wing experiments). These advantages are key for reconstructing wingbeat-resolved ‘aerodynamic workloops’ of flight muscles. However, whenever quick estimates are needed for aerodynamic power, the actuator disk model may provide a particularly simple estimate. Finally, our method not only provides the wingbeat-resolved aerodynamic power, it simultaneously provides the wingbeat-resolved wing kinematics and aerodynamic forces. All these parameters are informative for evaluating the *in vivo* external work performed by flight muscles in the context of the overall flight mechanics of the animal, which no model can currently offer with comparable accuracy and precision. For an optimal implementation to evaluate the aerodynamic performance of flying animals, we recommend following the AFP design guidelines in Hightower et al. (2017). The current AFP volume is a smaller design, meaning that the aerodynamic power exerted by the animal is lower than during free flight in the atmosphere as a result of wall effects (Hightower et al., 2017). However, when the main goal of the study is to evaluate muscle function using the AFP as an ergometer, our method succeeds in accurately measuring the aerodynamic power produced by the bird as it flies in the AFP, and this smaller volume design has the benefit of providing better time resolution.

In conclusion, our mathematical derivations based on first principles and experimental results for a slow flying dove show how the AFP, together with a simple kinematics analysis, can be used as an ergometer. This more direct measurement of aerodynamic power enables new investigations into how flying animals generate aerodynamic power by beating their wings *in vivo*. The aerodynamic power calculation algorithm is generally reliable and only requires numerical stabilization during stroke reversal (Fig. 2C), therefore this approach can be used to compute the time-resolved aerodynamic power generated by the flight muscles (which we will report elsewhere for more individuals). Further, the analysis presented here shows the promise of future applications of aerodynamic and hydrodynamic force platforms (Lentink, 2018) as ergometers for a range of locomotor behaviors, including flying, swimming and locomotion over a water surface. Finally, the mathematical framework we developed is equally useful for 3D aerodynamic force and power analyses based on flow field measurements using particle image velocimetry for aerial and aquatic locomotion.

## APPENDIX 1

### Measuring the surface of the wing in 3D

We automatically measured the 3D surface of the flying dove by using a structured-light system (Fig. 1A). Subsequently, by integrating manually clicked marker and tracker points, we fit 3D surfaces to the body, tail and left wing of the dove (Fig. 1B).

The 3D surface data were measured using an automated structured-light system designed to image rapidly deforming objects. We used five synchronized high-speed cameras (Phantom Miro: three M310, one R-311, one LC-310; Vision Research, Wayne, NJ, USA) and high-speed projectors (DLP^{®} LightCrafter™ E4500MKII™; EKB Technologies, Bat-Yam, Israel) operating at 1000 frames s^{−1}. All the cameras and projectors were simultaneously calibrated using a method we previously developed (Deetjen and Lentink, 2018). We then used the structured-light method we previously developed (Deetjen et al., 2017) to reconstruct the 3D shape of the dove. Changes in the stripe detection and stripe matching portions of the structured-light method were made to improve the 3D reconstruction robustness. We attached color filters (72.0 mm diameter custom filters; two blue: low-pass filter with 475 nm cutoff, two green: band-pass filter centered around 540 nm with 80 nm full width at half maximum, one red: high-pass filter with 600 nm cutoff; Andover Corporation, Salem, NH, USA) to each camera corresponding to the projecting color of its paired projector. Since there were only three different colors, some of the camera-projector pair colors overlapped. For these overlapping color pairs, we reduced the frame rate of each projector to 500 frames s^{−1} and alternated frames to distinguish between the patterns projected by the two projectors projecting in the same color.

In order to relate the 3D surface measurements of the dove to the force directions measured by the AFP and the direction of gravity, we identified a world reference frame (Fig. 1A). The direction of gravity was determined by tracking a falling sphere, using custom software developed in MATLAB R2017b (MathWorks, Inc., Natick, MA, USA). Next, the forward direction of the dove was determined by imaging a carbon fiber rod resting on the centers of both the takeoff and landing perch. Lastly, the perch positions were determined based on minimizing the reprojected error from visible points on each perch.

After acquiring 3D surface data, we fitted surfaces to the body, tail, and left wing of the dove (Fig. 1B). We assumed bilateral symmetry for the right wing because we focused our cameras and projectors on the position of the left wing. To begin, data points were manually categorized as body, tail, left wing or outlier points. This was achieved by first aligning the body in each frame using a weighted iterative closest point (ICP) algorithm (Kjer and Wilm, 2010) applied to the kinematic points on the body of the dove. Next, manually placed ellipsoids, cones and other 3D shapes were used to group the points that corresponded to each category. After categorizing the points, we fitted surfaces to the body, tail and left wing. The body surface, which we assumed did not change shape temporally, was fit using the Gridfit toolbox in MATLAB (https://www.mathworks.com/matlabcentral/fileexchange/8998-surface-fitting-using-gridfit), which was modified to fit a closed surface by using spherical coordinates. The tail surface was approximated as a flat plate with edges defined by a circular sector with a fixed 3D center and radius. The 3D data were fitted to this surface with the tail elevation, twist and spread angles varying smoothly in time. Finally, the left wing surface was fitted using the Gridfit toolbox modified to also vary smoothly in time. This fit was applied after the data for each frame was transformed to a ‘wing reference frame’, which we defined by using kinematics points on the wing. The edges of the wing were found by reprojecting the 3D surface back onto the five camera images to determine which locations were on the wing. The edges between frames were smoothed using morphological image processing and a modified version of the Gridfit toolbox to fit a closed edge using polar coordinates varying smoothly in time.

## APPENDIX 2

### Derivation of representative wing velocity for aerodynamic power

In order to calculate the time-resolved aerodynamic power a bird needs to exert in flight, we need to integrate measured body and wing kinematics (Fig. 1) together with aerodynamic forces (Fig. 2). The total aerodynamic power for both wings *P*_{aero} is the dot product of the wing force vector and the representative wing velocity vector times two if we assume symmetry between wings:
(A 1)

We measure two of the three components of the wing force directly, but we need to derive the representative wing velocity vector from first principles. To start, we can approximate Eqn A1 as a summation because there are contributions toward total power distributed across the entire wing: (A 2)

The summation is over *n* blade elements, which start at the shoulder and are spaced 1 mm apart. We expand the aerodynamic force of the *j*^{th} blade element which is made up of the lift vector and the drag vector:
(A 3)where ρ is air density, *S _{j}* is surface area,

*C*

_{L,j}and

*C*

_{D,j}are the lift and drag coefficients, and

**L̂**

_{j}and

**D̂**

_{j}are the lift and drag unit vectors. By definition, lift is perpendicular to velocity and drag points in the opposite direction as velocity , so Eqns A2,A3 together can be reduced to: (A 4)

Additionally, the same expansion of the aerodynamic force into lift and drag components can be applied directly to Eqn A1 to begin isolating the representative wing velocity: (A 5)

If we assume that the directions of wing lift and drag remain constant along the entire span of the wing, (**L̂**_{j} = **L̂** and **D̂**_{j} = **D̂** for all *j*), then we can reduce:
(A 6)

This simplification is possible because the *j*^{th} lift unit vector dotted with the representative wing velocity equals zero and the *j*^{th} drag unit vector dotted with the representative wing velocity equals the negative of the representative wing velocity. By combining Eqns A4 and A6 and solving for , we find that:
(A 7)Finally, if we assume that the drag coefficient is spanwise invariant (*C*_{D,j} = *C*_{D} for all *j*), and we recognize that velocity points in the opposite direction as drag, we can simplify:
(A 8)

We use this representative wing velocity vector to compute aerodynamic power based on this first-principles derivation.

## APPENDIX 3

### Derivation of body yaw angle effect on aerodynamic power

Thus far, we have assumed that the bird's flight is not only bilaterally symmetric, but also straight from one perch to the other: aligned with the world reference frame defined by the AFP. Here, we derive the needed adjustment for the case where the bird's flight direction is yawed with respect to the world reference frame. In this analysis we still make the reasonable assumption that the bird flies bilaterally symmetric along its flight path in an approximately straight line.

Because we still assume that the bird's flight is symmetric, the aerodynamic power produced by each wing is equivalent to each other. Hence, Eqn 1 remains valid, and we only need to calculate the representative wing velocity and aerodynamic force on the left wing. Furthermore, because the computation of the left wing velocity vector (Eqn 2) is based solely on kinematics, it is unaffected by the yaw angle. Similarly, the computation of the left wing drag and lift unit vectors are solely based on kinematics. Therefore, the key variables we need to compute in order to adjust the aerodynamic force (Eqn 6) for body yaw offset, are the adjusted magnitudes of drag, *D*_{A}, and lift, *L*_{A}.

In order to compute the adjusted magnitudes of drag and lift, we retrace the steps in Eqns 4–13 for both wings separately. The aerodynamic force from each wing adds up to the total force as: (A 9)where the subscripts L and R refer to the left and right wings respectively, and the x and z components of are computed in Eqns 4,5 using the AFP measurements. This can be expanded according to Eqn 6: (A 10)where the magnitudes of the drag and lift are the same for each wing for flights along an approximately straight line. From Eqns 10 and 11, we see that in order to expand further, we need to solve for the velocity vector of the right wing.

To solve for the representative wing velocity of the right wing, we utilize the bilateral symmetry of the wings in the bird body reference frame. We define the rotation matrix from the world reference frame to the bird body reference frame **R**_{ z} as:
(A 11)where the angle θ is the yaw angle of the bird. For this particular flight, θ = 6.8 deg. Next, to simplify the notation, we represent the x, y and z components of the left wing velocity as:
(A 12)where the velocity vector is in the world reference frame. To find the components of the right wing, we first transform the left wing velocity into the symmetric bird body reference frame, represented by a tilde:
(A 13)

Next, from the symmetry assumption, we can find the velocity vector in the right wing: (A 14)and transform back into the world reference frame: (A 15) (A 16)

Now that we have the representative wing velocity for each wing, we use these vectors to compute the aerodynamic force on each wing. Similar to the velocity vector, we derive the drag and lift vectors for the right wing using the left wing. Since the drag vector points in the opposite direction as the velocity vector, we can use the same transformation: (A 17) (A 18)

Computing the lift vector requires the radial vector, **r̂**_{span}, which has the same symmetry properties as the velocity vector:
(A 19)
(A 20)

For the left wing, the expanded lift vector equals:
(A 21)
(A 22)where the variables *A*, *B* and *C* are used to simplify the notation. For the right wing, the expanded lift vector equals:
(A 23)
(A 24)Substituting Eqns A17, A18, A22, and A24 into Eqn A10 results in:
(A 25)which follows the same format as Eqn 6. Hence, we can solve for the adjusted drag and lift magnitudes in the same manner as Eqns 12–13:
(A 26)
(A 27)
(A 28)

As expected, when the yaw angle is zero, the adjusted drag and lift magnitudes equal the non-adjusted values in Eqns 12,13. Further, at body yaw angles of ±90 deg, there is a singularity which should be avoided by sufficiently aligning the world reference frame with the bird's flight path as we do in our study. We used these body-yaw adjusted magnitudes for drag and lift to compute the aerodynamic force on each wing, which we subsequently used to compute the total aerodynamic power. The small influence of this body yaw adjustment on our aerodynamic power calculations is shown in Fig. S4.

**Acknowledgements**

We thank Bret Tobalske for helping us train the dove to fly from perch to perch and helping us perform the experiment.

## FOOTNOTES

**Competing interests**The authors declare no competing or financial interests.

**Author contributions**Conceptualization: M.E.D., D.D.C., D.L.; Methodology: M.E.D., D.D.C., D.L.; Software: M.E.D.; Validation: M.E.D.; Formal analysis: M.E.D., D.D.C., D.L.; Investigation: M.E.D., D.D.C., D.L.; Resources: D.L.; Data curation: M.E.D.; Writing - original draft: M.E.D., D.D.C., D.L.; Writing - review & editing: M.E.D., D.D.C., D.L.; Visualization: M.E.D., D.D.C., D.L.; Supervision: D.L.; Project administration: D.L.; Funding acquisition: D.L.

**Funding**This research was supported by the National Science Foundation (NSF) Graduate Research Fellowship under grant no. DGE-114747 to M.D., the Air Force Office of Scientific Research (AFOSR) under grant no. FA9550-16-1-0182, and NSF CAREER Award 1552419 to D.L., Stanford University Graduate Fellowship to M.D. and D.D.C., and National Defense Science and Engineering Graduate Fellowship to D.D.C.

**Data availability**Raw data from this study are available from the Dryad Digital Repository (Deetjen et al., 2020): 15dv41ntr.

**Supplementary information**Supplementary information available online at http://jeb.biologists.org/lookup/doi/10.1242/jeb.220475.supplemental

- Received June 21, 2019.
- Accepted March 26, 2020.

- © 2020. Published by The Company of Biologists Ltd