# to run this code use ctrl f to find and replace all of the instances of the joint name e.g. Shoulder with another joint name e.g. Elbow
# will then have to manually correct for the columns, never figured out a clean way to do it

#this code uses the ggplot2 library in the plotting section of the loop functions
library(ggplot2)

#All of the loops follow the same layout

#remember to use setwd () to set the directory to the location of the monte_carlo outputs

#loop function for the complete data set
graph.loop <- function(j, df, ti, ti2)
{
  ##Load and Setup##
  {
    
    buname <- (paste0("Burst/Outputs/"))
    biname <- (paste0("Biped/Outputs/"))
    quname <- (paste0("Quad/Outputs/"))
    buname2 <- (paste0(buname,j,df))
    biname2 <- (paste0(biname,j,df))
    quname2 <- (paste0(quname,j,df))
    
    # read in raw data
    
    mma_bu   <- read.csv((paste0(buname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_bi   <- read.csv((paste0(biname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_qu   <- read.csv((paste0(quname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdbu <- read.csv((paste0(buname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdbi <- read.csv((paste0(biname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdqu <- read.csv((paste0(quname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    zero <- read.csv("zero.csv", header = TRUE, sep = ",", quote =
                       "\"",dec = ".", fill = TRUE, comment.char = "")
    
    
    #make mean+/- 2 sd
    mmap_bu  <- mma_bu
    mmam_bu  <- mma_bu
    mmap_bi  <- mma_bi
    mmam_bi  <- mma_bi
    mmap_qu  <- mma_qu
    mmam_qu  <- mma_qu
    
    {
      if("BI" %in% colnames(mma_bi)){
        mmap_bi$BI      <- mma_bi$BI + (2*mma_sdbi$BI)
        mmam_bi$BI      <- mma_bi$BI - (2*mma_sdbi$BI)
        mmap_bu$BI      <- mma_bu$BI + (2*mma_sdbu$BI)
        mmam_bu$BI      <- mma_bu$BI - (2*mma_sdbu$BI)
        mmap_qu$BI      <- mma_qu$BI + (2*mma_sdqu$BI)
        mmam_qu$BI      <- mma_qu$BI - (2*mma_sdqu$BI)
      }
      if("CB" %in% colnames(mma_bi)){
        mmap_bi$CB      <- mma_bi$CB + (2*mma_sdbi$CB)
        mmam_bi$CB      <- mma_bi$CB - (2*mma_sdbi$CB)
        mmap_bu$CB      <- mma_bu$CB + (2*mma_sdbu$CB)
        mmam_bu$CB      <- mma_bu$CB - (2*mma_sdbu$CB)
        mmap_qu$CB      <- mma_qu$CB + (2*mma_sdqu$CB)
        mmam_qu$CB      <- mma_qu$CB - (2*mma_sdqu$CB)               
      }
      if("BR" %in% colnames(mma_bi)){
        mmap_bi$BR      <- mma_bi$BR + (2*mma_sdbi$BR)
        mmam_bi$BR      <- mma_bi$BR - (2*mma_sdbi$BR)
        mmap_bu$BR      <- mma_bu$BR + (2*mma_sdbu$BR)
        mmam_bu$BR      <- mma_bu$BR - (2*mma_sdbu$BR)
        mmap_qu$BR      <- mma_qu$BR + (2*mma_sdqu$BR)
        mmam_qu$BR      <- mma_qu$BR - (2*mma_sdqu$BR)               
      }
      if("DS" %in% colnames(mma_bi)){
        mmap_bi$DS      <- mma_bi$DS + (2*mma_sdbi$DS)
        mmam_bi$DS      <- mma_bi$DS - (2*mma_sdbi$DS)
        mmap_bu$DS      <- mma_bu$DS + (2*mma_sdbu$DS)
        mmam_bu$DS      <- mma_bu$DS - (2*mma_sdbu$DS)
        mmap_qu$DS      <- mma_qu$DS + (2*mma_sdqu$DS)
        mmam_qu$DS      <- mma_qu$DS - (2*mma_sdqu$DS)               
      }
      if("ECR" %in% colnames(mma_bi)){
        mmap_bi$ECR      <- mma_bi$ECR + (2*mma_sdbi$ECR)
        mmam_bi$ECR      <- mma_bi$ECR - (2*mma_sdbi$ECR)
        mmap_bu$ECR      <- mma_bu$ECR + (2*mma_sdbu$ECR)
        mmam_bu$ECR      <- mma_bu$ECR - (2*mma_sdbu$ECR)
        mmap_qu$ECR      <- mma_qu$ECR + (2*mma_sdqu$ECR)
        mmam_qu$ECR      <- mma_qu$ECR - (2*mma_sdqu$ECR)               
      }
      if("ECU" %in% colnames(mma_bi)){
        mmap_bi$ECU      <- mma_bi$ECU + (2*mma_sdbi$ECU)
        mmam_bi$ECU      <- mma_bi$ECU - (2*mma_sdbi$ECU)
        mmap_bu$ECU      <- mma_bu$ECU + (2*mma_sdbu$ECU)
        mmam_bu$ECU      <- mma_bu$ECU - (2*mma_sdbu$ECU)
        mmap_qu$ECU      <- mma_qu$ECU + (2*mma_sdqu$ECU)
        mmam_qu$ECU      <- mma_qu$ECU - (2*mma_sdqu$ECU)               
      }
      if("EDL" %in% colnames(mma_bi)){
        mmap_bi$EDL      <- mma_bi$EDL + (2*mma_sdbi$EDL)
        mmam_bi$EDL      <- mma_bi$EDL - (2*mma_sdbi$EDL)
        mmap_bu$EDL      <- mma_bu$EDL + (2*mma_sdbu$EDL)
        mmam_bu$EDL      <- mma_bu$EDL - (2*mma_sdbu$EDL)
        mmap_qu$EDL      <- mma_qu$EDL + (2*mma_sdqu$EDL)
        mmam_qu$EDL      <- mma_qu$EDL - (2*mma_sdqu$EDL)               
      }
      if("FCR" %in% colnames(mma_bi)){
        mmap_bi$FCR      <- mma_bi$FCR + (2*mma_sdbi$FCR)
        mmam_bi$FCR      <- mma_bi$FCR - (2*mma_sdbi$FCR)
        mmap_bu$FCR      <- mma_bu$FCR + (2*mma_sdbu$FCR)
        mmam_bu$FCR      <- mma_bu$FCR - (2*mma_sdbu$FCR)
        mmap_qu$FCR      <- mma_qu$FCR + (2*mma_sdqu$FCR)
        mmam_qu$FCR      <- mma_qu$FCR - (2*mma_sdqu$FCR)               
      }
      if("FCU" %in% colnames(mma_bi)){
        mmap_bi$FCU      <- mma_bi$FCU + (2*mma_sdbi$FCU)
        mmam_bi$FCU      <- mma_bi$FCU - (2*mma_sdbi$FCU)
        mmap_bu$FCU      <- mma_bu$FCU + (2*mma_sdbu$FCU)
        mmam_bu$FCU      <- mma_bu$FCU - (2*mma_sdbu$FCU)
        mmap_qu$FCU      <- mma_qu$FCU + (2*mma_sdqu$FCU)
        mmam_qu$FCU      <- mma_qu$FCU - (2*mma_sdqu$FCU)               
      }
      if("FDL" %in% colnames(mma_bi)){
        mmap_bi$FDL      <- mma_bi$FDL + (2*mma_sdbi$FDL)
        mmam_bi$FDL      <- mma_bi$FDL - (2*mma_sdbi$FDL)
        mmap_bu$FDL      <- mma_bu$FDL + (2*mma_sdbu$FDL)
        mmam_bu$FDL      <- mma_bu$FDL - (2*mma_sdbu$FDL)
        mmap_qu$FDL      <- mma_qu$FDL + (2*mma_sdqu$FDL)
        mmam_qu$FDL      <- mma_qu$FDL - (2*mma_sdqu$FDL)               
      }
      if("FDL.U" %in% colnames(mma_bi)){
        mmap_bi$FDL.U      <- mma_bi$FDL.U + (2*mma_sdbi$FDL.U)
        mmam_bi$FDL.U      <- mma_bi$FDL.U - (2*mma_sdbi$FDL.U)
        mmap_bu$FDL.U      <- mma_bu$FDL.U + (2*mma_sdbu$FDL.U)
        mmam_bu$FDL.U      <- mma_bu$FDL.U - (2*mma_sdbu$FDL.U)
        mmap_qu$FDL.U      <- mma_qu$FDL.U + (2*mma_sdqu$FDL.U)
        mmam_qu$FDL.U      <- mma_qu$FDL.U - (2*mma_sdqu$FDL.U)               
      }
      if("HR" %in% colnames(mma_bi)){
        mmap_bi$HR      <- mma_bi$HR + (2*mma_sdbi$HR)
        mmam_bi$HR      <- mma_bi$HR - (2*mma_sdbi$HR)
        mmap_bu$HR      <- mma_bu$HR + (2*mma_sdbu$HR)
        mmam_bu$HR      <- mma_bu$HR - (2*mma_sdbu$HR)
        mmap_qu$HR      <- mma_qu$HR + (2*mma_sdqu$HR)
        mmam_qu$HR      <- mma_qu$HR - (2*mma_sdqu$HR)               
      }
      if("LD" %in% colnames(mma_bi)){
        mmap_bi$LD      <- mma_bi$LD + (2*mma_sdbi$LD)
        mmam_bi$LD      <- mma_bi$LD - (2*mma_sdbi$LD)
        mmap_bu$LD      <- mma_bu$LD + (2*mma_sdbu$LD)
        mmam_bu$LD      <- mma_bu$LD - (2*mma_sdbu$LD)
        mmap_qu$LD      <- mma_qu$LD + (2*mma_sdqu$LD)
        mmam_qu$LD      <- mma_qu$LD - (2*mma_sdqu$LD)
      }
      if("PEC" %in% colnames(mma_bi)){
        mmap_bi$PEC      <- mma_bi$PEC + (2*mma_sdbi$PEC)
        mmam_bi$PEC      <- mma_bi$PEC - (2*mma_sdbi$PEC)
        mmap_bu$PEC      <- mma_bu$PEC + (2*mma_sdbu$PEC)
        mmam_bu$PEC      <- mma_bu$PEC - (2*mma_sdbu$PEC)
        mmap_qu$PEC      <- mma_qu$PEC + (2*mma_sdqu$PEC)
        mmam_qu$PEC      <- mma_qu$PEC - (2*mma_sdqu$PEC)
      }
      if("PT" %in% colnames(mma_bi)){
        mmap_bi$PT      <- mma_bi$PT + (2*mma_sdbi$PT)
        mmam_bi$PT      <- mma_bi$PT - (2*mma_sdbi$PT)
        mmap_bu$PT      <- mma_bu$PT + (2*mma_sdbu$PT)
        mmam_bu$PT      <- mma_bu$PT - (2*mma_sdbu$PT)
        mmap_qu$PT      <- mma_qu$PT + (2*mma_sdqu$PT)
        mmam_qu$PT      <- mma_qu$PT - (2*mma_sdqu$PT)
      }
      if("SHA" %in% colnames(mma_bi)){
        mmap_bi$SHA      <- mma_bi$SHA + (2*mma_sdbi$SHA)
        mmam_bi$SHA      <- mma_bi$SHA - (2*mma_sdbi$SHA)
        mmap_bu$SHA      <- mma_bu$SHA + (2*mma_sdbu$SHA)
        mmam_bu$SHA      <- mma_bu$SHA - (2*mma_sdbu$SHA)
        mmap_qu$SHA      <- mma_qu$SHA + (2*mma_sdqu$SHA)
        mmam_qu$SHA      <- mma_qu$SHA - (2*mma_sdqu$SHA)
      }
      if("SHP" %in% colnames(mma_bi)){
        mmap_bi$SHP      <- mma_bi$SHP + (2*mma_sdbi$SHP)
        mmam_bi$SHP      <- mma_bi$SHP - (2*mma_sdbi$SHP)
        mmap_bu$SHP      <- mma_bu$SHP + (2*mma_sdbu$SHP)
        mmam_bu$SHP      <- mma_bu$SHP - (2*mma_sdbu$SHP)
        mmap_qu$SHP      <- mma_qu$SHP + (2*mma_sdqu$SHP)
        mmam_qu$SHP      <- mma_qu$SHP - (2*mma_sdqu$SHP)
      }
      if("SUBS" %in% colnames(mma_bi)){
        mmap_bi$SUBS      <- mma_bi$SUBS + (2*mma_sdbi$SUBS)
        mmam_bi$SUBS      <- mma_bi$SUBS - (2*mma_sdbi$SUBS)
        mmap_bu$SUBS      <- mma_bu$SUBS + (2*mma_sdbu$SUBS)
        mmam_bu$SUBS      <- mma_bu$SUBS - (2*mma_sdbu$SUBS)
        mmap_qu$SUBS      <- mma_qu$SUBS + (2*mma_sdqu$SUBS)
        mmam_qu$SUBS      <- mma_qu$SUBS - (2*mma_sdqu$SUBS)
      }
      if("SUP" %in% colnames(mma_bi)){
        mmap_bi$SUP      <- mma_bi$SUP + (2*mma_sdbi$SUP)
        mmam_bi$SUP      <- mma_bi$SUP - (2*mma_sdbi$SUP)
        mmap_bu$SUP      <- mma_bu$SUP + (2*mma_sdbu$SUP)
        mmam_bu$SUP      <- mma_bu$SUP - (2*mma_sdbu$SUP)
        mmap_qu$SUP      <- mma_qu$SUP + (2*mma_sdqu$SUP)
        mmam_qu$SUP      <- mma_qu$SUP - (2*mma_sdqu$SUP)
      }
      if("SUPC" %in% colnames(mma_bi)){
        mmap_bi$SUPC      <- mma_bi$SUPC + (2*mma_sdbi$SUPC)
        mmam_bi$SUPC      <- mma_bi$SUPC - (2*mma_sdbi$SUPC)
        mmap_bu$SUPC      <- mma_bu$SUPC + (2*mma_sdbu$SUPC)
        mmam_bu$SUPC      <- mma_bu$SUPC - (2*mma_sdbu$SUPC)
        mmap_qu$SUPC      <- mma_qu$SUPC + (2*mma_sdqu$SUPC)
        mmam_qu$SUPC      <- mma_qu$SUPC - (2*mma_sdqu$SUPC)
      }
      if("TM" %in% colnames(mma_bi)){
        mmap_bi$TM      <- mma_bi$TM + (2*mma_sdbi$TM)
        mmam_bi$TM      <- mma_bi$TM - (2*mma_sdbi$TM)
        mmap_bu$TM      <- mma_bu$TM + (2*mma_sdbu$TM)
        mmam_bu$TM      <- mma_bu$TM - (2*mma_sdbu$TM)
        mmap_qu$TM      <- mma_qu$TM + (2*mma_sdqu$TM)
        mmam_qu$TM      <- mma_qu$TM - (2*mma_sdqu$TM)
      }
      if("TR.S" %in% colnames(mma_bi)){
        mmap_bi$TR.S      <- mma_bi$TR.S + (2*mma_sdbi$TR.S)
        mmam_bi$TR.S      <- mma_bi$TR.S - (2*mma_sdbi$TR.S)
        mmap_bu$TR.S      <- mma_bu$TR.S + (2*mma_sdbu$TR.S)
        mmam_bu$TR.S      <- mma_bu$TR.S - (2*mma_sdbu$TR.S)
        mmap_qu$TR.S      <- mma_qu$TR.S + (2*mma_sdqu$TR.S)
        mmam_qu$TR.S      <- mma_qu$TR.S - (2*mma_sdqu$TR.S)
      }
      if("TR.C" %in% colnames(mma_bi)){
        mmap_bi$TR.C      <- mma_bi$TR.C + (2*mma_sdbi$TR.C)
        mmam_bi$TR.C      <- mma_bi$TR.C - (2*mma_sdbi$TR.C)
        mmap_bu$TR.C      <- mma_bu$TR.C + (2*mma_sdbu$TR.C)
        mmam_bu$TR.C      <- mma_bu$TR.C - (2*mma_sdbu$TR.C)
        mmam_qu$TR.C      <- mma_qu$TR.C + (2*mma_sdqu$TR.C)
        mmap_qu$TR.C      <- mma_qu$TR.C - (2*mma_sdqu$TR.C)
      }
      if("TR.M" %in% colnames(mma_bi)){
        mmap_bi$TR.M      <- mma_bi$TR.M + (2*mma_sdbi$TR.M)
        mmam_bi$TR.M      <- mma_bi$TR.M - (2*mma_sdbi$TR.M)
        mmap_bu$TR.M      <- mma_bu$TR.M + (2*mma_sdbu$TR.M)
        mmam_bu$TR.M      <- mma_bu$TR.M - (2*mma_sdbu$TR.M)
        mmap_qu$TR.M      <- mma_qu$TR.M + (2*mma_sdqu$TR.M)
        mmam_qu$TR.M      <- mma_qu$TR.M - (2*mma_sdqu$TR.M)
      }
      if("TR.L" %in% colnames(mma_bi)){
        mmap_bi$TR.L      <- mma_bi$TR.L + (2*mma_sdbi$TR.L)
        mmam_bi$TR.L      <- mma_bi$TR.L - (2*mma_sdbi$TR.L)
        mmap_bu$TR.L      <- mma_bu$TR.L + (2*mma_sdbu$TR.L)
        mmam_bu$TR.L      <- mma_bu$TR.L - (2*mma_sdbu$TR.L)
        mmap_qu$TR.L      <- mma_qu$TR.L + (2*mma_sdqu$TR.L)
        mmam_qu$TR.L      <- mma_qu$TR.L - (2*mma_sdqu$TR.L)
      }
      if("ADD" %in% colnames(mma_bi)){
        mmap_bi$ADD      <- mma_bi$ADD + (2*mma_sdbi$ADD)
        mmam_bi$ADD      <- mma_bi$ADD - (2*mma_sdbi$ADD)
        mmap_bu$ADD      <- mma_bu$ADD + (2*mma_sdbu$ADD)
        mmam_bu$ADD      <- mma_bu$ADD - (2*mma_sdbu$ADD)
        mmap_qu$ADD      <- mma_qu$ADD + (2*mma_sdqu$ADD)
        mmam_qu$ADD      <- mma_qu$ADD - (2*mma_sdqu$ADD)
      }
      if("AMB" %in% colnames(mma_bi)){
        mmap_bi$AMB      <- mma_bi$AMB + (2*mma_sdbi$AMB)
        mmam_bi$AMB      <- mma_bi$AMB - (2*mma_sdbi$AMB)
        mmap_bu$AMB      <- mma_bu$AMB + (2*mma_sdbu$AMB)
        mmam_bu$AMB      <- mma_bu$AMB - (2*mma_sdbu$AMB)
        mmap_qu$AMB      <- mma_qu$AMB + (2*mma_sdqu$AMB)
        mmam_qu$AMB      <- mma_qu$AMB - (2*mma_sdqu$AMB)
      }
      if("CFB" %in% colnames(mma_bi)){
        mmap_bi$CFB      <- mma_bi$CFB + (2*mma_sdbi$CFB)
        mmam_bi$CFB      <- mma_bi$CFB - (2*mma_sdbi$CFB)
        mmap_bu$CFB      <- mma_bu$CFB + (2*mma_sdbu$CFB)
        mmam_bu$CFB      <- mma_bu$CFB - (2*mma_sdbu$CFB)
        mmap_qu$CFB      <- mma_qu$CFB + (2*mma_sdqu$CFB)
        mmam_qu$CFB      <- mma_qu$CFB - (2*mma_sdqu$CFB)
      }
      if("FMTE" %in% colnames(mma_bi)){
        mmap_bi$FMTE      <- mma_bi$FMTE + (2*mma_sdbi$FMTE)
        mmam_bi$FMTE      <- mma_bi$FMTE - (2*mma_sdbi$FMTE)
        mmap_bu$FMTE      <- mma_bu$FMTE + (2*mma_sdbu$FMTE)
        mmam_bu$FMTE      <- mma_bu$FMTE - (2*mma_sdbu$FMTE)
        mmap_qu$FMTE      <- mma_qu$FMTE + (2*mma_sdqu$FMTE)
        mmam_qu$FMTE      <- mma_qu$FMTE - (2*mma_sdqu$FMTE)
      }
      if("FMTI" %in% colnames(mma_bi)){
        mmap_bi$FMTI      <- mma_bi$FMTI + (2*mma_sdbi$FMTI)
        mmam_bi$FMTI      <- mma_bi$FMTI - (2*mma_sdbi$FMTI)
        mmap_bu$FMTI      <- mma_bu$FMTI + (2*mma_sdbu$FMTI)
        mmam_bu$FMTI      <- mma_bu$FMTI - (2*mma_sdbu$FMTI)
        mmap_qu$FMTI      <- mma_qu$FMTI + (2*mma_sdqu$FMTI)
        mmam_qu$FMTI      <- mma_qu$FMTI - (2*mma_sdqu$FMTI)
      }
      if("FTE" %in% colnames(mma_bi)){
        mmap_bi$FTE      <- mma_bi$FTE + (2*mma_sdbi$FTE)
        mmam_bi$FTE      <- mma_bi$FTE - (2*mma_sdbi$FTE)
        mmap_bu$FTE      <- mma_bu$FTE + (2*mma_sdbu$FTE)
        mmam_bu$FTE      <- mma_bu$FTE - (2*mma_sdbu$FTE)
        mmap_qu$FTE      <- mma_qu$FTE + (2*mma_sdqu$FTE)
        mmam_qu$FTE      <- mma_qu$FTE - (2*mma_sdqu$FTE)
      }
      if("FTI.Fan" %in% colnames(mma_bi)){
        mmap_bi$FTI.Fan      <- mma_bi$FTI.Fan + (2*mma_sdbi$FTI.Fan)
        mmam_bi$FTI.Fan      <- mma_bi$FTI.Fan - (2*mma_sdbi$FTI.Fan)
        mmap_bu$FTI.Fan      <- mma_bu$FTI.Fan + (2*mma_sdbu$FTI.Fan)
        mmam_bu$FTI.Fan      <- mma_bu$FTI.Fan - (2*mma_sdbu$FTI.Fan)
        mmap_qu$FTI.Fan      <- mma_qu$FTI.Fan + (2*mma_sdqu$FTI.Fan)
        mmam_qu$FTI.Fan      <- mma_qu$FTI.Fan - (2*mma_sdqu$FTI.Fan)
      }
      if("FTI2" %in% colnames(mma_bi)){
        mmap_bi$FTI2      <- mma_bi$FTI2 + (2*mma_sdbi$FTI2)
        mmam_bi$FTI2      <- mma_bi$FTI2 - (2*mma_sdbi$FTI2)
        mmap_bu$FTI2      <- mma_bu$FTI2 + (2*mma_sdbu$FTI2)
        mmam_bu$FTI2      <- mma_bu$FTI2 - (2*mma_sdbu$FTI2)
        mmap_qu$FTI2      <- mma_qu$FTI2 + (2*mma_sdqu$FTI2)
        mmam_qu$FTI2      <- mma_qu$FTI2 - (2*mma_sdqu$FTI2)
      }
      if("IFM" %in% colnames(mma_bi)){
        mmap_bi$IFM      <- mma_bi$IFM + (2*mma_sdbi$IFM)
        mmam_bi$IFM      <- mma_bi$IFM - (2*mma_sdbi$IFM)
        mmap_bu$IFM      <- mma_bu$IFM + (2*mma_sdbu$IFM)
        mmam_bu$IFM      <- mma_bu$IFM - (2*mma_sdbu$IFM)
        mmap_qu$IFM      <- mma_qu$IFM + (2*mma_sdqu$IFM)
        mmam_qu$IFM      <- mma_qu$IFM - (2*mma_sdqu$IFM)
      }
      if("ILF" %in% colnames(mma_bi)){
        mmap_bi$ILF      <- mma_bi$ILF + (2*mma_sdbi$ILF)
        mmam_bi$ILF      <- mma_bi$ILF - (2*mma_sdbi$ILF)
        mmap_bu$ILF      <- mma_bu$ILF + (2*mma_sdbu$ILF)
        mmam_bu$ILF      <- mma_bu$ILF - (2*mma_sdbu$ILF)
        mmap_qu$ILF      <- mma_qu$ILF + (2*mma_sdqu$ILF)
        mmam_qu$ILF      <- mma_qu$ILF - (2*mma_sdqu$ILF)
      }
      if("ITB" %in% colnames(mma_bi)){
        mmap_bi$ITB      <- mma_bi$ITB + (2*mma_sdbi$ITB)
        mmam_bi$ITB      <- mma_bi$ITB - (2*mma_sdbi$ITB)
        mmap_bu$ITB      <- mma_bu$ITB + (2*mma_sdbu$ITB)
        mmam_bu$ITB      <- mma_bu$ITB - (2*mma_sdbu$ITB)
        mmap_qu$ITB      <- mma_qu$ITB + (2*mma_sdqu$ITB)
        mmam_qu$ITB      <- mma_qu$ITB - (2*mma_sdqu$ITB)
      }
      if("PIFE" %in% colnames(mma_bi)){
        mmap_bi$PIFE      <- mma_bi$PIFE + (2*mma_sdbi$PIFE)
        mmam_bi$PIFE      <- mma_bi$PIFE - (2*mma_sdbi$PIFE)
        mmap_bu$PIFE      <- mma_bu$PIFE + (2*mma_sdbu$PIFE)
        mmam_bu$PIFE      <- mma_bu$PIFE - (2*mma_sdbu$PIFE)
        mmap_qu$PIFE      <- mma_qu$PIFE + (2*mma_sdqu$PIFE)
        mmam_qu$PIFE      <- mma_qu$PIFE - (2*mma_sdqu$PIFE)
      }
      if("PIFI" %in% colnames(mma_bi)){
        mmap_bi$PIFI      <- mma_bi$PIFI + (2*mma_sdbi$PIFI)
        mmam_bi$PIFI      <- mma_bi$PIFI - (2*mma_sdbi$PIFI)
        mmap_bu$PIFI      <- mma_bu$PIFI + (2*mma_sdbu$PIFI)
        mmam_bu$PIFI      <- mma_bu$PIFI - (2*mma_sdbu$PIFI)
        mmap_qu$PIFI      <- mma_qu$PIFI + (2*mma_sdqu$PIFI)
        mmam_qu$PIFI      <- mma_qu$PIFI - (2*mma_sdqu$PIFI)
      }
      if("SUM" %in% colnames(mma_bi)){
        mmap_bi$SUM      <- 0
        mmam_bi$SUM      <- 0
        mmap_bu$SUM      <- 0
        mmam_bu$SUM      <- 0
        mmap_qu$SUM      <- 0
        mmam_qu$SUM      <- 0
        mma_bi$SUM       <- 0
        mma_bu$SUM       <- 0
        mma_qu$SUM       <- 0
      }
    }
 
    #create opposite directions
    {
      pos_bi<- mma_bi
      pos_bi[pos_bi<0] <- 0
      pos_bi_p<- mmap_bi
      pos_bi_p[pos_bi_p<0] <- 0
      pos_bi_m<- mmam_bi
      pos_bi_m[pos_bi_m<0] <- 0
      neg_bi<- mma_bi
      neg_bi[neg_bi>0] <- 0
      neg_bi<-abs(neg_bi)
      neg_bi$SUM <- 0
      neg_bi_p<- mmap_bi
      neg_bi_p[neg_bi_p>0] <- 0
      neg_bi_p<-abs(neg_bi_p)
      neg_bi_p$SUM <- 0
      neg_bi_m<- mmam_bi
      neg_bi_m[neg_bi_m>0] <- 0
      neg_bi_m<-abs(neg_bi_m)
      neg_bi_m$SUM <- 0
      
      pos_bu<- mma_bu
      pos_bu[pos_bu<0] <- 0
      pos_bu_p<- mmap_bu
      pos_bu_p[pos_bu_p<0] <- 0
      pos_bu_m<- mmam_bu
      pos_bu_m[pos_bu_m<0] <- 0
      neg_bu<- mma_bu
      neg_bu[neg_bu>0] <- 0
      neg_bu<-abs(neg_bu)
      neg_bu$SUM <- 0
      neg_bu_p<- mmap_bu
      neg_bu_p[neg_bu_p>0] <- 0
      neg_bu_p<-abs(neg_bu_p)
      neg_bu_p$SUM <- 0
      neg_bu_m<- mmam_bu
      neg_bu_m[neg_bu_m>0] <- 0
      neg_bu_m<-abs(neg_bu_m)
      neg_bu_m$SUM <- 0
      
      pos_qu<- mma_qu
      pos_qu[pos_qu<0] <- 0
      pos_qu_p<- mmap_qu
      pos_qu_p[pos_qu_p<0] <- 0
      pos_qu_m<- mmam_qu
      pos_qu_m[pos_qu_m<0] <- 0
      neg_qu<- mma_qu
      neg_qu[pos_qu>0] <- 0
      neg_qu<-abs(neg_qu)
      neg_qu$SUM <- 0
      neg_qu_p<- mmap_qu
      neg_qu_p[pos_qu_p>0] <- 0
      neg_qu_p<-abs(neg_qu_p)
      neg_qu_p$SUM <- 0
      neg_qu_m<- mmam_qu
      neg_qu_m[pos_qu_m>0] <- 0
      neg_qu_m<-abs(neg_qu_m)
      neg_qu_m$SUM <- 0
    }
       
    # new sums
    {
    pos_bi$SUM <- rowSums(pos_bi[3:(length(pos_bi))])
    pos_bi_p$SUM <- rowSums(pos_bi_p[3:(length(pos_bi_p)-1)])
    pos_bi_m$SUM <- rowSums(pos_bi_m[3:(length(pos_bi_m)-1)])
    neg_bi$SUM <- rowSums(neg_bi[3:(length(neg_bi))])
    neg_bi_p$SUM <- rowSums(neg_bi_p[3:(length(neg_bi_p)-1)])
    neg_bi_m$SUM <- rowSums(neg_bi_m[3:(length(neg_bi_m)-1)])
    
    pos_bu$SUM <- rowSums(pos_bu[3:(length(pos_bu))])
    pos_bu_p$SUM <- rowSums(pos_bu_p[3:(length(pos_bu_p)-1)])
    pos_bu_m$SUM <- rowSums(pos_bu_m[3:(length(pos_bu_m)-1)])
    neg_bu$SUM <- rowSums(neg_bu[3:(length(neg_bu))])
    neg_bu_p$SUM <- rowSums(neg_bu_p[3:(length(neg_bu_p)-1)])
    neg_bu_m$SUM <- rowSums(neg_bu_m[3:(length(neg_bu_m)-1)])
    
    pos_qu$SUM <- rowSums(pos_qu[3:(length(pos_qu))])
    pos_qu_p$SUM <- rowSums(pos_qu_p[3:(length(pos_qu_p)-1)])
    pos_qu_m$SUM <- rowSums(pos_qu_m[3:(length(pos_qu_m)-1)])
    neg_qu$SUM <- rowSums(neg_qu[3:(length(neg_qu))])
    neg_qu_p$SUM <- rowSums(neg_qu_p[3:(length(neg_qu_p)-1)])
    neg_qu_m$SUM <- rowSums(neg_qu_m[3:(length(neg_qu_m)-1)])
    
    # new means
    pos_bi$mean <- rowMeans(pos_bi[3:(length(pos_bi)-1)])
    neg_bi$mean <- rowMeans(neg_bi[3:(length(neg_bi)-1)])
    
    pos_bu$mean <- rowMeans(pos_bu[3:(length(pos_bu)-1)])
    neg_bu$mean <- rowMeans(neg_bu[3:(length(neg_bu)-1)])
    
    pos_qu$mean <- rowMeans(pos_qu[3:(length(pos_qu)-1)])
    neg_qu$mean <- rowMeans(neg_qu[3:(length(neg_qu)-1)])
    
    }
  }
  
  ##Plotting##
  exname <- "Moment Arm Figures/New/"
  zero_line        <- head(zero$SUM,-4)
  data1 <- data.frame (pos_bu$time,pos_bu$SUM,pos_bu_p$SUM,pos_bu_m$SUM,pos_bu$mean,
                       pos_bi$SUM,pos_bi_p$SUM,pos_bi_m$SUM,pos_bi$mean,
                       pos_qu$SUM,pos_qu_p$SUM,pos_qu_m$SUM,pos_qu$mean)
  data2 <- data.frame (pos_bu$time,neg_bu$SUM,neg_bu_p$SUM,neg_bu_m$SUM,neg_bu$mean,
                       neg_bi$SUM,neg_bi_p$SUM,neg_bi_m$SUM,neg_bi$mean,
                       neg_qu$SUM,neg_qu_p$SUM,neg_qu_m$SUM,neg_qu$mean)
  #remove everything after launch
  data1 <- head(data1,-4)
  data2 <- head(data2,-4)
  
  ggplot(data1,aes(x=pos_bu.time)) +
    geom_line(aes(y=zero_line,colour="black"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_bu_m.SUM,ymax=pos_bu_p.SUM),colour="#D81B60",fill="#D81B60",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_bu.SUM,colour="Burst"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_bi_m.SUM,ymax=pos_bi_p.SUM),colour="#1E88E5",fill="#1E88E5",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_bi.SUM,colour="Counter"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_qu_m.SUM,ymax=pos_qu_p.SUM),colour="#FFC107",fill="#FFC107",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_qu.SUM,colour="Quad"), linewidth=1) +
    labs(x="Time", y="Summed Moment Arms") +
    ggtitle(paste0(ti))+
    theme(plot.title=element_text(hjust=0.5))+
    scale_fill_identity(name='Legend',guide='legend',labels=c('Burst','Counter','Quad')) +
    scale_color_manual(name='Launch',values=c('Burst'='#D81B60','Counter'='#1E88E5','Quad'='#FFC107'), 
                       labels=c('Burst','Counter','Quad'))
  #ggsave((paste0(exname,ti," Moment Arms.pdf")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  #ggsave((paste0(exname,ti," Moment Arms.png")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  write.csv(data1,paste0(exname,ti," values.csv"))
  
  ggplot(data2,aes(x=pos_bu.time)) +
    geom_line(aes(y=zero_line,colour="black"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_bu_m.SUM,ymax=neg_bu_p.SUM),colour="#D81B60",fill="#D81B60",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_bu.SUM,colour="Burst"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_bi_m.SUM,ymax=neg_bi_p.SUM),colour="#1E88E5",fill="#1E88E5",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_bi.SUM,colour="Counter"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_qu_m.SUM,ymax=neg_qu_p.SUM),colour="#FFC107",fill="#FFC107",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_qu.SUM,colour="Quad"), linewidth=1) +
    labs(x="Time", y="Summed Moment Arms") +
    ggtitle(paste0(ti2))+
    theme(plot.title=element_text(hjust=0.5))+
    scale_fill_identity(name='Legend',guide='legend',labels=c('Burst','Counter','Quad'))+
    scale_color_manual(name='Launch',values=c('Burst'='#D81B60','Counter'='#1E88E5','Quad'='#FFC107'), 
                       labels=c('Burst','Counter','Quad'))
  #ggsave((paste0(exname,ti2," Moment Arms.pdf")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  #ggsave((paste0(exname,ti2," Moment Arms.png")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  write.csv(data2,paste0(exname,ti2," values.csv"))
}
#loop function for the purely avian only data (b is for bird)
graph.loopb <- function(j, df, ti, ti2)
{
  ##Load and Setup##
  {
    
    buname <- (paste0("Burst/Outputs/"))
    biname <- (paste0("Biped/Outputs/"))
    quname <- (paste0("Quad/Outputs/"))
    buname2 <- (paste0(buname,j,df))
    biname2 <- (paste0(biname,j,df))
    quname2 <- (paste0(quname,j,df))
    
    # read in raw data
    
    mma_bu   <- read.csv((paste0(buname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_bi   <- read.csv((paste0(biname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_qu   <- read.csv((paste0(quname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdbu <- read.csv((paste0(buname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdbi <- read.csv((paste0(biname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdqu <- read.csv((paste0(quname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    zero <- read.csv("zero.csv", header = TRUE, sep = ",", quote =
                       "\"",dec = ".", fill = TRUE, comment.char = "")
    
    
    #make mean+/- 2 sd
    mmap_bu  <- mma_bu
    mmam_bu  <- mma_bu
    mmap_bi  <- mma_bi
    mmam_bi  <- mma_bi
    mmap_qu  <- mma_qu
    mmam_qu  <- mma_qu
    
    {
      if("BI" %in% colnames(mma_bi)){
        mmap_bi$BI      <- mma_bi$BI + (2*mma_sdbi$BI)
        mmam_bi$BI      <- mma_bi$BI - (2*mma_sdbi$BI)
        mmap_bu$BI      <- mma_bu$BI + (2*mma_sdbu$BI)
        mmam_bu$BI      <- mma_bu$BI - (2*mma_sdbu$BI)
        mmap_qu$BI      <- mma_qu$BI + (2*mma_sdqu$BI)
        mmam_qu$BI      <- mma_qu$BI - (2*mma_sdqu$BI)
      }
      if("CB" %in% colnames(mma_bi)){
        mmap_bi$CB      <- mma_bi$CB + (2*mma_sdbi$CB)
        mmam_bi$CB      <- mma_bi$CB - (2*mma_sdbi$CB)
        mmap_bu$CB      <- mma_bu$CB + (2*mma_sdbu$CB)
        mmam_bu$CB      <- mma_bu$CB - (2*mma_sdbu$CB)
        mmap_qu$CB      <- mma_qu$CB + (2*mma_sdqu$CB)
        mmam_qu$CB      <- mma_qu$CB - (2*mma_sdqu$CB)               
      }
      if("BR" %in% colnames(mma_bi)){
        mmap_bi$BR      <- mma_bi$BR + (2*mma_sdbi$BR)
        mmam_bi$BR      <- mma_bi$BR - (2*mma_sdbi$BR)
        mmap_bu$BR      <- mma_bu$BR + (2*mma_sdbu$BR)
        mmam_bu$BR      <- mma_bu$BR - (2*mma_sdbu$BR)
        mmap_qu$BR      <- mma_qu$BR + (2*mma_sdqu$BR)
        mmam_qu$BR      <- mma_qu$BR - (2*mma_sdqu$BR)               
      }
      if("DS" %in% colnames(mma_bi)){
        mmap_bi$DS      <- mma_bi$DS + (2*mma_sdbi$DS)
        mmam_bi$DS      <- mma_bi$DS - (2*mma_sdbi$DS)
        mmap_bu$DS      <- mma_bu$DS + (2*mma_sdbu$DS)
        mmam_bu$DS      <- mma_bu$DS - (2*mma_sdbu$DS)
        mmap_qu$DS      <- mma_qu$DS + (2*mma_sdqu$DS)
        mmam_qu$DS      <- mma_qu$DS - (2*mma_sdqu$DS)               
      }
      if("ECR" %in% colnames(mma_bi)){
        mmap_bi$ECR      <- mma_bi$ECR + (2*mma_sdbi$ECR)
        mmam_bi$ECR      <- mma_bi$ECR - (2*mma_sdbi$ECR)
        mmap_bu$ECR      <- mma_bu$ECR + (2*mma_sdbu$ECR)
        mmam_bu$ECR      <- mma_bu$ECR - (2*mma_sdbu$ECR)
        mmap_qu$ECR      <- mma_qu$ECR + (2*mma_sdqu$ECR)
        mmam_qu$ECR      <- mma_qu$ECR - (2*mma_sdqu$ECR)               
      }
      if("ECU" %in% colnames(mma_bi)){
        mmap_bi$ECU      <- mma_bi$ECU + (2*mma_sdbi$ECU)
        mmam_bi$ECU      <- mma_bi$ECU - (2*mma_sdbi$ECU)
        mmap_bu$ECU      <- mma_bu$ECU + (2*mma_sdbu$ECU)
        mmam_bu$ECU      <- mma_bu$ECU - (2*mma_sdbu$ECU)
        mmap_qu$ECU      <- mma_qu$ECU + (2*mma_sdqu$ECU)
        mmam_qu$ECU      <- mma_qu$ECU - (2*mma_sdqu$ECU)               
      }
      if("EDL" %in% colnames(mma_bi)){
        mmap_bi$EDL      <- mma_bi$EDL + (2*mma_sdbi$EDL)
        mmam_bi$EDL      <- mma_bi$EDL - (2*mma_sdbi$EDL)
        mmap_bu$EDL      <- mma_bu$EDL + (2*mma_sdbu$EDL)
        mmam_bu$EDL      <- mma_bu$EDL - (2*mma_sdbu$EDL)
        mmap_qu$EDL      <- mma_qu$EDL + (2*mma_sdqu$EDL)
        mmam_qu$EDL      <- mma_qu$EDL - (2*mma_sdqu$EDL)               
      }
      if("FCR" %in% colnames(mma_bi)){
        mmap_bi$FCR      <- mma_bi$FCR + (2*mma_sdbi$FCR)
        mmam_bi$FCR      <- mma_bi$FCR - (2*mma_sdbi$FCR)
        mmap_bu$FCR      <- mma_bu$FCR + (2*mma_sdbu$FCR)
        mmam_bu$FCR      <- mma_bu$FCR - (2*mma_sdbu$FCR)
        mmap_qu$FCR      <- mma_qu$FCR + (2*mma_sdqu$FCR)
        mmam_qu$FCR      <- mma_qu$FCR - (2*mma_sdqu$FCR)               
      }
      if("FCU" %in% colnames(mma_bi)){
        mmap_bi$FCU      <- mma_bi$FCU + (2*mma_sdbi$FCU)
        mmam_bi$FCU      <- mma_bi$FCU - (2*mma_sdbi$FCU)
        mmap_bu$FCU      <- mma_bu$FCU + (2*mma_sdbu$FCU)
        mmam_bu$FCU      <- mma_bu$FCU - (2*mma_sdbu$FCU)
        mmap_qu$FCU      <- mma_qu$FCU + (2*mma_sdqu$FCU)
        mmam_qu$FCU      <- mma_qu$FCU - (2*mma_sdqu$FCU)               
      }
      if("FDL" %in% colnames(mma_bi)){
        mmap_bi$FDL      <- mma_bi$FDL + (2*mma_sdbi$FDL)
        mmam_bi$FDL      <- mma_bi$FDL - (2*mma_sdbi$FDL)
        mmap_bu$FDL      <- mma_bu$FDL + (2*mma_sdbu$FDL)
        mmam_bu$FDL      <- mma_bu$FDL - (2*mma_sdbu$FDL)
        mmap_qu$FDL      <- mma_qu$FDL + (2*mma_sdqu$FDL)
        mmam_qu$FDL      <- mma_qu$FDL - (2*mma_sdqu$FDL)               
      }
      if("FDL.U" %in% colnames(mma_bi)){
        mmap_bi$FDL.U      <- (mma_bi$FDL.U + (2*mma_sdbi$FDL.U)) *0
        mmam_bi$FDL.U      <- (mma_bi$FDL.U - (2*mma_sdbi$FDL.U)) *0
        mmap_bu$FDL.U      <- (mma_bu$FDL.U + (2*mma_sdbu$FDL.U)) *0
        mmam_bu$FDL.U      <- (mma_bu$FDL.U - (2*mma_sdbu$FDL.U)) *0
        mmap_qu$FDL.U      <- (mma_qu$FDL.U + (2*mma_sdqu$FDL.U)) *0
        mmam_qu$FDL.U      <- (mma_qu$FDL.U - (2*mma_sdqu$FDL.U)) *0 
        mma_bi$FDL.U       <- mma_bi$FDL.U *0
        mma_bu$FDL.U       <- mma_bi$FDL.U *0
        mma_qu$FDL.U       <- mma_bi$FDL.U *0
      }
      if("HR" %in% colnames(mma_bi)){
        mmap_bi$HR      <- (mma_bi$HR + (2*mma_sdbi$HR)) *0
        mmam_bi$HR      <- (mma_bi$HR - (2*mma_sdbi$HR)) *0
        mmap_bu$HR      <- (mma_bu$HR + (2*mma_sdbu$HR)) *0
        mmam_bu$HR      <- (mma_bu$HR - (2*mma_sdbu$HR)) *0
        mmap_qu$HR      <- (mma_qu$HR + (2*mma_sdqu$HR)) *0
        mmam_qu$HR      <- (mma_qu$HR - (2*mma_sdqu$HR)) *0  
        mma_bi$HR       <- mma_bi$HR *0
        mma_bu$HR       <- mma_bi$HR *0
        mma_qu$HR       <- mma_bi$HR *0
      }
      if("LD" %in% colnames(mma_bi)){
        mmap_bi$LD      <- mma_bi$LD + (2*mma_sdbi$LD)
        mmam_bi$LD      <- mma_bi$LD - (2*mma_sdbi$LD)
        mmap_bu$LD      <- mma_bu$LD + (2*mma_sdbu$LD)
        mmam_bu$LD      <- mma_bu$LD - (2*mma_sdbu$LD)
        mmap_qu$LD      <- mma_qu$LD + (2*mma_sdqu$LD)
        mmam_qu$LD      <- mma_qu$LD - (2*mma_sdqu$LD)
      }
      if("PEC" %in% colnames(mma_bi)){
        mmap_bi$PEC      <- mma_bi$PEC + (2*mma_sdbi$PEC)
        mmam_bi$PEC      <- mma_bi$PEC - (2*mma_sdbi$PEC)
        mmap_bu$PEC      <- mma_bu$PEC + (2*mma_sdbu$PEC)
        mmam_bu$PEC      <- mma_bu$PEC - (2*mma_sdbu$PEC)
        mmap_qu$PEC      <- mma_qu$PEC + (2*mma_sdqu$PEC)
        mmam_qu$PEC      <- mma_qu$PEC - (2*mma_sdqu$PEC)
      }
      if("PT" %in% colnames(mma_bi)){
        mmap_bi$PT      <- mma_bi$PT + (2*mma_sdbi$PT)
        mmam_bi$PT      <- mma_bi$PT - (2*mma_sdbi$PT)
        mmap_bu$PT      <- mma_bu$PT + (2*mma_sdbu$PT)
        mmam_bu$PT      <- mma_bu$PT - (2*mma_sdbu$PT)
        mmap_qu$PT      <- mma_qu$PT + (2*mma_sdqu$PT)
        mmam_qu$PT      <- mma_qu$PT - (2*mma_sdqu$PT)
      }
      if("SHA" %in% colnames(mma_bi)){
        mmap_bi$SHA      <- mma_bi$SHA + (2*mma_sdbi$SHA)
        mmam_bi$SHA      <- mma_bi$SHA - (2*mma_sdbi$SHA)
        mmap_bu$SHA      <- mma_bu$SHA + (2*mma_sdbu$SHA)
        mmam_bu$SHA      <- mma_bu$SHA - (2*mma_sdbu$SHA)
        mmap_qu$SHA      <- mma_qu$SHA + (2*mma_sdqu$SHA)
        mmam_qu$SHA      <- mma_qu$SHA - (2*mma_sdqu$SHA)
      }
      if("SHP" %in% colnames(mma_bi)){
        mmap_bi$SHP      <- mma_bi$SHP + (2*mma_sdbi$SHP)
        mmam_bi$SHP      <- mma_bi$SHP - (2*mma_sdbi$SHP)
        mmap_bu$SHP      <- mma_bu$SHP + (2*mma_sdbu$SHP)
        mmam_bu$SHP      <- mma_bu$SHP - (2*mma_sdbu$SHP)
        mmap_qu$SHP      <- mma_qu$SHP + (2*mma_sdqu$SHP)
        mmam_qu$SHP      <- mma_qu$SHP - (2*mma_sdqu$SHP)
      }
      if("SUBS" %in% colnames(mma_bi)){
        mmap_bi$SUBS      <- mma_bi$SUBS + (2*mma_sdbi$SUBS)
        mmam_bi$SUBS      <- mma_bi$SUBS - (2*mma_sdbi$SUBS)
        mmap_bu$SUBS      <- mma_bu$SUBS + (2*mma_sdbu$SUBS)
        mmam_bu$SUBS      <- mma_bu$SUBS - (2*mma_sdbu$SUBS)
        mmap_qu$SUBS      <- mma_qu$SUBS + (2*mma_sdqu$SUBS)
        mmam_qu$SUBS      <- mma_qu$SUBS - (2*mma_sdqu$SUBS)
      }
      if("SUP" %in% colnames(mma_bi)){
        mmap_bi$SUP      <- mma_bi$SUP + (2*mma_sdbi$SUP)
        mmam_bi$SUP      <- mma_bi$SUP - (2*mma_sdbi$SUP)
        mmap_bu$SUP      <- mma_bu$SUP + (2*mma_sdbu$SUP)
        mmam_bu$SUP      <- mma_bu$SUP - (2*mma_sdbu$SUP)
        mmap_qu$SUP      <- mma_qu$SUP + (2*mma_sdqu$SUP)
        mmam_qu$SUP      <- mma_qu$SUP - (2*mma_sdqu$SUP)
      }
      if("SUPC" %in% colnames(mma_bi)){
        mmap_bi$SUPC      <- mma_bi$SUPC + (2*mma_sdbi$SUPC)
        mmam_bi$SUPC      <- mma_bi$SUPC - (2*mma_sdbi$SUPC)
        mmap_bu$SUPC      <- mma_bu$SUPC + (2*mma_sdbu$SUPC)
        mmam_bu$SUPC      <- mma_bu$SUPC - (2*mma_sdbu$SUPC)
        mmap_qu$SUPC      <- mma_qu$SUPC + (2*mma_sdqu$SUPC)
        mmam_qu$SUPC      <- mma_qu$SUPC - (2*mma_sdqu$SUPC)
      }
      if("TM" %in% colnames(mma_bi)){
        mmap_bi$TM      <- (mma_bi$TM + (2*mma_sdbi$TM)) *0
        mmam_bi$TM      <- (mma_bi$TM - (2*mma_sdbi$TM)) *0
        mmap_bu$TM      <- (mma_bu$TM + (2*mma_sdbu$TM)) *0
        mmam_bu$TM      <- (mma_bu$TM - (2*mma_sdbu$TM)) *0
        mmap_qu$TM      <- (mma_qu$TM + (2*mma_sdqu$TM)) *0
        mmam_qu$TM      <- (mma_qu$TM - (2*mma_sdqu$TM)) *0
        mma_bi$TM       <- mma_bi$TM *0
        mma_bu$TM       <- mma_bi$TM *0
        mma_qu$TM       <- mma_bi$TM *0
      }
      if("TR.S" %in% colnames(mma_bi)){
        mmap_bi$TR.S      <- mma_bi$TR.S + (2*mma_sdbi$TR.S)
        mmam_bi$TR.S      <- mma_bi$TR.S - (2*mma_sdbi$TR.S)
        mmap_bu$TR.S      <- mma_bu$TR.S + (2*mma_sdbu$TR.S)
        mmam_bu$TR.S      <- mma_bu$TR.S - (2*mma_sdbu$TR.S)
        mmap_qu$TR.S      <- mma_qu$TR.S + (2*mma_sdqu$TR.S)
        mmam_qu$TR.S      <- mma_qu$TR.S - (2*mma_sdqu$TR.S)
      }
      if("TR.C" %in% colnames(mma_bi)){
        mmap_bi$TR.C      <- (mma_bi$TR.C + (2*mma_sdbi$TR.C)) *0
        mmam_bi$TR.C      <- (mma_bi$TR.C - (2*mma_sdbi$TR.C)) *0
        mmap_bu$TR.C      <- (mma_bu$TR.C + (2*mma_sdbu$TR.C)) *0
        mmam_bu$TR.C      <- (mma_bu$TR.C - (2*mma_sdbu$TR.C)) *0
        mmam_qu$TR.C      <- (mma_qu$TR.C + (2*mma_sdqu$TR.C)) *0
        mmap_qu$TR.C      <- (mma_qu$TR.C - (2*mma_sdqu$TR.C)) *0
        mma_bi$TR.C       <- mma_bi$TR.C *0
        mma_bu$TR.C       <- mma_bi$TR.C *0
        mma_qu$TR.C       <- mma_bi$TR.C *0
      }
      if("TR.M" %in% colnames(mma_bi)){
        mmap_bi$TR.M      <- (mma_bi$TR.M + (2*mma_sdbi$TR.M)) *0
        mmam_bi$TR.M      <- (mma_bi$TR.M - (2*mma_sdbi$TR.M)) *0
        mmap_bu$TR.M      <- (mma_bu$TR.M + (2*mma_sdbu$TR.M)) *0
        mmam_bu$TR.M      <- (mma_bu$TR.M - (2*mma_sdbu$TR.M)) *0
        mmap_qu$TR.M      <- (mma_qu$TR.M + (2*mma_sdqu$TR.M)) *0
        mmam_qu$TR.M      <- (mma_qu$TR.M - (2*mma_sdqu$TR.M)) *0
        mma_bi$TR.M       <- mma_bi$TR.M *0
        mma_bu$TR.M       <- mma_bi$TR.M *0
        mma_qu$TR.M       <- mma_bi$TR.M *0
      }
      if("TR.L" %in% colnames(mma_bi)){
        mmap_bi$TR.L      <- mma_bi$TR.L + (2*mma_sdbi$TR.L)
        mmam_bi$TR.L      <- mma_bi$TR.L - (2*mma_sdbi$TR.L)
        mmap_bu$TR.L      <- mma_bu$TR.L + (2*mma_sdbu$TR.L)
        mmam_bu$TR.L      <- mma_bu$TR.L - (2*mma_sdbu$TR.L)
        mmap_qu$TR.L      <- mma_qu$TR.L + (2*mma_sdqu$TR.L)
        mmam_qu$TR.L      <- mma_qu$TR.L - (2*mma_sdqu$TR.L)
      }
      if("ADD" %in% colnames(mma_bi)){
        mmap_bi$ADD      <- mma_bi$ADD + (2*mma_sdbi$ADD)
        mmam_bi$ADD      <- mma_bi$ADD - (2*mma_sdbi$ADD)
        mmap_bu$ADD      <- mma_bu$ADD + (2*mma_sdbu$ADD)
        mmam_bu$ADD      <- mma_bu$ADD - (2*mma_sdbu$ADD)
        mmap_qu$ADD      <- mma_qu$ADD + (2*mma_sdqu$ADD)
        mmam_qu$ADD      <- mma_qu$ADD - (2*mma_sdqu$ADD)
      }
      if("AMB" %in% colnames(mma_bi)){
        mmap_bi$AMB      <- mma_bi$AMB + (2*mma_sdbi$AMB)
        mmam_bi$AMB      <- mma_bi$AMB - (2*mma_sdbi$AMB)
        mmap_bu$AMB      <- mma_bu$AMB + (2*mma_sdbu$AMB)
        mmam_bu$AMB      <- mma_bu$AMB - (2*mma_sdbu$AMB)
        mmap_qu$AMB      <- mma_qu$AMB + (2*mma_sdqu$AMB)
        mmam_qu$AMB      <- mma_qu$AMB - (2*mma_sdqu$AMB)
      }
      if("CFB" %in% colnames(mma_bi)){
        mmap_bi$CFB      <- mma_bi$CFB + (2*mma_sdbi$CFB)
        mmam_bi$CFB      <- mma_bi$CFB - (2*mma_sdbi$CFB)
        mmap_bu$CFB      <- mma_bu$CFB + (2*mma_sdbu$CFB)
        mmam_bu$CFB      <- mma_bu$CFB - (2*mma_sdbu$CFB)
        mmap_qu$CFB      <- mma_qu$CFB + (2*mma_sdqu$CFB)
        mmam_qu$CFB      <- mma_qu$CFB - (2*mma_sdqu$CFB)
      }
      if("FMTE" %in% colnames(mma_bi)){
        mmap_bi$FMTE      <- mma_bi$FMTE + (2*mma_sdbi$FMTE)
        mmam_bi$FMTE      <- mma_bi$FMTE - (2*mma_sdbi$FMTE)
        mmap_bu$FMTE      <- mma_bu$FMTE + (2*mma_sdbu$FMTE)
        mmam_bu$FMTE      <- mma_bu$FMTE - (2*mma_sdbu$FMTE)
        mmap_qu$FMTE      <- mma_qu$FMTE + (2*mma_sdqu$FMTE)
        mmam_qu$FMTE      <- mma_qu$FMTE - (2*mma_sdqu$FMTE)
      }
      if("FMTI" %in% colnames(mma_bi)){
        mmap_bi$FMTI      <- mma_bi$FMTI + (2*mma_sdbi$FMTI)
        mmam_bi$FMTI      <- mma_bi$FMTI - (2*mma_sdbi$FMTI)
        mmap_bu$FMTI      <- mma_bu$FMTI + (2*mma_sdbu$FMTI)
        mmam_bu$FMTI      <- mma_bu$FMTI - (2*mma_sdbu$FMTI)
        mmap_qu$FMTI      <- mma_qu$FMTI + (2*mma_sdqu$FMTI)
        mmam_qu$FMTI      <- mma_qu$FMTI - (2*mma_sdqu$FMTI)
      }
      if("FTE" %in% colnames(mma_bi)){
        mmap_bi$FTE      <- mma_bi$FTE + (2*mma_sdbi$FTE)
        mmam_bi$FTE      <- mma_bi$FTE - (2*mma_sdbi$FTE)
        mmap_bu$FTE      <- mma_bu$FTE + (2*mma_sdbu$FTE)
        mmam_bu$FTE      <- mma_bu$FTE - (2*mma_sdbu$FTE)
        mmap_qu$FTE      <- mma_qu$FTE + (2*mma_sdqu$FTE)
        mmam_qu$FTE      <- mma_qu$FTE - (2*mma_sdqu$FTE)
      }
      if("FTI.Fan" %in% colnames(mma_bi)){
        mmap_bi$FTI.Fan      <- mma_bi$FTI.Fan + (2*mma_sdbi$FTI.Fan)
        mmam_bi$FTI.Fan      <- mma_bi$FTI.Fan - (2*mma_sdbi$FTI.Fan)
        mmap_bu$FTI.Fan      <- mma_bu$FTI.Fan + (2*mma_sdbu$FTI.Fan)
        mmam_bu$FTI.Fan      <- mma_bu$FTI.Fan - (2*mma_sdbu$FTI.Fan)
        mmap_qu$FTI.Fan      <- mma_qu$FTI.Fan + (2*mma_sdqu$FTI.Fan)
        mmam_qu$FTI.Fan      <- mma_qu$FTI.Fan - (2*mma_sdqu$FTI.Fan)
      }
      if("FTI2" %in% colnames(mma_bi)){
        mmap_bi$FTI2      <- mma_bi$FTI2 + (2*mma_sdbi$FTI2)
        mmam_bi$FTI2      <- mma_bi$FTI2 - (2*mma_sdbi$FTI2)
        mmap_bu$FTI2      <- mma_bu$FTI2 + (2*mma_sdbu$FTI2)
        mmam_bu$FTI2      <- mma_bu$FTI2 - (2*mma_sdbu$FTI2)
        mmap_qu$FTI2      <- mma_qu$FTI2 + (2*mma_sdqu$FTI2)
        mmam_qu$FTI2      <- mma_qu$FTI2 - (2*mma_sdqu$FTI2)
      }
      if("IFM" %in% colnames(mma_bi)){
        mmap_bi$IFM      <- mma_bi$IFM + (2*mma_sdbi$IFM)
        mmam_bi$IFM      <- mma_bi$IFM - (2*mma_sdbi$IFM)
        mmap_bu$IFM      <- mma_bu$IFM + (2*mma_sdbu$IFM)
        mmam_bu$IFM      <- mma_bu$IFM - (2*mma_sdbu$IFM)
        mmap_qu$IFM      <- mma_qu$IFM + (2*mma_sdqu$IFM)
        mmam_qu$IFM      <- mma_qu$IFM - (2*mma_sdqu$IFM)
      }
      if("ILF" %in% colnames(mma_bi)){
        mmap_bi$ILF      <- mma_bi$ILF + (2*mma_sdbi$ILF)
        mmam_bi$ILF      <- mma_bi$ILF - (2*mma_sdbi$ILF)
        mmap_bu$ILF      <- mma_bu$ILF + (2*mma_sdbu$ILF)
        mmam_bu$ILF      <- mma_bu$ILF - (2*mma_sdbu$ILF)
        mmap_qu$ILF      <- mma_qu$ILF + (2*mma_sdqu$ILF)
        mmam_qu$ILF      <- mma_qu$ILF - (2*mma_sdqu$ILF)
      }
      if("ITB" %in% colnames(mma_bi)){
        mmap_bi$ITB      <- mma_bi$ITB + (2*mma_sdbi$ITB)
        mmam_bi$ITB      <- mma_bi$ITB - (2*mma_sdbi$ITB)
        mmap_bu$ITB      <- mma_bu$ITB + (2*mma_sdbu$ITB)
        mmam_bu$ITB      <- mma_bu$ITB - (2*mma_sdbu$ITB)
        mmap_qu$ITB      <- mma_qu$ITB + (2*mma_sdqu$ITB)
        mmam_qu$ITB      <- mma_qu$ITB - (2*mma_sdqu$ITB)
      }
      if("PIFE" %in% colnames(mma_bi)){
        mmap_bi$PIFE      <- mma_bi$PIFE + (2*mma_sdbi$PIFE)
        mmam_bi$PIFE      <- mma_bi$PIFE - (2*mma_sdbi$PIFE)
        mmap_bu$PIFE      <- mma_bu$PIFE + (2*mma_sdbu$PIFE)
        mmam_bu$PIFE      <- mma_bu$PIFE - (2*mma_sdbu$PIFE)
        mmap_qu$PIFE      <- mma_qu$PIFE + (2*mma_sdqu$PIFE)
        mmam_qu$PIFE      <- mma_qu$PIFE - (2*mma_sdqu$PIFE)
      }
      if("PIFI" %in% colnames(mma_bi)){
        mmap_bi$PIFI      <- mma_bi$PIFI + (2*mma_sdbi$PIFI)
        mmam_bi$PIFI      <- mma_bi$PIFI - (2*mma_sdbi$PIFI)
        mmap_bu$PIFI      <- mma_bu$PIFI + (2*mma_sdbu$PIFI)
        mmam_bu$PIFI      <- mma_bu$PIFI - (2*mma_sdbu$PIFI)
        mmap_qu$PIFI      <- mma_qu$PIFI + (2*mma_sdqu$PIFI)
        mmam_qu$PIFI      <- mma_qu$PIFI - (2*mma_sdqu$PIFI)
      }
      if("SUM" %in% colnames(mma_bi)){
        mmap_bi$SUM      <- 0
        mmam_bi$SUM      <- 0
        mmap_bu$SUM      <- 0
        mmam_bu$SUM      <- 0
        mmap_qu$SUM      <- 0
        mmam_qu$SUM      <- 0
        mma_bi$SUM       <- 0
        mma_bu$SUM       <- 0
        mma_qu$SUM       <- 0
      }
    }

    #create opposite directions
    {
      pos_bi<- mma_bi
      pos_bi[pos_bi<0] <- 0
      pos_bi_p<- mmap_bi
      pos_bi_p[pos_bi_p<0] <- 0
      pos_bi_m<- mmam_bi
      pos_bi_m[pos_bi_m<0] <- 0
      neg_bi<- mma_bi
      neg_bi[neg_bi>0] <- 0
      neg_bi<-abs(neg_bi)
      neg_bi$SUM <- 0
      neg_bi_p<- mmap_bi
      neg_bi_p[neg_bi_p>0] <- 0
      neg_bi_p<-abs(neg_bi_p)
      neg_bi_p$SUM <- 0
      neg_bi_m<- mmam_bi
      neg_bi_m[neg_bi_m>0] <- 0
      neg_bi_m<-abs(neg_bi_m)
      neg_bi_m$SUM <- 0
      
      pos_bu<- mma_bu
      pos_bu[pos_bu<0] <- 0
      pos_bu_p<- mmap_bu
      pos_bu_p[pos_bu_p<0] <- 0
      pos_bu_m<- mmam_bu
      pos_bu_m[pos_bu_m<0] <- 0
      neg_bu<- mma_bu
      neg_bu[neg_bu>0] <- 0
      neg_bu<-abs(neg_bu)
      neg_bu$SUM <- 0
      neg_bu_p<- mmap_bu
      neg_bu_p[neg_bu_p>0] <- 0
      neg_bu_p<-abs(neg_bu_p)
      neg_bu_p$SUM <- 0
      neg_bu_m<- mmam_bu
      neg_bu_m[neg_bu_m>0] <- 0
      neg_bu_m<-abs(neg_bu_m)
      neg_bu_m$SUM <- 0
      
      pos_qu<- mma_qu
      pos_qu[pos_qu<0] <- 0
      pos_qu_p<- mmap_qu
      pos_qu_p[pos_qu_p<0] <- 0
      pos_qu_m<- mmam_qu
      pos_qu_m[pos_qu_m<0] <- 0
      neg_qu<- mma_qu
      neg_qu[pos_qu>0] <- 0
      neg_qu<-abs(neg_qu)
      neg_qu$SUM <- 0
      neg_qu_p<- mmap_qu
      neg_qu_p[pos_qu_p>0] <- 0
      neg_qu_p<-abs(neg_qu_p)
      neg_qu_p$SUM <- 0
      neg_qu_m<- mmam_qu
      neg_qu_m[pos_qu_m>0] <- 0
      neg_qu_m<-abs(neg_qu_m)
      neg_qu_m$SUM <- 0
    }
    
    # new sums
    {
      pos_bi$SUM <- rowSums(pos_bi[3:(length(pos_bi))])
      pos_bi_p$SUM <- rowSums(pos_bi_p[3:(length(pos_bi_p)-1)])
      pos_bi_m$SUM <- rowSums(pos_bi_m[3:(length(pos_bi_m)-1)])
      neg_bi$SUM <- rowSums(neg_bi[3:(length(neg_bi))])
      neg_bi_p$SUM <- rowSums(neg_bi_p[3:(length(neg_bi_p)-1)])
      neg_bi_m$SUM <- rowSums(neg_bi_m[3:(length(neg_bi_m)-1)])
      
      pos_bu$SUM <- rowSums(pos_bu[3:(length(pos_bu))])
      pos_bu_p$SUM <- rowSums(pos_bu_p[3:(length(pos_bu_p)-1)])
      pos_bu_m$SUM <- rowSums(pos_bu_m[3:(length(pos_bu_m)-1)])
      neg_bu$SUM <- rowSums(neg_bu[3:(length(neg_bu))])
      neg_bu_p$SUM <- rowSums(neg_bu_p[3:(length(neg_bu_p)-1)])
      neg_bu_m$SUM <- rowSums(neg_bu_m[3:(length(neg_bu_m)-1)])
      
      pos_qu$SUM <- rowSums(pos_qu[3:(length(pos_qu))])
      pos_qu_p$SUM <- rowSums(pos_qu_p[3:(length(pos_qu_p)-1)])
      pos_qu_m$SUM <- rowSums(pos_qu_m[3:(length(pos_qu_m)-1)])
      neg_qu$SUM <- rowSums(neg_qu[3:(length(neg_qu))])
      neg_qu_p$SUM <- rowSums(neg_qu_p[3:(length(neg_qu_p)-1)])
      neg_qu_m$SUM <- rowSums(neg_qu_m[3:(length(neg_qu_m)-1)])
    }
  }
  
  ##Plotting##
  exname <- "Moment Arm Figures/New/"
  zero_line        <- head(zero$SUM,-4)
  data1 <- data.frame (pos_bu$time,pos_bu$SUM,pos_bu_p$SUM,pos_bu_m$SUM,
                       pos_bi$SUM,pos_bi_p$SUM,pos_bi_m$SUM,
                       pos_qu$SUM,pos_qu_p$SUM,pos_qu_m$SUM)
  data2 <- data.frame (pos_bu$time,neg_bu$SUM,neg_bu_p$SUM,neg_bu_m$SUM,
                       neg_bi$SUM,neg_bi_p$SUM,neg_bi_m$SUM,
                       neg_qu$SUM,neg_qu_p$SUM,neg_qu_m$SUM)
  #remove everything after launch
  data1 <- head(data1,-4)
  data2 <- head(data2,-4)
  
  ggplot(data1,aes(x=pos_bu.time)) +
    geom_line(aes(y=zero_line,colour="black"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_bu_m.SUM,ymax=pos_bu_p.SUM),colour="#D81B60",fill="#D81B60",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_bu.SUM,colour="Burst"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_bi_m.SUM,ymax=pos_bi_p.SUM),colour="#1E88E5",fill="#1E88E5",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_bi.SUM,colour="Counter"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_qu_m.SUM,ymax=pos_qu_p.SUM),colour="#FFC107",fill="#FFC107",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_qu.SUM,colour="Quad"), linewidth=1) +
    labs(x="Time", y="Summed Moment Arms") +
    ggtitle(paste0(ti))+
    theme(plot.title=element_text(hjust=0.5))+
    scale_fill_identity(name='Legend',guide='legend',labels=c('Burst','Counter','Quad')) +
    scale_color_manual(name='Launch',values=c('Burst'='#D81B60','Counter'='#1E88E5','Quad'='#FFC107'), 
                       labels=c('Burst','Counter','Quad'))
  #ggsave((paste0(exname,ti," Moment Arms.pdf")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  #ggsave((paste0(exname,ti," Moment Arms.png")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  write.csv(data1,paste0(exname,ti," values.csv"))
  
  ggplot(data2,aes(x=pos_bu.time)) +
    geom_line(aes(y=zero_line,colour="black"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_bu_m.SUM,ymax=neg_bu_p.SUM),colour="#D81B60",fill="#D81B60",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_bu.SUM,colour="Burst"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_bi_m.SUM,ymax=neg_bi_p.SUM),colour="#1E88E5",fill="#1E88E5",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_bi.SUM,colour="Counter"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_qu_m.SUM,ymax=neg_qu_p.SUM),colour="#FFC107",fill="#FFC107",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_qu.SUM,colour="Quad"), linewidth=1) +
    labs(x="Time", y="Summed Moment Arms") +
    ggtitle(paste0(ti2))+
    theme(plot.title=element_text(hjust=0.5))+
    scale_fill_identity(name='Legend',guide='legend',labels=c('Burst','Counter','Quad'))+
    scale_color_manual(name='Launch',values=c('Burst'='#D81B60','Counter'='#1E88E5','Quad'='#FFC107'), 
                       labels=c('Burst','Counter','Quad'))
  #ggsave((paste0(exname,ti2," Moment Arms.pdf")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  #ggsave((paste0(exname,ti2," Moment Arms.png")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  write.csv(data2,paste0(exname,ti2," values.csv"))
} 

#loop function for the crocodilian only data (c is for croc)
graph.loopc <- function(j, df, ti, ti2)
{
  ##Load and Setup##
  {
    
    buname <- (paste0("Burst/Outputs/"))
    biname <- (paste0("Biped/Outputs/"))
    quname <- (paste0("Quad/Outputs/"))
    buname2 <- (paste0(buname,j,df))
    biname2 <- (paste0(biname,j,df))
    quname2 <- (paste0(quname,j,df))
    
    # read in raw data
    
    mma_bu   <- read.csv((paste0(buname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_bi   <- read.csv((paste0(biname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_qu   <- read.csv((paste0(quname2,"_MMA_MC_mean.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdbu <- read.csv((paste0(buname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdbi <- read.csv((paste0(biname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    mma_sdqu <- read.csv((paste0(quname2,"_MMA_MC_stdev.csv")), header = TRUE, sep = ",", quote =
                           "\"",dec = ".", fill = TRUE, comment.char = "")
    zero <- read.csv("zero.csv", header = TRUE, sep = ",", quote =
                       "\"",dec = ".", fill = TRUE, comment.char = "")
    
    
    #make mean+/- 2 sd
    mmap_bu  <- mma_bu
    mmam_bu  <- mma_bu
    mmap_bi  <- mma_bi
    mmam_bi  <- mma_bi
    mmap_qu  <- mma_qu
    mmam_qu  <- mma_qu
    
    {
      if("BI" %in% colnames(mma_bi)){
        mmap_bi$BI      <- mma_bi$BI + (2*mma_sdbi$BI)
        mmam_bi$BI      <- mma_bi$BI - (2*mma_sdbi$BI)
        mmap_bu$BI      <- mma_bu$BI + (2*mma_sdbu$BI)
        mmam_bu$BI      <- mma_bu$BI - (2*mma_sdbu$BI)
        mmap_qu$BI      <- mma_qu$BI + (2*mma_sdqu$BI)
        mmam_qu$BI      <- mma_qu$BI - (2*mma_sdqu$BI)
      }
      if("CB" %in% colnames(mma_bi)){
        mmap_bi$CB      <- mma_bi$CB + (2*mma_sdbi$CB)
        mmam_bi$CB      <- mma_bi$CB - (2*mma_sdbi$CB)
        mmap_bu$CB      <- mma_bu$CB + (2*mma_sdbu$CB)
        mmam_bu$CB      <- mma_bu$CB - (2*mma_sdbu$CB)
        mmap_qu$CB      <- mma_qu$CB + (2*mma_sdqu$CB)
        mmam_qu$CB      <- mma_qu$CB - (2*mma_sdqu$CB)               
      }
      if("BR" %in% colnames(mma_bi)){
        mmap_bi$BR      <- mma_bi$BR + (2*mma_sdbi$BR)
        mmam_bi$BR      <- mma_bi$BR - (2*mma_sdbi$BR)
        mmap_bu$BR      <- mma_bu$BR + (2*mma_sdbu$BR)
        mmam_bu$BR      <- mma_bu$BR - (2*mma_sdbu$BR)
        mmap_qu$BR      <- mma_qu$BR + (2*mma_sdqu$BR)
        mmam_qu$BR      <- mma_qu$BR - (2*mma_sdqu$BR)               
      }
      if("DS" %in% colnames(mma_bi)){
        mmap_bi$DS      <- mma_bi$DS + (2*mma_sdbi$DS)
        mmam_bi$DS      <- mma_bi$DS - (2*mma_sdbi$DS)
        mmap_bu$DS      <- mma_bu$DS + (2*mma_sdbu$DS)
        mmam_bu$DS      <- mma_bu$DS - (2*mma_sdbu$DS)
        mmap_qu$DS      <- mma_qu$DS + (2*mma_sdqu$DS)
        mmam_qu$DS      <- mma_qu$DS - (2*mma_sdqu$DS)               
      }
      if("ECR" %in% colnames(mma_bi)){
        mmap_bi$ECR      <- mma_bi$ECR + (2*mma_sdbi$ECR)
        mmam_bi$ECR      <- mma_bi$ECR - (2*mma_sdbi$ECR)
        mmap_bu$ECR      <- mma_bu$ECR + (2*mma_sdbu$ECR)
        mmam_bu$ECR      <- mma_bu$ECR - (2*mma_sdbu$ECR)
        mmap_qu$ECR      <- mma_qu$ECR + (2*mma_sdqu$ECR)
        mmam_qu$ECR      <- mma_qu$ECR - (2*mma_sdqu$ECR)               
      }
      if("ECU" %in% colnames(mma_bi)){
        mmap_bi$ECU      <- mma_bi$ECU + (2*mma_sdbi$ECU)
        mmam_bi$ECU      <- mma_bi$ECU - (2*mma_sdbi$ECU)
        mmap_bu$ECU      <- mma_bu$ECU + (2*mma_sdbu$ECU)
        mmam_bu$ECU      <- mma_bu$ECU - (2*mma_sdbu$ECU)
        mmap_qu$ECU      <- mma_qu$ECU + (2*mma_sdqu$ECU)
        mmam_qu$ECU      <- mma_qu$ECU - (2*mma_sdqu$ECU)               
      }
      if("EDL" %in% colnames(mma_bi)){
        mmap_bi$EDL      <- mma_bi$EDL + (2*mma_sdbi$EDL)
        mmam_bi$EDL      <- mma_bi$EDL - (2*mma_sdbi$EDL)
        mmap_bu$EDL      <- mma_bu$EDL + (2*mma_sdbu$EDL)
        mmam_bu$EDL      <- mma_bu$EDL - (2*mma_sdbu$EDL)
        mmap_qu$EDL      <- mma_qu$EDL + (2*mma_sdqu$EDL)
        mmam_qu$EDL      <- mma_qu$EDL - (2*mma_sdqu$EDL)               
      }
      if("FCR" %in% colnames(mma_bi)){
        mmap_bi$FCR      <- (mma_bi$FCR + (2*mma_sdbi$FCR))
        mmam_bi$FCR      <- (mma_bi$FCR - (2*mma_sdbi$FCR))
        mmap_bu$FCR      <- (mma_bu$FCR + (2*mma_sdbu$FCR))
        mmam_bu$FCR      <- (mma_bu$FCR - (2*mma_sdbu$FCR))
        mmap_qu$FCR      <- (mma_qu$FCR + (2*mma_sdqu$FCR))
        mmam_qu$FCR      <- (mma_qu$FCR - (2*mma_sdqu$FCR))
        mma_bi$FCR       <- mma_bi$FCR *0
        mma_bu$FCR       <- mma_bu$FCR *0
        mma_qu$FCR       <- mma_qu$FCR *0
      }
      if("FCU" %in% colnames(mma_bi)){
        mmap_bi$FCU      <- mma_bi$FCU + (2*mma_sdbi$FCU)
        mmam_bi$FCU      <- mma_bi$FCU - (2*mma_sdbi$FCU)
        mmap_bu$FCU      <- mma_bu$FCU + (2*mma_sdbu$FCU)
        mmam_bu$FCU      <- mma_bu$FCU - (2*mma_sdbu$FCU)
        mmap_qu$FCU      <- mma_qu$FCU + (2*mma_sdqu$FCU)
        mmam_qu$FCU      <- mma_qu$FCU - (2*mma_sdqu$FCU)               
      }
      if("FDL" %in% colnames(mma_bi)){
        mmap_bi$FDL      <- mma_bi$FDL + (2*mma_sdbi$FDL)
        mmam_bi$FDL      <- mma_bi$FDL - (2*mma_sdbi$FDL)
        mmap_bu$FDL      <- mma_bu$FDL + (2*mma_sdbu$FDL)
        mmam_bu$FDL      <- mma_bu$FDL - (2*mma_sdbu$FDL)
        mmap_qu$FDL      <- mma_qu$FDL + (2*mma_sdqu$FDL)
        mmam_qu$FDL      <- mma_qu$FDL - (2*mma_sdqu$FDL)               
      }
      if("FDL.U" %in% colnames(mma_bi)){
        mmap_bi$FDL.U      <- mma_bi$FDL.U + (2*mma_sdbi$FDL.U)
        mmam_bi$FDL.U      <- mma_bi$FDL.U - (2*mma_sdbi$FDL.U)
        mmap_bu$FDL.U      <- mma_bu$FDL.U + (2*mma_sdbu$FDL.U)
        mmam_bu$FDL.U      <- mma_bu$FDL.U - (2*mma_sdbu$FDL.U)
        mmap_qu$FDL.U      <- mma_qu$FDL.U + (2*mma_sdqu$FDL.U)
        mmam_qu$FDL.U      <- mma_qu$FDL.U - (2*mma_sdqu$FDL.U)               
      }
      if("HR" %in% colnames(mma_bi)){
        mmap_bi$HR      <- mma_bi$HR + (2*mma_sdbi$HR)
        mmam_bi$HR      <- mma_bi$HR - (2*mma_sdbi$HR)
        mmap_bu$HR      <- mma_bu$HR + (2*mma_sdbu$HR)
        mmam_bu$HR      <- mma_bu$HR - (2*mma_sdbu$HR)
        mmap_qu$HR      <- mma_qu$HR + (2*mma_sdqu$HR)
        mmam_qu$HR      <- mma_qu$HR - (2*mma_sdqu$HR)               
      }
      if("LD" %in% colnames(mma_bi)){
        mmap_bi$LD      <- mma_bi$LD + (2*mma_sdbi$LD)
        mmam_bi$LD      <- mma_bi$LD - (2*mma_sdbi$LD)
        mmap_bu$LD      <- mma_bu$LD + (2*mma_sdbu$LD)
        mmam_bu$LD      <- mma_bu$LD - (2*mma_sdbu$LD)
        mmap_qu$LD      <- mma_qu$LD + (2*mma_sdqu$LD)
        mmam_qu$LD      <- mma_qu$LD - (2*mma_sdqu$LD)
      }
      if("PEC" %in% colnames(mma_bi)){
        mmap_bi$PEC      <- mma_bi$PEC + (2*mma_sdbi$PEC)
        mmam_bi$PEC      <- mma_bi$PEC - (2*mma_sdbi$PEC)
        mmap_bu$PEC      <- mma_bu$PEC + (2*mma_sdbu$PEC)
        mmam_bu$PEC      <- mma_bu$PEC - (2*mma_sdbu$PEC)
        mmap_qu$PEC      <- mma_qu$PEC + (2*mma_sdqu$PEC)
        mmam_qu$PEC      <- mma_qu$PEC - (2*mma_sdqu$PEC)
      }
      if("PT" %in% colnames(mma_bi)){
        mmap_bi$PT      <- mma_bi$PT + (2*mma_sdbi$PT)
        mmam_bi$PT      <- mma_bi$PT - (2*mma_sdbi$PT)
        mmap_bu$PT      <- mma_bu$PT + (2*mma_sdbu$PT)
        mmam_bu$PT      <- mma_bu$PT - (2*mma_sdbu$PT)
        mmap_qu$PT      <- mma_qu$PT + (2*mma_sdqu$PT)
        mmam_qu$PT      <- mma_qu$PT - (2*mma_sdqu$PT)
      }
      if("SHA" %in% colnames(mma_bi)){
        mmap_bi$SHA      <- mma_bi$SHA + (2*mma_sdbi$SHA)
        mmam_bi$SHA      <- mma_bi$SHA - (2*mma_sdbi$SHA)
        mmap_bu$SHA      <- mma_bu$SHA + (2*mma_sdbu$SHA)
        mmam_bu$SHA      <- mma_bu$SHA - (2*mma_sdbu$SHA)
        mmap_qu$SHA      <- mma_qu$SHA + (2*mma_sdqu$SHA)
        mmam_qu$SHA      <- mma_qu$SHA - (2*mma_sdqu$SHA)
      }
      if("SHP" %in% colnames(mma_bi)){
        mmap_bi$SHP      <- mma_bi$SHP + (2*mma_sdbi$SHP)
        mmam_bi$SHP      <- mma_bi$SHP - (2*mma_sdbi$SHP)
        mmap_bu$SHP      <- mma_bu$SHP + (2*mma_sdbu$SHP)
        mmam_bu$SHP      <- mma_bu$SHP - (2*mma_sdbu$SHP)
        mmap_qu$SHP      <- mma_qu$SHP + (2*mma_sdqu$SHP)
        mmam_qu$SHP      <- mma_qu$SHP - (2*mma_sdqu$SHP)
      }
      if("SUBS" %in% colnames(mma_bi)){
        mmap_bi$SUBS      <- mma_bi$SUBS + (2*mma_sdbi$SUBS)
        mmam_bi$SUBS      <- mma_bi$SUBS - (2*mma_sdbi$SUBS)
        mmap_bu$SUBS      <- mma_bu$SUBS + (2*mma_sdbu$SUBS)
        mmam_bu$SUBS      <- mma_bu$SUBS - (2*mma_sdbu$SUBS)
        mmap_qu$SUBS      <- mma_qu$SUBS + (2*mma_sdqu$SUBS)
        mmam_qu$SUBS      <- mma_qu$SUBS - (2*mma_sdqu$SUBS)
      }
      if("SUP" %in% colnames(mma_bi)){
        mmap_bi$SUP      <- mma_bi$SUP + (2*mma_sdbi$SUP)
        mmam_bi$SUP      <- mma_bi$SUP - (2*mma_sdbi$SUP)
        mmap_bu$SUP      <- mma_bu$SUP + (2*mma_sdbu$SUP)
        mmam_bu$SUP      <- mma_bu$SUP - (2*mma_sdbu$SUP)
        mmap_qu$SUP      <- mma_qu$SUP + (2*mma_sdqu$SUP)
        mmam_qu$SUP      <- mma_qu$SUP - (2*mma_sdqu$SUP)
      }
      if("SUPC" %in% colnames(mma_bi)){
        mmap_bi$SUPC      <- mma_bi$SUPC + (2*mma_sdbi$SUPC)
        mmam_bi$SUPC      <- mma_bi$SUPC - (2*mma_sdbi$SUPC)
        mmap_bu$SUPC      <- mma_bu$SUPC + (2*mma_sdbu$SUPC)
        mmam_bu$SUPC      <- mma_bu$SUPC - (2*mma_sdbu$SUPC)
        mmap_qu$SUPC      <- mma_qu$SUPC + (2*mma_sdqu$SUPC)
        mmam_qu$SUPC      <- mma_qu$SUPC - (2*mma_sdqu$SUPC)
      }
      if("TM" %in% colnames(mma_bi)){
        mmap_bi$TM      <- mma_bi$TM + (2*mma_sdbi$TM)
        mmam_bi$TM      <- mma_bi$TM - (2*mma_sdbi$TM)
        mmap_bu$TM      <- mma_bu$TM + (2*mma_sdbu$TM)
        mmam_bu$TM      <- mma_bu$TM - (2*mma_sdbu$TM)
        mmap_qu$TM      <- mma_qu$TM + (2*mma_sdqu$TM)
        mmam_qu$TM      <- mma_qu$TM - (2*mma_sdqu$TM)
      }
      if("TR.S" %in% colnames(mma_bi)){
        mmap_bi$TR.S      <- mma_bi$TR.S + (2*mma_sdbi$TR.S)
        mmam_bi$TR.S      <- mma_bi$TR.S - (2*mma_sdbi$TR.S)
        mmap_bu$TR.S      <- mma_bu$TR.S + (2*mma_sdbu$TR.S)
        mmam_bu$TR.S      <- mma_bu$TR.S - (2*mma_sdbu$TR.S)
        mmap_qu$TR.S      <- mma_qu$TR.S + (2*mma_sdqu$TR.S)
        mmam_qu$TR.S      <- mma_qu$TR.S - (2*mma_sdqu$TR.S)
      }
      if("TR.C" %in% colnames(mma_bi)){
        mmap_bi$TR.C      <- mma_bi$TR.C + (2*mma_sdbi$TR.C)
        mmam_bi$TR.C      <- mma_bi$TR.C - (2*mma_sdbi$TR.C)
        mmap_bu$TR.C      <- mma_bu$TR.C + (2*mma_sdbu$TR.C)
        mmam_bu$TR.C      <- mma_bu$TR.C - (2*mma_sdbu$TR.C)
        mmam_qu$TR.C      <- mma_qu$TR.C + (2*mma_sdqu$TR.C)
        mmap_qu$TR.C      <- mma_qu$TR.C - (2*mma_sdqu$TR.C)
      }
      if("TR.M" %in% colnames(mma_bi)){
        mmap_bi$TR.M      <- mma_bi$TR.M + (2*mma_sdbi$TR.M)
        mmam_bi$TR.M      <- mma_bi$TR.M - (2*mma_sdbi$TR.M)
        mmap_bu$TR.M      <- mma_bu$TR.M + (2*mma_sdbu$TR.M)
        mmam_bu$TR.M      <- mma_bu$TR.M - (2*mma_sdbu$TR.M)
        mmap_qu$TR.M      <- mma_qu$TR.M + (2*mma_sdqu$TR.M)
        mmam_qu$TR.M      <- mma_qu$TR.M - (2*mma_sdqu$TR.M)
      }
      if("TR.L" %in% colnames(mma_bi)){
        mmap_bi$TR.L      <- mma_bi$TR.L + (2*mma_sdbi$TR.L)
        mmam_bi$TR.L      <- mma_bi$TR.L - (2*mma_sdbi$TR.L)
        mmap_bu$TR.L      <- mma_bu$TR.L + (2*mma_sdbu$TR.L)
        mmam_bu$TR.L      <- mma_bu$TR.L - (2*mma_sdbu$TR.L)
        mmap_qu$TR.L      <- mma_qu$TR.L + (2*mma_sdqu$TR.L)
        mmam_qu$TR.L      <- mma_qu$TR.L - (2*mma_sdqu$TR.L)
      }
      if("ADD" %in% colnames(mma_bi)){
        mmap_bi$ADD      <- mma_bi$ADD + (2*mma_sdbi$ADD)
        mmam_bi$ADD      <- mma_bi$ADD - (2*mma_sdbi$ADD)
        mmap_bu$ADD      <- mma_bu$ADD + (2*mma_sdbu$ADD)
        mmam_bu$ADD      <- mma_bu$ADD - (2*mma_sdbu$ADD)
        mmap_qu$ADD      <- mma_qu$ADD + (2*mma_sdqu$ADD)
        mmam_qu$ADD      <- mma_qu$ADD - (2*mma_sdqu$ADD)
      }
      if("AMB" %in% colnames(mma_bi)){
        mmap_bi$AMB      <- mma_bi$AMB + (2*mma_sdbi$AMB)
        mmam_bi$AMB      <- mma_bi$AMB - (2*mma_sdbi$AMB)
        mmap_bu$AMB      <- mma_bu$AMB + (2*mma_sdbu$AMB)
        mmam_bu$AMB      <- mma_bu$AMB - (2*mma_sdbu$AMB)
        mmap_qu$AMB      <- mma_qu$AMB + (2*mma_sdqu$AMB)
        mmam_qu$AMB      <- mma_qu$AMB - (2*mma_sdqu$AMB)
      }
      if("CFB" %in% colnames(mma_bi)){
        mmap_bi$CFB      <- mma_bi$CFB + (2*mma_sdbi$CFB)
        mmam_bi$CFB      <- mma_bi$CFB - (2*mma_sdbi$CFB)
        mmap_bu$CFB      <- mma_bu$CFB + (2*mma_sdbu$CFB)
        mmam_bu$CFB      <- mma_bu$CFB - (2*mma_sdbu$CFB)
        mmap_qu$CFB      <- mma_qu$CFB + (2*mma_sdqu$CFB)
        mmam_qu$CFB      <- mma_qu$CFB - (2*mma_sdqu$CFB)
      }
      if("FMTE" %in% colnames(mma_bi)){
        mmap_bi$FMTE      <- mma_bi$FMTE + (2*mma_sdbi$FMTE)
        mmam_bi$FMTE      <- mma_bi$FMTE - (2*mma_sdbi$FMTE)
        mmap_bu$FMTE      <- mma_bu$FMTE + (2*mma_sdbu$FMTE)
        mmam_bu$FMTE      <- mma_bu$FMTE - (2*mma_sdbu$FMTE)
        mmap_qu$FMTE      <- mma_qu$FMTE + (2*mma_sdqu$FMTE)
        mmam_qu$FMTE      <- mma_qu$FMTE - (2*mma_sdqu$FMTE)
      }
      if("FMTI" %in% colnames(mma_bi)){
        mmap_bi$FMTI      <- mma_bi$FMTI + (2*mma_sdbi$FMTI)
        mmam_bi$FMTI      <- mma_bi$FMTI - (2*mma_sdbi$FMTI)
        mmap_bu$FMTI      <- mma_bu$FMTI + (2*mma_sdbu$FMTI)
        mmam_bu$FMTI      <- mma_bu$FMTI - (2*mma_sdbu$FMTI)
        mmap_qu$FMTI      <- mma_qu$FMTI + (2*mma_sdqu$FMTI)
        mmam_qu$FMTI      <- mma_qu$FMTI - (2*mma_sdqu$FMTI)
      }
      if("FTE" %in% colnames(mma_bi)){
        mmap_bi$FTE      <- mma_bi$FTE + (2*mma_sdbi$FTE)
        mmam_bi$FTE      <- mma_bi$FTE - (2*mma_sdbi$FTE)
        mmap_bu$FTE      <- mma_bu$FTE + (2*mma_sdbu$FTE)
        mmam_bu$FTE      <- mma_bu$FTE - (2*mma_sdbu$FTE)
        mmap_qu$FTE      <- mma_qu$FTE + (2*mma_sdqu$FTE)
        mmam_qu$FTE      <- mma_qu$FTE - (2*mma_sdqu$FTE)
      }
      if("FTI.Fan" %in% colnames(mma_bi)){
        mmap_bi$FTI.Fan      <- mma_bi$FTI.Fan + (2*mma_sdbi$FTI.Fan)
        mmam_bi$FTI.Fan      <- mma_bi$FTI.Fan - (2*mma_sdbi$FTI.Fan)
        mmap_bu$FTI.Fan      <- mma_bu$FTI.Fan + (2*mma_sdbu$FTI.Fan)
        mmam_bu$FTI.Fan      <- mma_bu$FTI.Fan - (2*mma_sdbu$FTI.Fan)
        mmap_qu$FTI.Fan      <- mma_qu$FTI.Fan + (2*mma_sdqu$FTI.Fan)
        mmam_qu$FTI.Fan      <- mma_qu$FTI.Fan - (2*mma_sdqu$FTI.Fan)
      }
      if("FTI2" %in% colnames(mma_bi)){
        mmap_bi$FTI2      <- mma_bi$FTI2 + (2*mma_sdbi$FTI2)
        mmam_bi$FTI2      <- mma_bi$FTI2 - (2*mma_sdbi$FTI2)
        mmap_bu$FTI2      <- mma_bu$FTI2 + (2*mma_sdbu$FTI2)
        mmam_bu$FTI2      <- mma_bu$FTI2 - (2*mma_sdbu$FTI2)
        mmap_qu$FTI2      <- mma_qu$FTI2 + (2*mma_sdqu$FTI2)
        mmam_qu$FTI2      <- mma_qu$FTI2 - (2*mma_sdqu$FTI2)
      }
      if("IFM" %in% colnames(mma_bi)){
        mmap_bi$IFM      <- mma_bi$IFM + (2*mma_sdbi$IFM)
        mmam_bi$IFM      <- mma_bi$IFM - (2*mma_sdbi$IFM)
        mmap_bu$IFM      <- mma_bu$IFM + (2*mma_sdbu$IFM)
        mmam_bu$IFM      <- mma_bu$IFM - (2*mma_sdbu$IFM)
        mmap_qu$IFM      <- mma_qu$IFM + (2*mma_sdqu$IFM)
        mmam_qu$IFM      <- mma_qu$IFM - (2*mma_sdqu$IFM)
      }
      if("ILF" %in% colnames(mma_bi)){
        mmap_bi$ILF      <- mma_bi$ILF + (2*mma_sdbi$ILF)
        mmam_bi$ILF      <- mma_bi$ILF - (2*mma_sdbi$ILF)
        mmap_bu$ILF      <- mma_bu$ILF + (2*mma_sdbu$ILF)
        mmam_bu$ILF      <- mma_bu$ILF - (2*mma_sdbu$ILF)
        mmap_qu$ILF      <- mma_qu$ILF + (2*mma_sdqu$ILF)
        mmam_qu$ILF      <- mma_qu$ILF - (2*mma_sdqu$ILF)
      }
      if("ITB" %in% colnames(mma_bi)){
        mmap_bi$ITB      <- mma_bi$ITB + (2*mma_sdbi$ITB)
        mmam_bi$ITB      <- mma_bi$ITB - (2*mma_sdbi$ITB)
        mmap_bu$ITB      <- mma_bu$ITB + (2*mma_sdbu$ITB)
        mmam_bu$ITB      <- mma_bu$ITB - (2*mma_sdbu$ITB)
        mmap_qu$ITB      <- mma_qu$ITB + (2*mma_sdqu$ITB)
        mmam_qu$ITB      <- mma_qu$ITB - (2*mma_sdqu$ITB)
      }
      if("PIFE" %in% colnames(mma_bi)){
        mmap_bi$PIFE      <- mma_bi$PIFE + (2*mma_sdbi$PIFE)
        mmam_bi$PIFE      <- mma_bi$PIFE - (2*mma_sdbi$PIFE)
        mmap_bu$PIFE      <- mma_bu$PIFE + (2*mma_sdbu$PIFE)
        mmam_bu$PIFE      <- mma_bu$PIFE - (2*mma_sdbu$PIFE)
        mmap_qu$PIFE      <- mma_qu$PIFE + (2*mma_sdqu$PIFE)
        mmam_qu$PIFE      <- mma_qu$PIFE - (2*mma_sdqu$PIFE)
      }
      if("PIFI" %in% colnames(mma_bi)){
        mmap_bi$PIFI      <- mma_bi$PIFI + (2*mma_sdbi$PIFI)
        mmam_bi$PIFI      <- mma_bi$PIFI - (2*mma_sdbi$PIFI)
        mmap_bu$PIFI      <- mma_bu$PIFI + (2*mma_sdbu$PIFI)
        mmam_bu$PIFI      <- mma_bu$PIFI - (2*mma_sdbu$PIFI)
        mmap_qu$PIFI      <- mma_qu$PIFI + (2*mma_sdqu$PIFI)
        mmam_qu$PIFI      <- mma_qu$PIFI - (2*mma_sdqu$PIFI)
      }
      if("SUM" %in% colnames(mma_bi)){
        mmap_bi$SUM      <- 0
        mmam_bi$SUM      <- 0
        mmap_bu$SUM      <- 0
        mmam_bu$SUM      <- 0
        mmap_qu$SUM      <- 0
        mmam_qu$SUM      <- 0
        mma_bi$SUM       <- 0
        mma_bu$SUM       <- 0
        mma_qu$SUM       <- 0
      }
    }
    
    #create opposite directions
    {
      pos_bi<- mma_bi
      pos_bi[pos_bi<0] <- 0
      pos_bi_p<- mmap_bi
      pos_bi_p[pos_bi_p<0] <- 0
      pos_bi_m<- mmam_bi
      pos_bi_m[pos_bi_m<0] <- 0
      neg_bi<- mma_bi
      neg_bi[neg_bi>0] <- 0
      neg_bi<-abs(neg_bi)
      neg_bi$SUM <- 0
      neg_bi_p<- mmap_bi
      neg_bi_p[neg_bi_p>0] <- 0
      neg_bi_p<-abs(neg_bi_p)
      neg_bi_p$SUM <- 0
      neg_bi_m<- mmam_bi
      neg_bi_m[neg_bi_m>0] <- 0
      neg_bi_m<-abs(neg_bi_m)
      neg_bi_m$SUM <- 0
      
      pos_bu<- mma_bu
      pos_bu[pos_bu<0] <- 0
      pos_bu_p<- mmap_bu
      pos_bu_p[pos_bu_p<0] <- 0
      pos_bu_m<- mmam_bu
      pos_bu_m[pos_bu_m<0] <- 0
      neg_bu<- mma_bu
      neg_bu[neg_bu>0] <- 0
      neg_bu<-abs(neg_bu)
      neg_bu$SUM <- 0
      neg_bu_p<- mmap_bu
      neg_bu_p[neg_bu_p>0] <- 0
      neg_bu_p<-abs(neg_bu_p)
      neg_bu_p$SUM <- 0
      neg_bu_m<- mmam_bu
      neg_bu_m[neg_bu_m>0] <- 0
      neg_bu_m<-abs(neg_bu_m)
      neg_bu_m$SUM <- 0
      
      pos_qu<- mma_qu
      pos_qu[pos_qu<0] <- 0
      pos_qu_p<- mmap_qu
      pos_qu_p[pos_qu_p<0] <- 0
      pos_qu_m<- mmam_qu
      pos_qu_m[pos_qu_m<0] <- 0
      neg_qu<- mma_qu
      neg_qu[pos_qu>0] <- 0
      neg_qu<-abs(neg_qu)
      neg_qu$SUM <- 0
      neg_qu_p<- mmap_qu
      neg_qu_p[pos_qu_p>0] <- 0
      neg_qu_p<-abs(neg_qu_p)
      neg_qu_p$SUM <- 0
      neg_qu_m<- mmam_qu
      neg_qu_m[pos_qu_m>0] <- 0
      neg_qu_m<-abs(neg_qu_m)
      neg_qu_m$SUM <- 0
    }
    
    # new sums
    {
      pos_bi$SUM <- rowSums(pos_bi[3:(length(pos_bi))])
      pos_bi_p$SUM <- rowSums(pos_bi_p[3:(length(pos_bi_p)-1)])
      pos_bi_m$SUM <- rowSums(pos_bi_m[3:(length(pos_bi_m)-1)])
      neg_bi$SUM <- rowSums(neg_bi[3:(length(neg_bi))])
      neg_bi_p$SUM <- rowSums(neg_bi_p[3:(length(neg_bi_p)-1)])
      neg_bi_m$SUM <- rowSums(neg_bi_m[3:(length(neg_bi_m)-1)])
      
      pos_bu$SUM <- rowSums(pos_bu[3:(length(pos_bu))])
      pos_bu_p$SUM <- rowSums(pos_bu_p[3:(length(pos_bu_p)-1)])
      pos_bu_m$SUM <- rowSums(pos_bu_m[3:(length(pos_bu_m)-1)])
      neg_bu$SUM <- rowSums(neg_bu[3:(length(neg_bu))])
      neg_bu_p$SUM <- rowSums(neg_bu_p[3:(length(neg_bu_p)-1)])
      neg_bu_m$SUM <- rowSums(neg_bu_m[3:(length(neg_bu_m)-1)])
      
      pos_qu$SUM <- rowSums(pos_qu[3:(length(pos_qu))])
      pos_qu_p$SUM <- rowSums(pos_qu_p[3:(length(pos_qu_p)-1)])
      pos_qu_m$SUM <- rowSums(pos_qu_m[3:(length(pos_qu_m)-1)])
      neg_qu$SUM <- rowSums(neg_qu[3:(length(neg_qu))])
      neg_qu_p$SUM <- rowSums(neg_qu_p[3:(length(neg_qu_p)-1)])
      neg_qu_m$SUM <- rowSums(neg_qu_m[3:(length(neg_qu_m)-1)])
    }
  }
  
  ##Plotting##
  exname <- "Moment Arm Figures/New/"
  zero_line        <- head(zero$SUM,-4)
  data1 <- data.frame (pos_bu$time,pos_bu$SUM,pos_bu_p$SUM,pos_bu_m$SUM,
                       pos_bi$SUM,pos_bi_p$SUM,pos_bi_m$SUM,
                       pos_qu$SUM,pos_qu_p$SUM,pos_qu_m$SUM)
  data2 <- data.frame (pos_bu$time,neg_bu$SUM,neg_bu_p$SUM,neg_bu_m$SUM,
                       neg_bi$SUM,neg_bi_p$SUM,neg_bi_m$SUM,
                       neg_qu$SUM,neg_qu_p$SUM,neg_qu_m$SUM)
  #remove everything after launch
  data1 <- head(data1,-4)
  data2 <- head(data2,-4)
  
  ggplot(data1,aes(x=pos_bu.time)) +
    geom_line(aes(y=zero_line,colour="black"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_bu_m.SUM,ymax=pos_bu_p.SUM),colour="#D81B60",fill="#D81B60",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_bu.SUM,colour="Burst"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_bi_m.SUM,ymax=pos_bi_p.SUM),colour="#1E88E5",fill="#1E88E5",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_bi.SUM,colour="Counter"), linewidth=1) +
    geom_ribbon(aes(ymin=pos_qu_m.SUM,ymax=pos_qu_p.SUM),colour="#FFC107",fill="#FFC107",linetype=2,alpha=0.25) +
    geom_line(aes(y=pos_qu.SUM,colour="Quad"), linewidth=1) +
    labs(x="Time", y="Summed Moment Arms") +
    ggtitle(paste0(ti))+
    theme(plot.title=element_text(hjust=0.5))+
    scale_fill_identity(name='Legend',guide='legend',labels=c('Burst','Counter','Quad')) +
    scale_color_manual(name='Launch',values=c('Burst'='#D81B60','Counter'='#1E88E5','Quad'='#FFC107'), 
                       labels=c('Burst','Counter','Quad'))
  #ggsave((paste0(exname,ti," Moment Arms.pdf")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  #ggsave((paste0(exname,ti," Moment Arms.png")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  write.csv(data1,paste0(exname,ti," values.csv"))
  
  ggplot(data2,aes(x=pos_bu.time)) +
    geom_line(aes(y=zero_line,colour="black"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_bu_m.SUM,ymax=neg_bu_p.SUM),colour="#D81B60",fill="#D81B60",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_bu.SUM,colour="Burst"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_bi_m.SUM,ymax=neg_bi_p.SUM),colour="#1E88E5",fill="#1E88E5",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_bi.SUM,colour="Counter"), linewidth=1) +
    geom_ribbon(aes(ymin=neg_qu_m.SUM,ymax=neg_qu_p.SUM),colour="#FFC107",fill="#FFC107",linetype=2,alpha=0.25) +
    geom_line(aes(y=neg_qu.SUM,colour="Quad"), linewidth=1) +
    labs(x="Time", y="Summed Moment Arms") +
    ggtitle(paste0(ti2))+
    theme(plot.title=element_text(hjust=0.5))+
    scale_fill_identity(name='Legend',guide='legend',labels=c('Burst','Counter','Quad'))+
    scale_color_manual(name='Launch',values=c('Burst'='#D81B60','Counter'='#1E88E5','Quad'='#FFC107'), 
                       labels=c('Burst','Counter','Quad'))
  #ggsave((paste0(exname,ti2," Moment Arms.pdf")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  #ggsave((paste0(exname,ti2," Moment Arms.png")), plot = last_plot(), device = NULL, path = NULL, scale = 1, dpi = 600)
  write.csv(data2,paste0(exname,ti2," values.csv"))
}  

graph.loop (j = "Shoulder", df = "_ABAD", ti= "Shoulder Abduction",         ti2="Shoulder Adduction") 
graph.loop (j = "Shoulder", df = "_LAR",  ti= "Shoulder External Rotation", ti2="Shoulder Internal Rotation")   
graph.loop (j = "Shoulder", df = "_FE",   ti= "Shoulder Extension",         ti2= "Shoulder Flexion")
graph.loop (j = "Elbow",    df = "_ABAD", ti= "Elbow Abduction",            ti2="Elbow Adduction") 
graph.loop (j = "Elbow",    df = "_LAR",  ti= "Elbow External Rotation",    ti2="Elbow Internal Rotation")   
graph.loop (j = "Elbow",    df = "_FE",   ti= "Elbow Extension",            ti2="Elbow Flexion")
graph.loop (j = "Wrist",    df = "_ABAD", ti= "Wrist Abduction",            ti2="Wrist Adduction") 
graph.loop (j = "Wrist",    df = "_LAR",  ti= "Wrist External Rotation",    ti2="Wrist Internal Rotation")   
graph.loop (j = "Wrist",    df = "_FE",   ti= "Wrist Extension",            ti2="Wrist Flexion")
graph.loop (j = "WMC",      df = "_ABAD", ti= "WMC Abduction",              ti2="WMC Adduction") 
graph.loop (j = "WMC",      df = "_LAR",  ti= "WMC External Rotation",      ti2="WMC Internal Rotation")   
graph.loop (j = "WMC",      df = "_FE",   ti= "WMC Extension",              ti2="WMC Flexion")
graph.loop (j = "WP1",      df = "_ABAD", ti= "WP1 Abduction",              ti2="WP1 Adduction") 
graph.loop (j = "WP1",      df = "_LAR",  ti= "WP1 External Rotation",      ti2="WP1 Internal Rotation")   
graph.loop (j = "WP1",      df = "_FE",   ti= "WP1 Extension",              ti2="WP1 Flexion")
graph.loop (j = "Hip",      df = "_ABAD", ti= "Hip Abduction",              ti2="Hip Adduction") 
graph.loop (j = "Hip",      df = "_LAR",  ti= "Hip External Rotation",      ti2="Hip Internal Rotation")   
graph.loop (j = "Hip",      df = "_FE",   ti= "Hip Extension",              ti2="Hip Flexion")
graph.loop (j = "Knee",     df = "_ABAD", ti= "Knee Abduction",             ti2="Knee Adduction") 
graph.loop (j = "Knee",     df = "_LAR",  ti= "Knee External Rotation",     ti2="Knee Internal Rotation")   
graph.loop (j = "Knee",     df = "_FE",   ti= "Knee Extension",             ti2="Knee Flexion")

#Bird
graph.loopb (j = "Shoulder", df = "_ABAD", ti= "Shoulder Abduction b",         ti2="Shoulder Adduction b") 
graph.loopb (j = "Shoulder", df = "_LAR",  ti= "Shoulder External Rotation b", ti2="Shoulder Internal Rotation b")   
graph.loopb (j = "Shoulder", df = "_FE",   ti= "Shoulder Extension b",         ti2= "Shoulder Flexion b")
graph.loopb (j = "Elbow",    df = "_ABAD", ti= "Elbow Abduction b",            ti2="Elbow Adduction b") 
graph.loopb (j = "Elbow",    df = "_LAR",  ti= "Elbow External Rotation b",    ti2="Elbow Internal Rotation b")   
graph.loopb (j = "Elbow",    df = "_FE",   ti= "Elbow Extension b",            ti2="Elbow Flexion b")
graph.loopb (j = "Wrist",    df = "_ABAD", ti= "Wrist Abduction b",            ti2="Wrist Adduction b") 
graph.loopb (j = "Wrist",    df = "_LAR",  ti= "Wrist External Rotation b",    ti2="Wrist Internal Rotation b")   
graph.loopb (j = "Wrist",    df = "_FE",   ti= "Wrist Extension b",            ti2="Wrist Flexion b")
graph.loopb (j = "WMC",      df = "_ABAD", ti= "WMC Abduction b",              ti2="WMC Adduction b") 
graph.loopb (j = "WMC",      df = "_LAR",  ti= "WMC External Rotation b",      ti2="WMC Internal Rotation b")   
graph.loopb (j = "WMC",      df = "_FE",   ti= "WMC Extension b",              ti2="WMC Flexion b")
graph.loopb (j = "WP1",      df = "_ABAD", ti= "WP1 Abduction b",              ti2="WP1 Adduction b") 
graph.loopb (j = "WP1",      df = "_LAR",  ti= "WP1 External Rotation b",      ti2="WP1 Internal Rotation b")   
graph.loopb (j = "WP1",      df = "_FE",   ti= "WP1 Extension b",              ti2="WP1 Flexion b")
graph.loopb (j = "Hip",      df = "_ABAD", ti= "Hip Abduction b",              ti2="Hip Adduction b") 
graph.loopb (j = "Hip",      df = "_LAR",  ti= "Hip External Rotation b",      ti2="Hip Internal Rotation b")   
graph.loopb (j = "Hip",      df = "_FE",   ti= "Hip Extension b",              ti2="Hip Flexion b")
graph.loopb (j = "Knee",     df = "_ABAD", ti= "Knee Abduction b",             ti2="Knee Adduction b") 
graph.loopb (j = "Knee",     df = "_LAR",  ti= "Knee External Rotation b",     ti2="Knee Internal Rotation b")   
graph.loopb (j = "Knee",     df = "_FE",   ti= "Knee Extension b",             ti2="Knee Flexion b")

#croc
graph.loopc (j = "Shoulder", df = "_ABAD", ti= "Shoulder Abduction c",         ti2="Shoulder Adduction c") 
graph.loopc (j = "Shoulder", df = "_LAR",  ti= "Shoulder External Rotation c", ti2="Shoulder Internal Rotation c")   
graph.loopc (j = "Shoulder", df = "_FE",   ti= "Shoulder Extension c",         ti2= "Shoulder Flexion c")
graph.loopc (j = "Elbow",    df = "_ABAD", ti= "Elbow Abduction c",            ti2="Elbow Adduction c") 
graph.loopc (j = "Elbow",    df = "_LAR",  ti= "Elbow External Rotation c",    ti2="Elbow Internal Rotation c")   
graph.loopc (j = "Elbow",    df = "_FE",   ti= "Elbow Extension c",            ti2="Elbow Flexion c")
graph.loopc (j = "Wrist",    df = "_ABAD", ti= "Wrist Abduction c",            ti2="Wrist Adduction c") 
graph.loopc (j = "Wrist",    df = "_LAR",  ti= "Wrist External Rotation c",    ti2="Wrist Internal Rotation c")   
graph.loopc (j = "Wrist",    df = "_FE",   ti= "Wrist Extension c",            ti2="Wrist Flexion c")
graph.loopc (j = "WMC",      df = "_ABAD", ti= "WMC Abduction c",              ti2="WMC Adduction c") 
graph.loopc (j = "WMC",      df = "_LAR",  ti= "WMC External Rotation c",      ti2="WMC Internal Rotation c")   
graph.loopc (j = "WMC",      df = "_FE",   ti= "WMC Extension c",              ti2="WMC Flexion c")
graph.loopc (j = "WP1",      df = "_ABAD", ti= "WP1 Abduction c",              ti2="WP1 Adduction c") 
graph.loopc (j = "WP1",      df = "_LAR",  ti= "WP1 External Rotation c",      ti2="WP1 Internal Rotation c")   
graph.loopc (j = "WP1",      df = "_FE",   ti= "WP1 Extension c",              ti2="WP1 Flexion c")
graph.loopc (j = "Hip",      df = "_ABAD", ti= "Hip Abduction c",              ti2="Hip Adduction c") 
graph.loopc (j = "Hip",      df = "_LAR",  ti= "Hip External Rotation c",      ti2="Hip Internal Rotation c")   
graph.loopc (j = "Hip",      df = "_FE",   ti= "Hip Extension c",              ti2="Hip Flexion c")
graph.loopc (j = "Knee",     df = "_ABAD", ti= "Knee Abduction c",             ti2="Knee Adduction c") 
graph.loopc (j = "Knee",     df = "_LAR",  ti= "Knee External Rotation c",     ti2="Knee Internal Rotation c")   
graph.loopc (j = "Knee",     df = "_FE",   ti= "Knee Extension c",             ti2="Knee Flexion c")