R/mainrunfunctions.R
In charvey23/AvInertia: Calculate the Inertial Properties of a Flying Bird
Defines functions
store_data
massprop_birdwing
massprop_restbody
combine_inertialprop
Documented in
combine_inertialprop
massprop_birdwing
massprop_restbody
store_data
# Script containing the overarching function that calculates the moment of
# inertia and CG for an entire bird
# ------------------------------------------------------------------------------
# ------------------ Mass Properties - Combine all body components -------------
# ------------------------------------------------------------------------------
#' Combine body and wing inertial components.
#'
#' @description Combines data exported from massprop_restbody
#' and massprop_birdwing.
#'
#' @param curr_torsotail_data {Output dataframe from massprop_restbody.}
#' @param left_wing_data {Output dataframe from massprop_birdwing.}
#' @param right_wing_data {Output dataframe from massprop_birdwing.}
#' @param dat_id_curr Dataframe related to the current bird wing ID
#' info that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string.}
#' \item{BirdID}{Bird ID code as a string.}
#' \item{TestID}{Test ID code as a string.}
#' \item{frameID}{Video frame ID code as a string.}
#' }
#' @param dat_bird_curr Dataframe related to the current bird wing that must
#' include the following columns:
#' \itemize{
#' \item{total_bird_mass}{Mass of full bird for the current wing (kg).}
#' }
#' @param symmetric {Logical indicating if the input wings are symmetric or not.
#' If True than left_wing_data = right_wing_data.}
#'
#' @return a dataframe containing all of the inertial properties for each wing
#' component and the full bird about it's center of gravity and the vehicle
#' reference point (VRP)
#'
#' @examples
#' # refer to the vignette
#' library(AvInertia)
#'
#' # load data
#' data(dat_id_curr, package = "AvInertia")
#' data(dat_bird_curr, package = "AvInertia")
#' data(dat_feat_curr, package = "AvInertia")
#' data(dat_bone_curr, package = "AvInertia")
#' data(dat_mat, package = "AvInertia")
#' data(clean_pts, package = "AvInertia")
#'
#' # 1. Determine the center of gravity of the bird's torso (including the legs)
#' dat_torsotail_out = massprop_restbody(dat_id_curr, dat_bird_curr)
#' # 2. Calculate the inertia of the flight feathers about the tip of the calamus
#' feather_inertia <- compute_feat_inertia(dat_mat, dat_feat_curr, dat_bird_curr)
#' # 3. Determine the center of gravity of one of the bird's wings
#' dat_wing_out = massprop_birdwing(dat_id_curr, dat_bird_curr,
#' dat_bone_curr, dat_feat_curr, dat_mat, clean_pts,
#' feather_inertia, plot_var = 0)
#' # Visualize the center of gravity of each wing component in the x and y axis
#' dat_wing_out = massprop_birdwing(dat_id_curr, dat_bird_curr,
#' dat_bone_curr, dat_feat_curr, dat_mat, clean_pts,
#' feather_inertia, plot_var = "yx")
#' # or the y and z axis
#' dat_wing_out = massprop_birdwing(dat_id_curr, dat_bird_curr,
#' dat_bone_curr, dat_feat_curr, dat_mat, clean_pts,
#' feather_inertia, plot_var = "yz")
#' # 4. Combine all data and obtain the center of gravity, moment of inertia
#' # and principal axes of the bird
#' curr_full_bird = combine_inertialprop(dat_torsotail_out,dat_wing_out,
#' dat_wing_out, dat_id_curr, dat_bird_curr, symmetric=TRUE)
#'
#' @export
combine_inertialprop <- function(curr_torsotail_data,left_wing_data,
right_wing_data, dat_id_curr, dat_bird_curr, symmetric){
# --------------------- Initialize variables -----------------------
# pre-define storage matrices
mass_properties = as.data.frame(matrix(0, nrow = 0, ncol = 7)) # overall data
colnames(mass_properties) = c("species","BirdID","TestID","FrameID",
"component","object","value")
# Compute the full bird results
fullbird = list()
fullbird$I = matrix(0, nrow = 3, ncol = 3)
fullbird$CG = matrix(0, nrow = 3, ncol = 1)
# Set variables to NULL to avoid having to define as a global variable as
# CRAN can't see a binding for feather within the dataframe dat_feat_curr
component=NULL
object=NULL
# --------------------- Combine results -----------------------
# --- Mass ---
fullbird$m = sum(subset(curr_torsotail_data, object == "m")$value,
subset(left_wing_data, object == "m" &
component == "wing")$value,
subset(right_wing_data, object == "m" &
component == "wing")$value)
# --- Moment of Inertia tensor --- **Origin is about the VRP**
# Diagonal elements
fullbird$I[1,1] = sum(subset(curr_torsotail_data, object == "Ixx")$value) +
subset(left_wing_data, object == "Ixx" &
component == "wing")$value +
subset(right_wing_data, object == "Ixx" &
component == "wing")$value
fullbird$I[2,2] = sum(subset(curr_torsotail_data, object == "Iyy")$value) +
subset(left_wing_data, object == "Iyy" &
component == "wing")$value +
subset(right_wing_data, object == "Iyy" &
component == "wing")$value
fullbird$I[3,3] = sum(subset(curr_torsotail_data, object == "Izz")$value) +
subset(left_wing_data, object == "Izz" &
component == "wing")$value +
subset(right_wing_data, object == "Izz" &
component == "wing")$value
# ------ Off-diagonal elements ----
# only compute xz because the symmetry of the wing will cancel out xy and yz
fullbird$I[1,3] = sum(subset(curr_torsotail_data, object == "Ixz")$value) +
subset(left_wing_data, object == "Ixz" &
component == "wing")$value +
subset(right_wing_data, object == "Ixz" &
component == "wing")$value
fullbird$I[3,1] = fullbird$I[1,3]
# --- Center of gravity vector ---
# include neck contribution if it was calculated seperate from the head
if (length(subset(curr_torsotail_data, object == "m" &
component == "neck")$value) == 0) {
fullbird$CG[1] = (
subset(curr_torsotail_data, object == "CGx" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGx" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGx" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
fullbird$CG[3] = (
subset(curr_torsotail_data, object == "CGz" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGz" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGz" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
} else {
fullbird$CG[1] = (
subset(curr_torsotail_data, object == "CGx" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "neck")$value * subset(curr_torsotail_data, object == "m" &
component == "neck")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGx" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGx" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
fullbird$CG[3] = (
subset(curr_torsotail_data, object == "CGz" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "neck")$value * subset(curr_torsotail_data, object == "m" &
component == "neck")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGz" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGz" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
}
if (!symmetric){
fullbird$I[2,3] = sum(subset(curr_torsotail_data, object == "Iyz")$value) +
subset(left_wing_data, object == "Iyz" &
component == "wing")$value +
subset(right_wing_data, object == "Iyz" &
component == "wing")$value
fullbird$I[3,2] = fullbird$I[2,3]
fullbird$I[1,2] = sum(subset(curr_torsotail_data, object == "Ixy")$value) +
subset(left_wing_data, object == "Ixy" &
component == "wing")$value +
subset(right_wing_data, object == "Ixy" &
component == "wing")$value
fullbird$I[2,1] = fullbird$I[1,2]
if (length(subset(curr_torsotail_data, object == "m" &
component == "neck")$value) == 0) {
fullbird$CG[2] = (
subset(curr_torsotail_data, object == "CGy" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGy" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGy" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
} else{
fullbird$CG[2] = (
subset(curr_torsotail_data, object == "CGy" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "neck")$value * subset(curr_torsotail_data, object == "m" &
component == "neck")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGy" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGy" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
}
}
# Save the error between the measured bird mass and the final output mass
err_mass = fullbird$m - dat_bird_curr$total_bird_mass
dat_err = data.frame(species = dat_id_curr$species,
BirdID = dat_id_curr$BirdID,
TestID = dat_id_curr$TestID,
FrameID = dat_id_curr$FrameID,
component = "full", object = "m_err",
value = err_mass)
# ------------------------- Save the full data -----------------------------
# origin about the VRP
curr_full_bird_vrp = store_data(dat_id_curr,fullbird,
mass_properties,"full_VRP")
# ---------- Adjust the final moment of inertia tensor to be about the CG ----
fullbird$I = parallelaxis(fullbird$I,-fullbird$CG,fullbird$m,"A")
# store data so far
curr_full_bird = store_data(dat_id_curr,fullbird,
mass_properties,"full")
# calculate the principal axes
pri_axes = eigen(fullbird$I)
# --- saves the principal axes -----
pri_axes = eigen(fullbird$I)
new_row1 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "maj_axis_x",
value = pri_axes$vectors[1,1])
new_row2 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "maj_axis_y",
value = pri_axes$vectors[2,1])
new_row3 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "maj_axis_z",
value = pri_axes$vectors[3,1])
new_row4 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_axis_x",
value = pri_axes$vectors[1,2])
new_row5 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_axis_y",
value = pri_axes$vectors[2,2])
new_row6 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_axis_z",
value = pri_axes$vectors[3,2])
new_row7 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_axis_x",
value = pri_axes$vectors[1,3])
new_row8 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_axis_y",
value = pri_axes$vectors[2,3])
new_row9 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_axis_z",
value = pri_axes$vectors[3,3])
new_row10 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "max_eign",
value = pri_axes$values[1])
new_row11 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_eign",
value = pri_axes$values[2])
new_row12 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_eign",
value = pri_axes$values[3])
curr_full_bird = rbind(curr_full_bird,curr_full_bird_vrp[1:6,],dat_err,
new_row1,new_row2,new_row3,new_row4,new_row5,
new_row6,new_row7,new_row8,new_row9,new_row10,
new_row11,new_row12)
return(curr_full_bird)
}
# -------------------- Mass Properties - Body less wings -----------------------
#' Calculate the center of gravity and moment of inertia for the head, neck,
#' torso and tail.
#'
#' Function that reads in anatomical data and returns the moment of inertia
#' tensor and center of gravity for the head, neck, tail and torso.
#'
#' @param dat_wingID_curr Dataframe related to the current bird wing ID
#' info that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string.}
#' \item{BirdID}{Bird ID code as a string.}
#' \item{TestID}{Test ID code as a string.}
#' \item{frameID}{Video frame ID code as a string.}
#' }
#'
#' @param dat_bird_curr Dataframe related to the current bird wing that must
#' include the following columns:
#' \itemize{
#' \item{extend_neck}{Logical input defining whether the neck should be
#' modeled as extended or not.}
#' \item{head_length}{Length of the head from base to tip (m).}
#' \item{head_mass}{Mass of the head (kg).}
#' \item{head_height}{Height of the head at the base (m).}
#' \item{neck_mass}{Mass of the neck (kg).}
#' \item{neck_width}{OPTIONAL - Average width of the stretched neck (m).}
#' \item{neck_length}{OPTIONAL - Length of the stretched neck (m).}
#' \item{torsotail_length}{Length from the beginning of the torso to the tip of the tail (m).}
#' \item{torsotail_mass}{Mass of the torso and tail (kg).}
#' \item{tail_length}{Length of the tail (m).}
#' \item{tail_mass}{Mass of the tail (kg).}
#' \item{tail_width}{Average width of the furled tail (m).}
#' \item{right_leg_mass}{Mass of the right leg (kg).}
#' \item{left_leg_mass}{Mass of the left leg (kg).}
#' \item{body_width_max}{Maximum width of the torso (m).}
#' \item{body_width_at_leg_insert}{Width of the body at the point
#' where the legs are inserted (m).}
#' \item{x_loc_of_body_max}{x coordinate from the VRP in the
#' full bird frame of reference of the maximum body width (m).}
#' \item{x_loc_leg_insertion}{x coordinate from the VRP in the
#' full bird frame of reference of the leg insertion location (m).}
#' \item{x_loc_TorsotailCoG}{x coordinate from the VRP in the
#' full bird frame of reference of the center of gravity of the torso and tail (m).}
#' \item{z_loc_TorsotailCoG}{x coordinate from the VRP in the
#' full bird frame of reference of the the center of gravity of the torso and tail (m).}
#' }
#'
#' @return Function returns a dataframe that includes the moment of inertia and
#' center of gravity of head, neck, torso and tail.
#'
#' @inherit combine_inertialprop examples
#'
#' @export
#'
massprop_restbody <- function(dat_wingID_curr, dat_bird_curr){
# --------------------- Initialize variables -----------------------
# pre-define storage matrices
mass_properties = as.data.frame(matrix(0, nrow = 0, ncol = 7)) # overall data
colnames(mass_properties) = c("species","BirdID","TestID","FrameID",
"component","object","value")
# -------------------------------------------------------------
# ----------------- Head and neck data ------------------------
# -------------------------------------------------------------
neck_start = c(0,0,0)
if (dat_bird_curr$extend_neck){
neck_end = c(dat_bird_curr$neck_length,0,0)
head_end = neck_end + c(dat_bird_curr$head_length,0,0)
# Calculate the effects of the head
head = massprop_head(dat_bird_curr$head_mass,
0.5*dat_bird_curr$head_height,
dat_bird_curr$head_length,neck_end,head_end)
# Calculate the effects of the neck
neck = massprop_neck(dat_bird_curr$neck_mass,
0.5*dat_bird_curr$neck_width,
dat_bird_curr$neck_length,neck_start,neck_end)
} else{
head_end = c(dat_bird_curr$head_length,0,0)
# Calculate the effects of the head + neck
head = massprop_head((dat_bird_curr$head_mass+dat_bird_curr$neck_mass),
0.5*dat_bird_curr$head_height,
dat_bird_curr$head_length,neck_start,head_end)
}
# -------------------------------------------------------------
# ------------------- Torso and tail data -------------------------
# -------------------------------------------------------------
tail_start = c(-(dat_bird_curr$torsotail_length-dat_bird_curr$tail_length),0,0)
tail_end = c(-dat_bird_curr$torsotail_length,0,0)
m_legs = dat_bird_curr$right_leg_mass + dat_bird_curr$left_leg_mass
tail = massprop_tail(dat_bird_curr$tail_mass,
dat_bird_curr$tail_length,
dat_bird_curr$tail_width,tail_start,tail_end)
if(abs(tail$CG[1]) > abs(tail_end[1]) | abs(tail$CG[1]) < abs(tail_start[1])){
warning("Tail CG is incorrect")
}
# adjust the COG from the torso tail to be torso only
m_torso = dat_bird_curr$torsotail_mass - dat_bird_curr$tail_mass
l_torso = dat_bird_curr$torsotail_length - dat_bird_curr$tail_length
CG_x_torso = ((dat_bird_curr$torsotail_mass*dat_bird_curr$x_loc_TorsotailCoG) -
(dat_bird_curr$tail_mass*tail$CG[1]))/m_torso
CG_z_torso = ((dat_bird_curr$torsotail_mass*dat_bird_curr$z_loc_TorsotailCoG) -
(dat_bird_curr$tail_mass*tail$CG[3]))/m_torso
# CAUTION: INGOING CG_x must be positive as it calculates a total distance
# along the known axis
torso = massprop_torso(m_torso, m_legs,
dat_bird_curr$body_width_max,
dat_bird_curr$body_height_max,
dat_bird_curr$x_loc_of_body_max,
dat_bird_curr$body_width_at_leg_insert,
dat_bird_curr$x_loc_leg_insertion,
l_torso, abs(CG_x_torso),
CG_z_torso, neck_start, tail_start)
# ----------------------------------------------------
# ----------------- Save Data ------------------------
# ----------------------------------------------------
# ---- Head ----
mass_properties = store_data(dat_wingID_curr,head,mass_properties,"head")
# ---- Neck ----
if (dat_bird_curr$extend_neck){
mass_properties = store_data(dat_wingID_curr,neck,mass_properties,"neck")
}
# ---- Torso/Legs ----
mass_properties = store_data(dat_wingID_curr,torso,mass_properties,"torso")
# ---- Tail ----
mass_properties = store_data(dat_wingID_curr,tail,mass_properties,"tail")
return(mass_properties)
}
# ----------------- Mass Properties - Halfspan bird wing ----------------------
#' Calculate the center of gravity and moment of inertia for a
#' halfspan wing.
#'
#' Function that reads in anatomical data and returns the moment of inertia
#' tensor and center of gravity of a wing one side of the bird.
#'
#' @param dat_wingID_curr Dataframe related to the current bird wing ID info
#' that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string.}
#' \item{BirdID}{Bird ID code as a string.}
#' \item{TestID}{Test ID code as a string.}
#' \item{frameID}{Video frame ID code as a string.}
#' }
#'
#' @param dat_bird_curr Dataframe related to the current bird wing that must
#' include the following columns:
#' \itemize{
#' \item{total_bird_mass}{Mass of full bird for the current wing (kg).}
#' \item{wing_mass}{Mass of one wing, should be the current wing (kg).}
#' \item{barb_radius}{Radius of feather barb for current species (m).}
#' \item{barb_distance}{Distance between feather barbs for current species (m).}
#' \item{brachial_muscle_mass}{Mass of all muscles in the brachial region
#' of the wing (kg).}
#' \item{antebrachial_muscle_mass}{Mass of all muscles in the antebrachial region
#' of the wing (kg).}
#' \item{manus_muscle_mass}{Mass of all muscles in the manus region
#' of the wing (kg).}
#' \item{all_skin_coverts_mass}{Mass of all skin and covert feathers (kg).}
#' \item{tertiary_mass}{Mass of tertiary feathers (kg).}
#' }
#'
#' @param dat_bone_curr Dataframe (6 row x 5 column) related to the current bird
#' wing bones that must include the following columns:
#' \itemize{
#' \item{bone}{Bone ID code. Must include:
#' "Humerus","Ulna","Radius","Carpometacarpus","Ulnare" and "Radiale".}
#' \item{bone_mass}{Mass of bone in the same row as the appropriate
#' bone ID code (kg).}
#' \item{bone_len}{Length of bone in the same row as the appropriate
#' bone ID code (m).}
#' \item{bone_out_rad}{Outer radius of bone in the same row as the appropriate
#' bone ID code (m).}
#' \item{bone_in_rad}{Inner radius of bone in the same row as the appropriate
#' bone ID code (m).}
#' }
#'
#' @param dat_feat_curr Dataframe related to the current bird wing feathers
#' input as a dataframe with the following structure:
#' \itemize{
#' \item{feather}{Feather ID code. Must be in standard format i.e.
#' 1st primary is "P1", third secondary is "S3", etc.
#' Alula feathers should be grouped and named "alula".}
#' \item{m_f}{Mass of feather in the same row as the
#' appropriate feather ID code (kg).}
#' \item{l_cal}{Length of calamus in the same row as the
#' appropriate feather ID code (m).}
#' \item{l_vane}{Length of rachis/vane in the same row as the
#' appropriate feather ID code (m).}
#' \item{w_cal}{Width (diameter) of calamus in the same row as the
#' appropriate feather ID code (m).}
#' \item{w_vp}{Width of proximal vane (average value) in the same row as the
#' appropriate feather ID code (m).}
#' \item{w_vd}{Width of distal vane (average value) in the same row as the
#' appropriate feather ID code (m).}
#' \item{vane_angle}{Interior angle between the rachis and calamus (degrees).}
#' }
#' NOTE: Alula feathers will be treated as point mass so only the mass of the
#' feathers is required. Other columns can be left blank.
#'
#' @param dat_mat_curr Dataframe related to the current species input as a
#' dataframe with the following structure:
#' \itemize{
#' \item{material}{Material information. Must include the following:
#' "Bone","Skin","Muscle","Cortex", "Medullary"}
#' \item{density}{Density of each material (kg/m^3).}
#' }
#'
#' @param clean_pts A data frame of the key positions of the bird as follows:
#' \itemize{
#' \item{pt1x, pt1y, pt1z}{Point on the shoulder joint (m).}
#' \item{pt2x, pt1y, pt2z}{Point on the elbow joint (m).}
#' \item{pt3x, pt3y, pt3z}{Point on the wrist joint (m).}
#' \item{pt4x, pt4y, pt4z}{Point on the end of carpometacarpus (m).}
#' \item{pt6x, pt6y, pt6z}{Point on the leading edge of the wing in front of the
#' wrist joint (m).}
#' \item{pt8x, pt8y, pt8z}{Point on tip of most distal primary (m).}
#' \item{pt9x, pt9y, pt9z}{Point on the tip of the last primary to model as if
#' it is on the end of the carpometacarpus (m).}
#' \item{pt10x, pt10y, pt10z}{Point on tip of last primary to model as if
#' it was distributed along the carpometacarpus (m).}
#' \item{pt11x, pt11y, pt11z}{Point on tip of most proximal feather (m).}
#' \item{pt12x, pt12y, pt12z}{Point on exterior shoulder position
#' (wing root leading edge) (m).}
#' }
#'
#' @param feather_inertia A list with one entry per flight feather. Each primary feather includes the following variables:
#' \itemize{
#' \item{I_pri}{a 3x3 matrix representing the moment of inertia about each feather calamus tip (kg-m^2).}
#' \item{CG_pri}{a 1x3 vector (x,y,z) representing the center of gravity of the primary feather (m).}
#' \item{m_pri}{a double representing the mass of the primary feather (kg).}
#' }
#' Each secondary feather includes the following variables:
#' \itemize{
#' \item{I_sec}{a 3x3 matrix representing the moment of inertia about each feather calamus tip (kg-m^2).}
#' \item{CG_sec}{a 1x3 vector (x,y,z) representing the center of gravity of the primary feather (m).}
#' \item{m_sec}{a double representing the mass of the primary feather (kg).}
#' }
#'
#' @param plot_var A string that defines the x-axis and y-axis of the output plot.
#' Can either equal "yx" or "yz".
#'
#' @section CAUTION:
#' All points must all have the vehicle reference point (VRP) as their
#' origin and the vehicle major axes as their frame of reference. This
#' is normally selected so that the VRP is in line with the body center
#' of gravity. Ensure the axes used represent a right-handed axis system.
#'
#' @author Christina Harvey
#'
#' @return Function returns a dataframe that includes the moment of inertia and
#' center of gravity of one wing about the VRP in the VRP frame and that of each
#' major anatomical group i.e. skin, feathers, bones, muscles.
#'
#' @inherit combine_inertialprop examples
#'
#' @export
#'
massprop_birdwing <- function(dat_wingID_curr, dat_bird_curr, dat_bone_curr,
dat_feat_curr, dat_mat_curr, clean_pts,
feather_inertia, plot_var){
# --------------------- Initialize variables -----------------------
mass_properties = as.data.frame(matrix(0, nrow = 0, ncol = 7)) # overall data
column_names = c("species","BirdID","TestID","FrameID",
"component","object","value")
colnames(mass_properties) = column_names
# Set variables to NULL to avoid having to define as a global variable as
# CRAN can't see a binding for feather within the dataframe dat_feat_curr
feather=NULL
bone=NULL
mass_properties_bone = mass_properties # specific bone data
mass_properties_muscle = mass_properties # specific muscle data
mass_properties_skin = mass_properties # specific skin data
mass_properties_feathers = mass_properties # individual feather mass data
# define incoming points
Pt1 = clean_pts[1,] # shoulder
Pt2 = clean_pts[2,] # elbow
Pt3 = clean_pts[3,] # wrist
Pt4 = clean_pts[4,] # end of carpometacarpus (cmc)
Pt6 = clean_pts[5,] # leading edge in front of wrist
Pt8 = clean_pts[6,] # tip of most distal primary
Pt9 = clean_pts[7,] # tip of last primary as if on the end of the cmc
Pt10 = clean_pts[8,] # tip of S1
Pt11 = clean_pts[9,] # tip of last secondary feather at wing root
Pt12 = clean_pts[10,] # leading edge of wing root
# --------------------------------------------------
# --------------- Bone Data ------------------------
# --------------------------------------------------
rho_bone = dat_mat_curr$density[which(dat_mat_curr$material == "Bone")]
dat_bone_hum = subset(dat_bone_curr, bone == "Humerus")
dat_bone_uln = subset(dat_bone_curr, bone == "Ulna")
dat_bone_rad = subset(dat_bone_curr, bone == "Radius")
dat_bone_car = subset(dat_bone_curr, bone == "Carpometacarpus")
hum = massprop_bones(dat_bone_hum$bone_mass,dat_bone_hum$bone_len,
dat_bone_hum$bone_out_rad,dat_bone_hum$bone_in_rad,
rho_bone, Pt1, Pt2)
ulna = massprop_bones(dat_bone_uln$bone_mass,dat_bone_uln$bone_len,
dat_bone_uln$bone_out_rad,dat_bone_uln$bone_in_rad,
rho_bone, Pt2, Pt3)
radius = massprop_bones(dat_bone_rad$bone_mass,dat_bone_rad$bone_len,
dat_bone_rad$bone_out_rad,dat_bone_rad$bone_in_rad,
rho_bone, Pt2, Pt3)
car = massprop_bones(dat_bone_car$bone_mass,dat_bone_car$bone_len,
dat_bone_car$bone_out_rad,dat_bone_car$bone_in_rad,
rho_bone, Pt3, Pt4)
wristbone = massprop_pm((subset(dat_bone_curr,
bone == "Ulnare")$bone_mass +
subset(dat_bone_curr,
bone == "Radiale")$bone_mass), Pt3)
# --- All Bones ---
prop_bone = list()
# simply addition as long as about the same origin in the same frame of
# reference
#(Frame of reference: VRP | Origin: VRP)
prop_bone$I = hum$I + ulna$I + radius$I + car$I + wristbone$I
# weighted average of the individual center of mass
#(Frame of reference: VRP | Origin: VRP)
prop_bone$CG = (
dat_bone_hum$bone_mass * hum$CG +
dat_bone_uln$bone_mass * ulna$CG +
dat_bone_rad$bone_mass * radius$CG +
dat_bone_car$bone_mass * car$CG +
(
subset(dat_bone_curr,
bone == "Ulnare")$bone_mass +
subset(dat_bone_curr,
bone == "Radiale")$bone_mass
) * wristbone$CG
) / sum(dat_bone_curr$bone_mass)
prop_bone$m = sum(dat_bone_curr$bone_mass)
# ----------------------------------------------------
# --------------- Muscle Data ------------------------
# ----------------------------------------------------
rho_muscle = dat_mat_curr$density[which(dat_mat_curr$material == "Muscle")]
mass_muscles = c()
mass_muscles[1] = dat_bird_curr$brachial_muscle_mass
mass_muscles[2] = dat_bird_curr$antebrachial_muscle_mass
mass_muscles[3] = dat_bird_curr$manus_muscle_mass
brach = massprop_muscles(mass_muscles[1],rho_muscle,dat_bone_hum$bone_len,
Pt1,Pt2)
abrach = massprop_muscles(mass_muscles[2],rho_muscle,dat_bone_uln$bone_len,
Pt2,Pt3)
manus = massprop_muscles(mass_muscles[3],rho_muscle,dat_bone_car$bone_len,
Pt3,Pt4)
# --- All Muscles ---
prop_muscles = list()
# simply addition as long as about same origin in the same frame of reference
# (Frame of reference: VRP | Origin: VRP)
prop_muscles$I = brach$I + abrach$I + manus$I
# weighted average of the individual center of mass
# (Frame of reference: VRP | Origin: VRP)
prop_muscles$CG = (mass_muscles[1]*brach$CG + mass_muscles[2]*abrach$CG +
mass_muscles[3]*manus$CG)/sum(mass_muscles)
prop_muscles$m = sum(mass_muscles)
# -------------------------------------------------------
# ----------------- Feather Data ------------------------
# -------------------------------------------------------
# density information
rho_cor = dat_mat_curr$density[which(dat_mat_curr$material == "Cortex")]
rho_med = dat_mat_curr$density[which(dat_mat_curr$material == "Medullary")]
# separate out primaries and secondaries
primaries = dat_feat_curr[grep("P",dat_feat_curr$feather),]
secondaries = dat_feat_curr[grep("S",dat_feat_curr$feather),]
no_sec = length(secondaries$feather)
no_pri = length(primaries$feather)
#pre-define storage matrices
res_pri = list()
res_pri$I = array(dim = c(3,3,no_pri))
res_pri$CG = array(dim = c(no_pri,3))
res_pri$CGm = array(dim = c(no_pri,3))
res_pri$m = array(dim = c(no_pri))
res_sec = list()
res_sec$I = array(dim = c(3,3,no_sec))
res_sec$CG = array(dim = c(no_sec,3))
res_sec$CGm = array(dim = c(no_sec,3))
res_sec$m = array(dim = c(no_sec))
# determine the orientation and normal of each feather
feather_info = orient_feather(no_pri,no_sec,Pt1,Pt2,Pt3,Pt4,
Pt8,Pt9,Pt10,Pt11,Pt12)
# --------------------------- Primaries --------------------------------------
# P1 -> P10
for (i in 1:no_pri){
# Adjust MOI and CG to the current orientation
tmp = structural2VRP_feat(feather_inertia$m_pri[i],
feather_inertia$I_pri[,,i],
feather_inertia$CG_pri[i,],
feather_info$loc_start[i,],
feather_info$loc_end[i,],
feather_info$normal[i,])
# Save MOI, CG and CG*mass
res_pri$I[,,i] = tmp$I
res_pri$CG[i,] = tmp$CG
res_pri$CGm[i,] = tmp$CG*tmp$m
res_pri$m[i] = tmp$m
}
# --------------------------Secondaries -------------------------------------
# S1 -> last secondary
for (i in 1:no_sec){
# Adjust MOI and CG to the current orientation
tmp = structural2VRP_feat(feather_inertia$m_sec[i],
feather_inertia$I_sec[,,i],
feather_inertia$CG_sec[i,],
feather_info$loc_start[i+no_pri,],
feather_info$loc_end[i+no_pri,],
feather_info$normal[i+no_pri,])
# Save MOI, CG and CG*mass
res_sec$I[,,i] = tmp$I
res_sec$CG[i,] = tmp$CG
res_sec$CGm[i,] = tmp$CG*tmp$m
res_sec$m[i] = tmp$m
}
# -------------------------------- Alula -------------------------------------
m_alula = subset(dat_feat_curr,feather == "alula")$m_f
alula = massprop_pm(m_alula, Pt6)
# ------------------------------ Tertiaries ----------------------------------
# position where the teritaries likely encounter the body
edge_tert = c(Pt11[1],min(Pt12[2],Pt11[2]),Pt12[3])
prop_tertiary1 = massprop_skin(0.5*dat_bird_curr$tertiary_mass,
rho_cor,rbind(Pt12,edge_tert,Pt2))
prop_tertiary2 = massprop_skin(0.5*dat_bird_curr$tertiary_mass,
rho_cor,rbind(Pt11,Pt2,edge_tert))
# --- All Feathers ---
prop_feathers = list()
I_feathers_pri = rbind(rowSums(res_pri$I[,1,]),
rowSums(res_pri$I[,2,]),
rowSums(res_pri$I[,3,]))
I_feathers_sec = rbind(rowSums(res_sec$I[,1,]),
rowSums(res_sec$I[,2,]),
rowSums(res_sec$I[,3,]))
prop_feathers$I = I_feathers_pri + I_feathers_sec + alula$I +
prop_tertiary1$I + prop_tertiary2$I
prop_feathers$m = sum(dat_feat_curr$m_f) + dat_bird_curr$tertiary_mass
prop_feathers$CG = (
colSums(res_pri$CGm) + colSums(res_sec$CGm) +
alula$CG * m_alula +
prop_tertiary1$CG * 0.5 * dat_bird_curr$tertiary_mass +
prop_tertiary2$CG * 0.5 * dat_bird_curr$tertiary_mass
) / prop_feathers$m
# ---------------------------------------------------------
# --------------- Skin/Covert Data ------------------------
# ---------------------------------------------------------
rho_skin = dat_mat_curr$density[which(dat_mat_curr$material == "Skin")]
mass_skin = dat_bird_curr$wing_mass - prop_feathers$m - prop_bone$m - prop_muscles$m
prop_skin = massprop_skin(mass_skin,rho_skin,rbind(Pt12,Pt6,Pt2))
# ----------------------------------------------------
# ----------------- Save Data ------------------------
# ----------------------------------------------------
# add data to bone specific data frame
mass_properties_bone = store_data(dat_wingID_curr,hum,
mass_properties_bone,"humerus")
mass_properties_bone = store_data(dat_wingID_curr,ulna,
mass_properties_bone,"ulna")
mass_properties_bone = store_data(dat_wingID_curr,radius,
mass_properties_bone,"radius")
mass_properties_bone = store_data(dat_wingID_curr,car,
mass_properties_bone,"carpo")
mass_properties_bone = store_data(dat_wingID_curr,wristbone,
mass_properties_bone,"wristbones")
# add data to muscle specific data frame
mass_properties_muscle = store_data(dat_wingID_curr,brach,
mass_properties_muscle,"brach")
mass_properties_muscle = store_data(dat_wingID_curr,abrach,
mass_properties_muscle,"abrach")
mass_properties_muscle = store_data(dat_wingID_curr,manus,
mass_properties_muscle,"manus")
# add data to skin specific data frame
mass_properties_skin = store_data(dat_wingID_curr,prop_skin,
mass_properties_skin,"skin_prop")
# add data to feather specific data frame
# --- Primaries ---
for (i in 1:no_pri){
feather_name = paste("P",i,sep = "")
curr_res_pri = list()
curr_res_pri$I = res_pri$I[,,i]
curr_res_pri$CG = res_pri$CG[i,]
curr_res_pri$m = res_pri$m[i]
mass_properties_feathers = store_data(dat_wingID_curr,curr_res_pri,
mass_properties_feathers,feather_name)
}
# --- Secondaries ---
for (i in 1:no_sec){
feather_name = paste("S",i,sep = "")
curr_res_sec = list()
curr_res_sec$I = res_sec$I[,,i]
curr_res_sec$CG = res_sec$CG[i,]
curr_res_sec$m = res_sec$m[i]
mass_properties_feathers = store_data(dat_wingID_curr,curr_res_sec,
mass_properties_feathers,feather_name)
}
# save all combined data groups to the master list
mass_properties = store_data(dat_wingID_curr,prop_bone,
mass_properties,"bones")
mass_properties = store_data(dat_wingID_curr,prop_muscles,
mass_properties,"muscles")
mass_properties = store_data(dat_wingID_curr,prop_skin,
mass_properties,"skin")
mass_properties = store_data(dat_wingID_curr,prop_feathers,
mass_properties,"feathers")
# save all wing data
prop_wing = list()
prop_wing$I = prop_bone$I + prop_muscles$I + prop_skin$I + prop_feathers$I
prop_wing$m = prop_bone$m + prop_muscles$m + prop_skin$m + prop_feathers$m
prop_wing$CG = (prop_bone$CG*prop_bone$m + prop_muscles$CG*prop_muscles$m +
prop_skin$CG*prop_skin$m +
prop_feathers$CG*prop_feathers$m)/prop_wing$m
mass_properties = store_data(dat_wingID_curr,prop_wing,mass_properties,"wing")
# Plot to verify correct outputs
if (plot_var != 0){
CGplot = plot_CGloc(clean_pts,mass_properties,mass_properties_skin,
mass_properties_bone,mass_properties_feathers,
mass_properties_muscle, prop_tertiary1,
prop_tertiary2, plot_var)
}
return(mass_properties)
}
# -------------------- Store Data -------------------------------
#' Store data from the inertia calculations in long format
#'
#' Function to store moment of inertia tensor and center of gravity vector
#' components in long format
#'
#' @param dat_wingID_curr Dataframe related to the current bird wing ID
#' info that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string}
#' \item{BirdID}{Bird ID code as a string}
#' \item{TestID}{Test ID code as a string}
#' \item{frameID}{Video frame ID code as a string}
#' }
#'
#' @param dat_mass Dataframe containing the new MOI and CG data to add
#' to mass_properties as new rows. Must include:
#' \itemize{
#' \item{I}{Moment of inertia tensor (kg-m^2)}
#' \item{CG}{Center of gravity with three location components (m)}
#' }
#'
#' @param mass_properties Dataframe containing any previously saved data.
#' Must have the following columns: "species","BirdID","TestID","FrameID",
#' "prop_type","component","value".
#'
#' @param name Name of the component for which the moment of inertia and
#' center of gravity were computed.
#'
#' @return This function returns mass_properties as an updated dataframe
#' with a new row corresponding to the dat_mass information
#'
#' @inherit combine_inertialprop examples
#'
#' @export
store_data <- function(dat_wingID_curr,dat_mass,mass_properties,name){
prop_type_list = c("Ixx","Iyy","Izz","Ixy","Iyz","Ixz","CGx","CGy","CGz")
prop_type_ind = cbind(c(1,2,3,1,2,1),c(1,2,3,2,3,3))
species = dat_wingID_curr$species
BirdID = dat_wingID_curr$BirdID
testID = dat_wingID_curr$TestID
frameID = dat_wingID_curr$FrameID
# Moment of inertia tensor
for (i in 1:6){
# saves the name and value of the tensor component
new_row = data.frame(species = species,
BirdID = BirdID,
TestID = testID,
FrameID = frameID,
component = name,
object = prop_type_list[i],
value = dat_mass$I[prop_type_ind[i,1],
prop_type_ind[i,2]])
mass_properties = rbind(mass_properties,new_row)
}
# Center of gravity
for (i in 1:3){
# saves the name and value of the CG component
new_row = data.frame(species = species,BirdID = BirdID,
TestID = testID,
FrameID = frameID,
component = name,
object = prop_type_list[6+i],
value = dat_mass$CG[i])
mass_properties = rbind(mass_properties,new_row)
}
# Mass
# saves the name and value of the mass
new_row = data.frame(species = species,
BirdID = BirdID,
TestID = testID,
FrameID = frameID,
component = name,
object = "m",
value = dat_mass$m)
mass_properties = rbind(mass_properties,new_row)
return(mass_properties)
}
charvey23/AvInertia documentation built on April 14, 2022, 1:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions store_data massprop_birdwing massprop_restbody combine_inertialprop
Documented in combine_inertialprop massprop_birdwing massprop_restbody store_data
# Script containing the overarching function that calculates the moment of
# inertia and CG for an entire bird
# ------------------------------------------------------------------------------
# ------------------ Mass Properties - Combine all body components -------------
# ------------------------------------------------------------------------------
#' Combine body and wing inertial components.
#'
#' @description Combines data exported from massprop_restbody
#' and massprop_birdwing.
#'
#' @param curr_torsotail_data {Output dataframe from massprop_restbody.}
#' @param left_wing_data {Output dataframe from massprop_birdwing.}
#' @param right_wing_data {Output dataframe from massprop_birdwing.}
#' @param dat_id_curr Dataframe related to the current bird wing ID
#' info that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string.}
#' \item{BirdID}{Bird ID code as a string.}
#' \item{TestID}{Test ID code as a string.}
#' \item{frameID}{Video frame ID code as a string.}
#' }
#' @param dat_bird_curr Dataframe related to the current bird wing that must
#' include the following columns:
#' \itemize{
#' \item{total_bird_mass}{Mass of full bird for the current wing (kg).}
#' }
#' @param symmetric {Logical indicating if the input wings are symmetric or not.
#' If True than left_wing_data = right_wing_data.}
#'
#' @return a dataframe containing all of the inertial properties for each wing
#' component and the full bird about it's center of gravity and the vehicle
#' reference point (VRP)
#'
#' @examples
#' # refer to the vignette
#' library(AvInertia)
#'
#' # load data
#' data(dat_id_curr, package = "AvInertia")
#' data(dat_bird_curr, package = "AvInertia")
#' data(dat_feat_curr, package = "AvInertia")
#' data(dat_bone_curr, package = "AvInertia")
#' data(dat_mat, package = "AvInertia")
#' data(clean_pts, package = "AvInertia")
#'
#' # 1. Determine the center of gravity of the bird's torso (including the legs)
#' dat_torsotail_out = massprop_restbody(dat_id_curr, dat_bird_curr)
#' # 2. Calculate the inertia of the flight feathers about the tip of the calamus
#' feather_inertia <- compute_feat_inertia(dat_mat, dat_feat_curr, dat_bird_curr)
#' # 3. Determine the center of gravity of one of the bird's wings
#' dat_wing_out = massprop_birdwing(dat_id_curr, dat_bird_curr,
#' dat_bone_curr, dat_feat_curr, dat_mat, clean_pts,
#' feather_inertia, plot_var = 0)
#' # Visualize the center of gravity of each wing component in the x and y axis
#' dat_wing_out = massprop_birdwing(dat_id_curr, dat_bird_curr,
#' dat_bone_curr, dat_feat_curr, dat_mat, clean_pts,
#' feather_inertia, plot_var = "yx")
#' # or the y and z axis
#' dat_wing_out = massprop_birdwing(dat_id_curr, dat_bird_curr,
#' dat_bone_curr, dat_feat_curr, dat_mat, clean_pts,
#' feather_inertia, plot_var = "yz")
#' # 4. Combine all data and obtain the center of gravity, moment of inertia
#' # and principal axes of the bird
#' curr_full_bird = combine_inertialprop(dat_torsotail_out,dat_wing_out,
#' dat_wing_out, dat_id_curr, dat_bird_curr, symmetric=TRUE)
#'
#' @export
combine_inertialprop <- function(curr_torsotail_data,left_wing_data,
right_wing_data, dat_id_curr, dat_bird_curr, symmetric){
# --------------------- Initialize variables -----------------------
# pre-define storage matrices
mass_properties = as.data.frame(matrix(0, nrow = 0, ncol = 7)) # overall data
colnames(mass_properties) = c("species","BirdID","TestID","FrameID",
"component","object","value")
# Compute the full bird results
fullbird = list()
fullbird$I = matrix(0, nrow = 3, ncol = 3)
fullbird$CG = matrix(0, nrow = 3, ncol = 1)
# Set variables to NULL to avoid having to define as a global variable as
# CRAN can't see a binding for feather within the dataframe dat_feat_curr
component=NULL
object=NULL
# --------------------- Combine results -----------------------
# --- Mass ---
fullbird$m = sum(subset(curr_torsotail_data, object == "m")$value,
subset(left_wing_data, object == "m" &
component == "wing")$value,
subset(right_wing_data, object == "m" &
component == "wing")$value)
# --- Moment of Inertia tensor --- **Origin is about the VRP**
# Diagonal elements
fullbird$I[1,1] = sum(subset(curr_torsotail_data, object == "Ixx")$value) +
subset(left_wing_data, object == "Ixx" &
component == "wing")$value +
subset(right_wing_data, object == "Ixx" &
component == "wing")$value
fullbird$I[2,2] = sum(subset(curr_torsotail_data, object == "Iyy")$value) +
subset(left_wing_data, object == "Iyy" &
component == "wing")$value +
subset(right_wing_data, object == "Iyy" &
component == "wing")$value
fullbird$I[3,3] = sum(subset(curr_torsotail_data, object == "Izz")$value) +
subset(left_wing_data, object == "Izz" &
component == "wing")$value +
subset(right_wing_data, object == "Izz" &
component == "wing")$value
# ------ Off-diagonal elements ----
# only compute xz because the symmetry of the wing will cancel out xy and yz
fullbird$I[1,3] = sum(subset(curr_torsotail_data, object == "Ixz")$value) +
subset(left_wing_data, object == "Ixz" &
component == "wing")$value +
subset(right_wing_data, object == "Ixz" &
component == "wing")$value
fullbird$I[3,1] = fullbird$I[1,3]
# --- Center of gravity vector ---
# include neck contribution if it was calculated seperate from the head
if (length(subset(curr_torsotail_data, object == "m" &
component == "neck")$value) == 0) {
fullbird$CG[1] = (
subset(curr_torsotail_data, object == "CGx" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGx" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGx" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
fullbird$CG[3] = (
subset(curr_torsotail_data, object == "CGz" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGz" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGz" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
} else {
fullbird$CG[1] = (
subset(curr_torsotail_data, object == "CGx" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "neck")$value * subset(curr_torsotail_data, object == "m" &
component == "neck")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGx" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGx" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGx" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
fullbird$CG[3] = (
subset(curr_torsotail_data, object == "CGz" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "neck")$value * subset(curr_torsotail_data, object == "m" &
component == "neck")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGz" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGz" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGz" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
}
if (!symmetric){
fullbird$I[2,3] = sum(subset(curr_torsotail_data, object == "Iyz")$value) +
subset(left_wing_data, object == "Iyz" &
component == "wing")$value +
subset(right_wing_data, object == "Iyz" &
component == "wing")$value
fullbird$I[3,2] = fullbird$I[2,3]
fullbird$I[1,2] = sum(subset(curr_torsotail_data, object == "Ixy")$value) +
subset(left_wing_data, object == "Ixy" &
component == "wing")$value +
subset(right_wing_data, object == "Ixy" &
component == "wing")$value
fullbird$I[2,1] = fullbird$I[1,2]
if (length(subset(curr_torsotail_data, object == "m" &
component == "neck")$value) == 0) {
fullbird$CG[2] = (
subset(curr_torsotail_data, object == "CGy" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGy" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGy" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
} else{
fullbird$CG[2] = (
subset(curr_torsotail_data, object == "CGy" &
component == "head")$value * subset(curr_torsotail_data, object == "m" &
component == "head")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "neck")$value * subset(curr_torsotail_data, object == "m" &
component == "neck")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "torso")$value * subset(curr_torsotail_data, object == "m" &
component == "torso")$value +
subset(curr_torsotail_data, object == "CGy" &
component == "tail")$value * subset(curr_torsotail_data, object == "m" &
component == "tail")$value +
subset(left_wing_data, object == "CGy" &
component == "wing")$value * subset(left_wing_data, object == "m" &
component == "wing")$value +
subset(right_wing_data, object == "CGy" &
component == "wing")$value * subset(right_wing_data, object == "m" &
component == "wing")$value
) / fullbird$m
}
}
# Save the error between the measured bird mass and the final output mass
err_mass = fullbird$m - dat_bird_curr$total_bird_mass
dat_err = data.frame(species = dat_id_curr$species,
BirdID = dat_id_curr$BirdID,
TestID = dat_id_curr$TestID,
FrameID = dat_id_curr$FrameID,
component = "full", object = "m_err",
value = err_mass)
# ------------------------- Save the full data -----------------------------
# origin about the VRP
curr_full_bird_vrp = store_data(dat_id_curr,fullbird,
mass_properties,"full_VRP")
# ---------- Adjust the final moment of inertia tensor to be about the CG ----
fullbird$I = parallelaxis(fullbird$I,-fullbird$CG,fullbird$m,"A")
# store data so far
curr_full_bird = store_data(dat_id_curr,fullbird,
mass_properties,"full")
# calculate the principal axes
pri_axes = eigen(fullbird$I)
# --- saves the principal axes -----
pri_axes = eigen(fullbird$I)
new_row1 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "maj_axis_x",
value = pri_axes$vectors[1,1])
new_row2 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "maj_axis_y",
value = pri_axes$vectors[2,1])
new_row3 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "maj_axis_z",
value = pri_axes$vectors[3,1])
new_row4 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_axis_x",
value = pri_axes$vectors[1,2])
new_row5 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_axis_y",
value = pri_axes$vectors[2,2])
new_row6 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_axis_z",
value = pri_axes$vectors[3,2])
new_row7 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_axis_x",
value = pri_axes$vectors[1,3])
new_row8 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_axis_y",
value = pri_axes$vectors[2,3])
new_row9 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_axis_z",
value = pri_axes$vectors[3,3])
new_row10 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "max_eign",
value = pri_axes$values[1])
new_row11 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "int_eign",
value = pri_axes$values[2])
new_row12 = data.frame(species = curr_full_bird$species[1],
BirdID = curr_full_bird$BirdID[1],
TestID = curr_full_bird$TestID[1],
FrameID = curr_full_bird$FrameID[1],
component = "full", object = "min_eign",
value = pri_axes$values[3])
curr_full_bird = rbind(curr_full_bird,curr_full_bird_vrp[1:6,],dat_err,
new_row1,new_row2,new_row3,new_row4,new_row5,
new_row6,new_row7,new_row8,new_row9,new_row10,
new_row11,new_row12)
return(curr_full_bird)
}
# -------------------- Mass Properties - Body less wings -----------------------
#' Calculate the center of gravity and moment of inertia for the head, neck,
#' torso and tail.
#'
#' Function that reads in anatomical data and returns the moment of inertia
#' tensor and center of gravity for the head, neck, tail and torso.
#'
#' @param dat_wingID_curr Dataframe related to the current bird wing ID
#' info that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string.}
#' \item{BirdID}{Bird ID code as a string.}
#' \item{TestID}{Test ID code as a string.}
#' \item{frameID}{Video frame ID code as a string.}
#' }
#'
#' @param dat_bird_curr Dataframe related to the current bird wing that must
#' include the following columns:
#' \itemize{
#' \item{extend_neck}{Logical input defining whether the neck should be
#' modeled as extended or not.}
#' \item{head_length}{Length of the head from base to tip (m).}
#' \item{head_mass}{Mass of the head (kg).}
#' \item{head_height}{Height of the head at the base (m).}
#' \item{neck_mass}{Mass of the neck (kg).}
#' \item{neck_width}{OPTIONAL - Average width of the stretched neck (m).}
#' \item{neck_length}{OPTIONAL - Length of the stretched neck (m).}
#' \item{torsotail_length}{Length from the beginning of the torso to the tip of the tail (m).}
#' \item{torsotail_mass}{Mass of the torso and tail (kg).}
#' \item{tail_length}{Length of the tail (m).}
#' \item{tail_mass}{Mass of the tail (kg).}
#' \item{tail_width}{Average width of the furled tail (m).}
#' \item{right_leg_mass}{Mass of the right leg (kg).}
#' \item{left_leg_mass}{Mass of the left leg (kg).}
#' \item{body_width_max}{Maximum width of the torso (m).}
#' \item{body_width_at_leg_insert}{Width of the body at the point
#' where the legs are inserted (m).}
#' \item{x_loc_of_body_max}{x coordinate from the VRP in the
#' full bird frame of reference of the maximum body width (m).}
#' \item{x_loc_leg_insertion}{x coordinate from the VRP in the
#' full bird frame of reference of the leg insertion location (m).}
#' \item{x_loc_TorsotailCoG}{x coordinate from the VRP in the
#' full bird frame of reference of the center of gravity of the torso and tail (m).}
#' \item{z_loc_TorsotailCoG}{x coordinate from the VRP in the
#' full bird frame of reference of the the center of gravity of the torso and tail (m).}
#' }
#'
#' @return Function returns a dataframe that includes the moment of inertia and
#' center of gravity of head, neck, torso and tail.
#'
#' @inherit combine_inertialprop examples
#'
#' @export
#'
massprop_restbody <- function(dat_wingID_curr, dat_bird_curr){
# --------------------- Initialize variables -----------------------
# pre-define storage matrices
mass_properties = as.data.frame(matrix(0, nrow = 0, ncol = 7)) # overall data
colnames(mass_properties) = c("species","BirdID","TestID","FrameID",
"component","object","value")
# -------------------------------------------------------------
# ----------------- Head and neck data ------------------------
# -------------------------------------------------------------
neck_start = c(0,0,0)
if (dat_bird_curr$extend_neck){
neck_end = c(dat_bird_curr$neck_length,0,0)
head_end = neck_end + c(dat_bird_curr$head_length,0,0)
# Calculate the effects of the head
head = massprop_head(dat_bird_curr$head_mass,
0.5*dat_bird_curr$head_height,
dat_bird_curr$head_length,neck_end,head_end)
# Calculate the effects of the neck
neck = massprop_neck(dat_bird_curr$neck_mass,
0.5*dat_bird_curr$neck_width,
dat_bird_curr$neck_length,neck_start,neck_end)
} else{
head_end = c(dat_bird_curr$head_length,0,0)
# Calculate the effects of the head + neck
head = massprop_head((dat_bird_curr$head_mass+dat_bird_curr$neck_mass),
0.5*dat_bird_curr$head_height,
dat_bird_curr$head_length,neck_start,head_end)
}
# -------------------------------------------------------------
# ------------------- Torso and tail data -------------------------
# -------------------------------------------------------------
tail_start = c(-(dat_bird_curr$torsotail_length-dat_bird_curr$tail_length),0,0)
tail_end = c(-dat_bird_curr$torsotail_length,0,0)
m_legs = dat_bird_curr$right_leg_mass + dat_bird_curr$left_leg_mass
tail = massprop_tail(dat_bird_curr$tail_mass,
dat_bird_curr$tail_length,
dat_bird_curr$tail_width,tail_start,tail_end)
if(abs(tail$CG[1]) > abs(tail_end[1]) | abs(tail$CG[1]) < abs(tail_start[1])){
warning("Tail CG is incorrect")
}
# adjust the COG from the torso tail to be torso only
m_torso = dat_bird_curr$torsotail_mass - dat_bird_curr$tail_mass
l_torso = dat_bird_curr$torsotail_length - dat_bird_curr$tail_length
CG_x_torso = ((dat_bird_curr$torsotail_mass*dat_bird_curr$x_loc_TorsotailCoG) -
(dat_bird_curr$tail_mass*tail$CG[1]))/m_torso
CG_z_torso = ((dat_bird_curr$torsotail_mass*dat_bird_curr$z_loc_TorsotailCoG) -
(dat_bird_curr$tail_mass*tail$CG[3]))/m_torso
# CAUTION: INGOING CG_x must be positive as it calculates a total distance
# along the known axis
torso = massprop_torso(m_torso, m_legs,
dat_bird_curr$body_width_max,
dat_bird_curr$body_height_max,
dat_bird_curr$x_loc_of_body_max,
dat_bird_curr$body_width_at_leg_insert,
dat_bird_curr$x_loc_leg_insertion,
l_torso, abs(CG_x_torso),
CG_z_torso, neck_start, tail_start)
# ----------------------------------------------------
# ----------------- Save Data ------------------------
# ----------------------------------------------------
# ---- Head ----
mass_properties = store_data(dat_wingID_curr,head,mass_properties,"head")
# ---- Neck ----
if (dat_bird_curr$extend_neck){
mass_properties = store_data(dat_wingID_curr,neck,mass_properties,"neck")
}
# ---- Torso/Legs ----
mass_properties = store_data(dat_wingID_curr,torso,mass_properties,"torso")
# ---- Tail ----
mass_properties = store_data(dat_wingID_curr,tail,mass_properties,"tail")
return(mass_properties)
}
# ----------------- Mass Properties - Halfspan bird wing ----------------------
#' Calculate the center of gravity and moment of inertia for a
#' halfspan wing.
#'
#' Function that reads in anatomical data and returns the moment of inertia
#' tensor and center of gravity of a wing one side of the bird.
#'
#' @param dat_wingID_curr Dataframe related to the current bird wing ID info
#' that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string.}
#' \item{BirdID}{Bird ID code as a string.}
#' \item{TestID}{Test ID code as a string.}
#' \item{frameID}{Video frame ID code as a string.}
#' }
#'
#' @param dat_bird_curr Dataframe related to the current bird wing that must
#' include the following columns:
#' \itemize{
#' \item{total_bird_mass}{Mass of full bird for the current wing (kg).}
#' \item{wing_mass}{Mass of one wing, should be the current wing (kg).}
#' \item{barb_radius}{Radius of feather barb for current species (m).}
#' \item{barb_distance}{Distance between feather barbs for current species (m).}
#' \item{brachial_muscle_mass}{Mass of all muscles in the brachial region
#' of the wing (kg).}
#' \item{antebrachial_muscle_mass}{Mass of all muscles in the antebrachial region
#' of the wing (kg).}
#' \item{manus_muscle_mass}{Mass of all muscles in the manus region
#' of the wing (kg).}
#' \item{all_skin_coverts_mass}{Mass of all skin and covert feathers (kg).}
#' \item{tertiary_mass}{Mass of tertiary feathers (kg).}
#' }
#'
#' @param dat_bone_curr Dataframe (6 row x 5 column) related to the current bird
#' wing bones that must include the following columns:
#' \itemize{
#' \item{bone}{Bone ID code. Must include:
#' "Humerus","Ulna","Radius","Carpometacarpus","Ulnare" and "Radiale".}
#' \item{bone_mass}{Mass of bone in the same row as the appropriate
#' bone ID code (kg).}
#' \item{bone_len}{Length of bone in the same row as the appropriate
#' bone ID code (m).}
#' \item{bone_out_rad}{Outer radius of bone in the same row as the appropriate
#' bone ID code (m).}
#' \item{bone_in_rad}{Inner radius of bone in the same row as the appropriate
#' bone ID code (m).}
#' }
#'
#' @param dat_feat_curr Dataframe related to the current bird wing feathers
#' input as a dataframe with the following structure:
#' \itemize{
#' \item{feather}{Feather ID code. Must be in standard format i.e.
#' 1st primary is "P1", third secondary is "S3", etc.
#' Alula feathers should be grouped and named "alula".}
#' \item{m_f}{Mass of feather in the same row as the
#' appropriate feather ID code (kg).}
#' \item{l_cal}{Length of calamus in the same row as the
#' appropriate feather ID code (m).}
#' \item{l_vane}{Length of rachis/vane in the same row as the
#' appropriate feather ID code (m).}
#' \item{w_cal}{Width (diameter) of calamus in the same row as the
#' appropriate feather ID code (m).}
#' \item{w_vp}{Width of proximal vane (average value) in the same row as the
#' appropriate feather ID code (m).}
#' \item{w_vd}{Width of distal vane (average value) in the same row as the
#' appropriate feather ID code (m).}
#' \item{vane_angle}{Interior angle between the rachis and calamus (degrees).}
#' }
#' NOTE: Alula feathers will be treated as point mass so only the mass of the
#' feathers is required. Other columns can be left blank.
#'
#' @param dat_mat_curr Dataframe related to the current species input as a
#' dataframe with the following structure:
#' \itemize{
#' \item{material}{Material information. Must include the following:
#' "Bone","Skin","Muscle","Cortex", "Medullary"}
#' \item{density}{Density of each material (kg/m^3).}
#' }
#'
#' @param clean_pts A data frame of the key positions of the bird as follows:
#' \itemize{
#' \item{pt1x, pt1y, pt1z}{Point on the shoulder joint (m).}
#' \item{pt2x, pt1y, pt2z}{Point on the elbow joint (m).}
#' \item{pt3x, pt3y, pt3z}{Point on the wrist joint (m).}
#' \item{pt4x, pt4y, pt4z}{Point on the end of carpometacarpus (m).}
#' \item{pt6x, pt6y, pt6z}{Point on the leading edge of the wing in front of the
#' wrist joint (m).}
#' \item{pt8x, pt8y, pt8z}{Point on tip of most distal primary (m).}
#' \item{pt9x, pt9y, pt9z}{Point on the tip of the last primary to model as if
#' it is on the end of the carpometacarpus (m).}
#' \item{pt10x, pt10y, pt10z}{Point on tip of last primary to model as if
#' it was distributed along the carpometacarpus (m).}
#' \item{pt11x, pt11y, pt11z}{Point on tip of most proximal feather (m).}
#' \item{pt12x, pt12y, pt12z}{Point on exterior shoulder position
#' (wing root leading edge) (m).}
#' }
#'
#' @param feather_inertia A list with one entry per flight feather. Each primary feather includes the following variables:
#' \itemize{
#' \item{I_pri}{a 3x3 matrix representing the moment of inertia about each feather calamus tip (kg-m^2).}
#' \item{CG_pri}{a 1x3 vector (x,y,z) representing the center of gravity of the primary feather (m).}
#' \item{m_pri}{a double representing the mass of the primary feather (kg).}
#' }
#' Each secondary feather includes the following variables:
#' \itemize{
#' \item{I_sec}{a 3x3 matrix representing the moment of inertia about each feather calamus tip (kg-m^2).}
#' \item{CG_sec}{a 1x3 vector (x,y,z) representing the center of gravity of the primary feather (m).}
#' \item{m_sec}{a double representing the mass of the primary feather (kg).}
#' }
#'
#' @param plot_var A string that defines the x-axis and y-axis of the output plot.
#' Can either equal "yx" or "yz".
#'
#' @section CAUTION:
#' All points must all have the vehicle reference point (VRP) as their
#' origin and the vehicle major axes as their frame of reference. This
#' is normally selected so that the VRP is in line with the body center
#' of gravity. Ensure the axes used represent a right-handed axis system.
#'
#' @author Christina Harvey
#'
#' @return Function returns a dataframe that includes the moment of inertia and
#' center of gravity of one wing about the VRP in the VRP frame and that of each
#' major anatomical group i.e. skin, feathers, bones, muscles.
#'
#' @inherit combine_inertialprop examples
#'
#' @export
#'
massprop_birdwing <- function(dat_wingID_curr, dat_bird_curr, dat_bone_curr,
dat_feat_curr, dat_mat_curr, clean_pts,
feather_inertia, plot_var){
# --------------------- Initialize variables -----------------------
mass_properties = as.data.frame(matrix(0, nrow = 0, ncol = 7)) # overall data
column_names = c("species","BirdID","TestID","FrameID",
"component","object","value")
colnames(mass_properties) = column_names
# Set variables to NULL to avoid having to define as a global variable as
# CRAN can't see a binding for feather within the dataframe dat_feat_curr
feather=NULL
bone=NULL
mass_properties_bone = mass_properties # specific bone data
mass_properties_muscle = mass_properties # specific muscle data
mass_properties_skin = mass_properties # specific skin data
mass_properties_feathers = mass_properties # individual feather mass data
# define incoming points
Pt1 = clean_pts[1,] # shoulder
Pt2 = clean_pts[2,] # elbow
Pt3 = clean_pts[3,] # wrist
Pt4 = clean_pts[4,] # end of carpometacarpus (cmc)
Pt6 = clean_pts[5,] # leading edge in front of wrist
Pt8 = clean_pts[6,] # tip of most distal primary
Pt9 = clean_pts[7,] # tip of last primary as if on the end of the cmc
Pt10 = clean_pts[8,] # tip of S1
Pt11 = clean_pts[9,] # tip of last secondary feather at wing root
Pt12 = clean_pts[10,] # leading edge of wing root
# --------------------------------------------------
# --------------- Bone Data ------------------------
# --------------------------------------------------
rho_bone = dat_mat_curr$density[which(dat_mat_curr$material == "Bone")]
dat_bone_hum = subset(dat_bone_curr, bone == "Humerus")
dat_bone_uln = subset(dat_bone_curr, bone == "Ulna")
dat_bone_rad = subset(dat_bone_curr, bone == "Radius")
dat_bone_car = subset(dat_bone_curr, bone == "Carpometacarpus")
hum = massprop_bones(dat_bone_hum$bone_mass,dat_bone_hum$bone_len,
dat_bone_hum$bone_out_rad,dat_bone_hum$bone_in_rad,
rho_bone, Pt1, Pt2)
ulna = massprop_bones(dat_bone_uln$bone_mass,dat_bone_uln$bone_len,
dat_bone_uln$bone_out_rad,dat_bone_uln$bone_in_rad,
rho_bone, Pt2, Pt3)
radius = massprop_bones(dat_bone_rad$bone_mass,dat_bone_rad$bone_len,
dat_bone_rad$bone_out_rad,dat_bone_rad$bone_in_rad,
rho_bone, Pt2, Pt3)
car = massprop_bones(dat_bone_car$bone_mass,dat_bone_car$bone_len,
dat_bone_car$bone_out_rad,dat_bone_car$bone_in_rad,
rho_bone, Pt3, Pt4)
wristbone = massprop_pm((subset(dat_bone_curr,
bone == "Ulnare")$bone_mass +
subset(dat_bone_curr,
bone == "Radiale")$bone_mass), Pt3)
# --- All Bones ---
prop_bone = list()
# simply addition as long as about the same origin in the same frame of
# reference
#(Frame of reference: VRP | Origin: VRP)
prop_bone$I = hum$I + ulna$I + radius$I + car$I + wristbone$I
# weighted average of the individual center of mass
#(Frame of reference: VRP | Origin: VRP)
prop_bone$CG = (
dat_bone_hum$bone_mass * hum$CG +
dat_bone_uln$bone_mass * ulna$CG +
dat_bone_rad$bone_mass * radius$CG +
dat_bone_car$bone_mass * car$CG +
(
subset(dat_bone_curr,
bone == "Ulnare")$bone_mass +
subset(dat_bone_curr,
bone == "Radiale")$bone_mass
) * wristbone$CG
) / sum(dat_bone_curr$bone_mass)
prop_bone$m = sum(dat_bone_curr$bone_mass)
# ----------------------------------------------------
# --------------- Muscle Data ------------------------
# ----------------------------------------------------
rho_muscle = dat_mat_curr$density[which(dat_mat_curr$material == "Muscle")]
mass_muscles = c()
mass_muscles[1] = dat_bird_curr$brachial_muscle_mass
mass_muscles[2] = dat_bird_curr$antebrachial_muscle_mass
mass_muscles[3] = dat_bird_curr$manus_muscle_mass
brach = massprop_muscles(mass_muscles[1],rho_muscle,dat_bone_hum$bone_len,
Pt1,Pt2)
abrach = massprop_muscles(mass_muscles[2],rho_muscle,dat_bone_uln$bone_len,
Pt2,Pt3)
manus = massprop_muscles(mass_muscles[3],rho_muscle,dat_bone_car$bone_len,
Pt3,Pt4)
# --- All Muscles ---
prop_muscles = list()
# simply addition as long as about same origin in the same frame of reference
# (Frame of reference: VRP | Origin: VRP)
prop_muscles$I = brach$I + abrach$I + manus$I
# weighted average of the individual center of mass
# (Frame of reference: VRP | Origin: VRP)
prop_muscles$CG = (mass_muscles[1]*brach$CG + mass_muscles[2]*abrach$CG +
mass_muscles[3]*manus$CG)/sum(mass_muscles)
prop_muscles$m = sum(mass_muscles)
# -------------------------------------------------------
# ----------------- Feather Data ------------------------
# -------------------------------------------------------
# density information
rho_cor = dat_mat_curr$density[which(dat_mat_curr$material == "Cortex")]
rho_med = dat_mat_curr$density[which(dat_mat_curr$material == "Medullary")]
# separate out primaries and secondaries
primaries = dat_feat_curr[grep("P",dat_feat_curr$feather),]
secondaries = dat_feat_curr[grep("S",dat_feat_curr$feather),]
no_sec = length(secondaries$feather)
no_pri = length(primaries$feather)
#pre-define storage matrices
res_pri = list()
res_pri$I = array(dim = c(3,3,no_pri))
res_pri$CG = array(dim = c(no_pri,3))
res_pri$CGm = array(dim = c(no_pri,3))
res_pri$m = array(dim = c(no_pri))
res_sec = list()
res_sec$I = array(dim = c(3,3,no_sec))
res_sec$CG = array(dim = c(no_sec,3))
res_sec$CGm = array(dim = c(no_sec,3))
res_sec$m = array(dim = c(no_sec))
# determine the orientation and normal of each feather
feather_info = orient_feather(no_pri,no_sec,Pt1,Pt2,Pt3,Pt4,
Pt8,Pt9,Pt10,Pt11,Pt12)
# --------------------------- Primaries --------------------------------------
# P1 -> P10
for (i in 1:no_pri){
# Adjust MOI and CG to the current orientation
tmp = structural2VRP_feat(feather_inertia$m_pri[i],
feather_inertia$I_pri[,,i],
feather_inertia$CG_pri[i,],
feather_info$loc_start[i,],
feather_info$loc_end[i,],
feather_info$normal[i,])
# Save MOI, CG and CG*mass
res_pri$I[,,i] = tmp$I
res_pri$CG[i,] = tmp$CG
res_pri$CGm[i,] = tmp$CG*tmp$m
res_pri$m[i] = tmp$m
}
# --------------------------Secondaries -------------------------------------
# S1 -> last secondary
for (i in 1:no_sec){
# Adjust MOI and CG to the current orientation
tmp = structural2VRP_feat(feather_inertia$m_sec[i],
feather_inertia$I_sec[,,i],
feather_inertia$CG_sec[i,],
feather_info$loc_start[i+no_pri,],
feather_info$loc_end[i+no_pri,],
feather_info$normal[i+no_pri,])
# Save MOI, CG and CG*mass
res_sec$I[,,i] = tmp$I
res_sec$CG[i,] = tmp$CG
res_sec$CGm[i,] = tmp$CG*tmp$m
res_sec$m[i] = tmp$m
}
# -------------------------------- Alula -------------------------------------
m_alula = subset(dat_feat_curr,feather == "alula")$m_f
alula = massprop_pm(m_alula, Pt6)
# ------------------------------ Tertiaries ----------------------------------
# position where the teritaries likely encounter the body
edge_tert = c(Pt11[1],min(Pt12[2],Pt11[2]),Pt12[3])
prop_tertiary1 = massprop_skin(0.5*dat_bird_curr$tertiary_mass,
rho_cor,rbind(Pt12,edge_tert,Pt2))
prop_tertiary2 = massprop_skin(0.5*dat_bird_curr$tertiary_mass,
rho_cor,rbind(Pt11,Pt2,edge_tert))
# --- All Feathers ---
prop_feathers = list()
I_feathers_pri = rbind(rowSums(res_pri$I[,1,]),
rowSums(res_pri$I[,2,]),
rowSums(res_pri$I[,3,]))
I_feathers_sec = rbind(rowSums(res_sec$I[,1,]),
rowSums(res_sec$I[,2,]),
rowSums(res_sec$I[,3,]))
prop_feathers$I = I_feathers_pri + I_feathers_sec + alula$I +
prop_tertiary1$I + prop_tertiary2$I
prop_feathers$m = sum(dat_feat_curr$m_f) + dat_bird_curr$tertiary_mass
prop_feathers$CG = (
colSums(res_pri$CGm) + colSums(res_sec$CGm) +
alula$CG * m_alula +
prop_tertiary1$CG * 0.5 * dat_bird_curr$tertiary_mass +
prop_tertiary2$CG * 0.5 * dat_bird_curr$tertiary_mass
) / prop_feathers$m
# ---------------------------------------------------------
# --------------- Skin/Covert Data ------------------------
# ---------------------------------------------------------
rho_skin = dat_mat_curr$density[which(dat_mat_curr$material == "Skin")]
mass_skin = dat_bird_curr$wing_mass - prop_feathers$m - prop_bone$m - prop_muscles$m
prop_skin = massprop_skin(mass_skin,rho_skin,rbind(Pt12,Pt6,Pt2))
# ----------------------------------------------------
# ----------------- Save Data ------------------------
# ----------------------------------------------------
# add data to bone specific data frame
mass_properties_bone = store_data(dat_wingID_curr,hum,
mass_properties_bone,"humerus")
mass_properties_bone = store_data(dat_wingID_curr,ulna,
mass_properties_bone,"ulna")
mass_properties_bone = store_data(dat_wingID_curr,radius,
mass_properties_bone,"radius")
mass_properties_bone = store_data(dat_wingID_curr,car,
mass_properties_bone,"carpo")
mass_properties_bone = store_data(dat_wingID_curr,wristbone,
mass_properties_bone,"wristbones")
# add data to muscle specific data frame
mass_properties_muscle = store_data(dat_wingID_curr,brach,
mass_properties_muscle,"brach")
mass_properties_muscle = store_data(dat_wingID_curr,abrach,
mass_properties_muscle,"abrach")
mass_properties_muscle = store_data(dat_wingID_curr,manus,
mass_properties_muscle,"manus")
# add data to skin specific data frame
mass_properties_skin = store_data(dat_wingID_curr,prop_skin,
mass_properties_skin,"skin_prop")
# add data to feather specific data frame
# --- Primaries ---
for (i in 1:no_pri){
feather_name = paste("P",i,sep = "")
curr_res_pri = list()
curr_res_pri$I = res_pri$I[,,i]
curr_res_pri$CG = res_pri$CG[i,]
curr_res_pri$m = res_pri$m[i]
mass_properties_feathers = store_data(dat_wingID_curr,curr_res_pri,
mass_properties_feathers,feather_name)
}
# --- Secondaries ---
for (i in 1:no_sec){
feather_name = paste("S",i,sep = "")
curr_res_sec = list()
curr_res_sec$I = res_sec$I[,,i]
curr_res_sec$CG = res_sec$CG[i,]
curr_res_sec$m = res_sec$m[i]
mass_properties_feathers = store_data(dat_wingID_curr,curr_res_sec,
mass_properties_feathers,feather_name)
}
# save all combined data groups to the master list
mass_properties = store_data(dat_wingID_curr,prop_bone,
mass_properties,"bones")
mass_properties = store_data(dat_wingID_curr,prop_muscles,
mass_properties,"muscles")
mass_properties = store_data(dat_wingID_curr,prop_skin,
mass_properties,"skin")
mass_properties = store_data(dat_wingID_curr,prop_feathers,
mass_properties,"feathers")
# save all wing data
prop_wing = list()
prop_wing$I = prop_bone$I + prop_muscles$I + prop_skin$I + prop_feathers$I
prop_wing$m = prop_bone$m + prop_muscles$m + prop_skin$m + prop_feathers$m
prop_wing$CG = (prop_bone$CG*prop_bone$m + prop_muscles$CG*prop_muscles$m +
prop_skin$CG*prop_skin$m +
prop_feathers$CG*prop_feathers$m)/prop_wing$m
mass_properties = store_data(dat_wingID_curr,prop_wing,mass_properties,"wing")
# Plot to verify correct outputs
if (plot_var != 0){
CGplot = plot_CGloc(clean_pts,mass_properties,mass_properties_skin,
mass_properties_bone,mass_properties_feathers,
mass_properties_muscle, prop_tertiary1,
prop_tertiary2, plot_var)
}
return(mass_properties)
}
# -------------------- Store Data -------------------------------
#' Store data from the inertia calculations in long format
#'
#' Function to store moment of inertia tensor and center of gravity vector
#' components in long format
#'
#' @param dat_wingID_curr Dataframe related to the current bird wing ID
#' info that must include the following columns:
#' \itemize{
#' \item{species}{Species ID code as a string}
#' \item{BirdID}{Bird ID code as a string}
#' \item{TestID}{Test ID code as a string}
#' \item{frameID}{Video frame ID code as a string}
#' }
#'
#' @param dat_mass Dataframe containing the new MOI and CG data to add
#' to mass_properties as new rows. Must include:
#' \itemize{
#' \item{I}{Moment of inertia tensor (kg-m^2)}
#' \item{CG}{Center of gravity with three location components (m)}
#' }
#'
#' @param mass_properties Dataframe containing any previously saved data.
#' Must have the following columns: "species","BirdID","TestID","FrameID",
#' "prop_type","component","value".
#'
#' @param name Name of the component for which the moment of inertia and
#' center of gravity were computed.
#'
#' @return This function returns mass_properties as an updated dataframe
#' with a new row corresponding to the dat_mass information
#'
#' @inherit combine_inertialprop examples
#'
#' @export
store_data <- function(dat_wingID_curr,dat_mass,mass_properties,name){
prop_type_list = c("Ixx","Iyy","Izz","Ixy","Iyz","Ixz","CGx","CGy","CGz")
prop_type_ind = cbind(c(1,2,3,1,2,1),c(1,2,3,2,3,3))
species = dat_wingID_curr$species
BirdID = dat_wingID_curr$BirdID
testID = dat_wingID_curr$TestID
frameID = dat_wingID_curr$FrameID
# Moment of inertia tensor
for (i in 1:6){
# saves the name and value of the tensor component
new_row = data.frame(species = species,
BirdID = BirdID,
TestID = testID,
FrameID = frameID,
component = name,
object = prop_type_list[i],
value = dat_mass$I[prop_type_ind[i,1],
prop_type_ind[i,2]])
mass_properties = rbind(mass_properties,new_row)
}
# Center of gravity
for (i in 1:3){
# saves the name and value of the CG component
new_row = data.frame(species = species,BirdID = BirdID,
TestID = testID,
FrameID = frameID,
component = name,
object = prop_type_list[6+i],
value = dat_mass$CG[i])
mass_properties = rbind(mass_properties,new_row)
}
# Mass
# saves the name and value of the mass
new_row = data.frame(species = species,
BirdID = BirdID,
TestID = testID,
FrameID = frameID,
component = name,
object = "m",
value = dat_mass$m)
mass_properties = rbind(mass_properties,new_row)
return(mass_properties)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.