[chronojump] RaceAnalyzer. Introduced force sensor calculations

[Xavier Padullés <xpadulles src gnome org>]

commit e389e7a39d23cfd3dd865a29712634261fe90796
Author: Xavier Padullés <x padulles gmail com>
Date:   Tue Mar 12 14:17:30 2019 +0100

    RaceAnalyzer. Introduced force sensor calculations

 r-scripts/scripts-util.R  | 66 ++++++++++++++++++++++++++++++++++++++++-------
 r-scripts/sprintEncoder.R | 44 +++++++++++++++++++++++++------
 2 files changed, 93 insertions(+), 17 deletions(-)
---
diff --git a/r-scripts/scripts-util.R b/r-scripts/scripts-util.R
index 05db2963..c94aff55 100644
--- a/r-scripts/scripts-util.R
+++ b/r-scripts/scripts-util.R
@@ -7,12 +7,8 @@ interpolateXAtY <- function(X, Y, desiredY){
                 
         #find the closest sample
         nextSample = 1
-        print("X:")
-        print(X)
-        print("Y:")
-        print(Y)
-        print("desiredY:")
-        print(desiredY)
+        print("Calculating the interpolation")
+        print(paste("desiredY:", desiredY))
         while (Y[nextSample] < desiredY){
                 nextSample = nextSample +1
         }
@@ -20,11 +16,63 @@ interpolateXAtY <- function(X, Y, desiredY){
         previousSample = nextSample - 1
         
         if(Y[nextSample] == desiredY){
-                desiredX = X[nextSample]
+                correspondingX = X[nextSample]
         } else {
-                desiredX = X[previousSample] + (desiredY  - Y[previousSample]) * (X[nextSample] - 
X[previousSample]) / (Y[nextSample] - Y[previousSample])
+                correspondingX = X[previousSample] + (desiredY  - Y[previousSample]) * (X[nextSample] - 
X[previousSample]) / (Y[nextSample] - Y[previousSample])
+        }
+        print(paste("CorrespondingX:", correspondingX))
+        return(correspondingX)
+}
+
+
+#Calculate the area under a curve using the cross product of consecutive vectors with the origin in the 
first point
+getAreaUnderCurve <- function(x, y)
+{
+        print("Calculating Area")
+        print("X:")
+        print(x)
+        print("Y:")
+        print(y)
+        x = c(x, x[length(x)], x[1])
+        y = c(y, 0, 0)
+        totalArea = 0
+        print(paste("V",1," = ", "(",x[1 + 1] - x[1],",", y[1 + 1] - y[1], ")", sep = ""))
+        for(i in 2:(length(x) -1))
+        {
+                parallelogramArea = ((x[i + 1] - x[1])* (y[i] - y[1]) - (x[i] - x[1]) * (y[i+1] - y[1]))
+                
+                # print(paste("V",i," = ", "(",x[i + 1] - x[1],",", y[i + 1] - y[1], ")", sep = ""))
+                # print(parallelogramArea/2)
+                
+                totalArea = totalArea + parallelogramArea
         }
-        return(desiredX)
+        print(paste("toalArea:", totalArea/2))
+        return(totalArea/2) #The area of the parallelograms are twice the triangles areas
+}
+
+#Calculates the mean of a curve interval
+getMeanValue <- function(X, Y, startX, endX)
+{
+        print(paste("Calculating mean In the X range of [", startX, ",", endX, "]"))
+        
+        
+        #Calculating the value of Y corresponding at startX
+        print("Calculating the first Y value")
+        startY = interpolateXAtY(X = Y, Y = X, desiredY = startX) #The order are changed because we are 
looking for the Y value instead of X value
+        
+        #Calculating the last value using the interpolation
+        print("Calculating the last Y value")
+        endY = interpolateXAtY(X = Y, Y = X, desiredY = endX) #The order are changed because we are looking 
for the Y value instead of X value
+
+        X = X[which(X > startX & X < endX)]
+        X = c(startX, X, endX)
+        
+        Y = Y[which(X > startX & X < endX)]
+        Y = c(startY, Y, endY)
+        
+        #calculating the area under the curve (integral)
+        area = getAreaUnderCurve(X , Y)
+        return(area / (X[length(X)] - X[1]))
 }
 
 prepareGraph <- function(os, pngFile, width, height)
diff --git a/r-scripts/sprintEncoder.R b/r-scripts/sprintEncoder.R
index 77dcb48e..5a5b49a6 100644
--- a/r-scripts/sprintEncoder.R
+++ b/r-scripts/sprintEncoder.R
@@ -61,6 +61,7 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
         
         encoderCarrera = read.csv2(file = filename, sep = ";")
         colnames(encoderCarrera) = c("displacement", "time", "force")
+        encoderCarrera$force = encoderCarrera$force * 0.140142 /2
         totalTime = encoderCarrera$time/1E6     #Converting microseconds to seconds
         elapsedTime = diff(c(0,totalTime))      #The elapsed time between each sample
         
@@ -77,13 +78,24 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
         speed = encoderCarrera$displacement / elapsedTime
         #accel = c(0,diff(speed)/elapsedTime[2:length(elapsedTime)])
         accel = (speed[3] - speed[1]) / (totalTime[3] - totalTime[1])
-        for(i in 3:length(speed))
+        for(i in 3:(length(speed) -1))
         {
                 accel = c(accel, (speed[i+1] - speed[i-1]) / (totalTime[i +1] - totalTime[i -1]))
         }
-        accel = c(accel, accel[length(accel)])
-        force = accel * Mass + Ka*(speed - Vw)^2
-        power = force * speed
+        accel = c(accel[1],accel, accel[length(accel)])
+        forceBody = accel * Mass + Ka*(speed - Vw)^2
+        totalForce = forceBody + encoderCarrera$force
+        power = totalForce * speed
+        # print("speed:")
+        # print(speed)
+        # print("accel:")
+        # print(accel)
+        # print("forceBody:")
+        # print(forceBody)
+        # print("forceRope:")
+        # print(encoderCarrera$force)
+        print("power:")
+        print(power)
         
         #Finding when the sprint starts
         trimmingSamples = getTrimmingSamples(totalTime, position, speed, accel, testLength)
@@ -98,14 +110,14 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
         Vmax =summary(model)$coeff[1,1]
         K = summary(model)$coeff[2,1]
         return(list(Vmax = Vmax, K = K,
-                    time = time, rawPosition = position, rawSpeed = speed, rawAccel = accel, rawForce = 
force, rawPower = power,
-                    rawVmax = max(speed[trimmingSamples$start:trimmingSamples$end]), rawAmax = 
max(accel[trimmingSamples$start:trimmingSamples$end]), rawFmax = 
max(force[trimmingSamples$start:trimmingSamples$end]), rawPmax = 
max(power[trimmingSamples$start:trimmingSamples$end]),
+                    time = time, rawPosition = position, rawSpeed = speed, rawAccel = accel, rawForce = 
totalForce, rawPower = power,
+                    rawVmax = max(speed[trimmingSamples$start:trimmingSamples$end]), rawAmax = 
max(accel[trimmingSamples$start:trimmingSamples$end]), rawFmax = 
max(totalForce[trimmingSamples$start:trimmingSamples$end]), rawPmax = 
max(power[trimmingSamples$start:trimmingSamples$end]),
                     startSample = trimmingSamples$start, endSample = trimmingSamples$end, testLength = 
testLength))
 }
 
 plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title = "Test graph",
-                                  plotRawMeanSpeed = TRUE, plotRawSpeed = TRUE, plotRawAccel = TRUE, 
plotRawForce = TRUE, plotRawPower = FALSE,
-                                  plotFittedSpeed = TRUE, plotFittedAccel = FALSE, plotFittedForce = TRUE, 
plotFittedPower = FALSE)
+                                  plotRawMeanSpeed = TRUE, plotRawSpeed = TRUE, plotRawAccel = FALSE, 
plotRawForce = FALSE, plotMeanRawForce = TRUE, plotRawPower = FALSE, plotRawMeanPower = TRUE,
+                                  plotFittedSpeed = TRUE, plotFittedAccel = FALSE, plotFittedForce = FALSE, 
plotFittedPower = FALSE)
 {
         #Plotting position
         # plot(sprintRawDynamics$time[sprintRawDynamics$startSample:sprintRawDynamics$endSample], 
sprintRawDynamics$rawPosition[sprintRawDynamics$startSample:sprintRawDynamics$endSample],
@@ -145,6 +157,9 @@ plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title
                                   
sprintRawDynamics$rawPosition[sprintRawDynamics$startSample:sprintRawDynamics$endSample],
                                   splitPosition)
         meanSpeed = splitPosition / splitTime
+        meanForce =getMeanValue(sprintRawDynamics$time, sprintRawDynamics$rawForce, 
sprintRawDynamics$time[sprintRawDynamics$startSample], splitTime)
+        meanPower =getMeanValue(sprintRawDynamics$time, sprintRawDynamics$rawPower, 
sprintRawDynamics$time[sprintRawDynamics$startSample], splitTime)
+
         while(splitPosition[length(splitPosition)] + 5 < sprintRawDynamics$testLength)
         {
                 splitPosition = c(splitPosition, splitPosition[length(splitPosition)] + 5)
@@ -153,6 +168,10 @@ plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title
                                                      splitPosition[length(splitPosition)]))
                 meanSpeed = c(meanSpeed, (splitPosition[length(splitPosition)] - 
splitPosition[length(splitPosition) -1]) /
                                       (splitTime[length(splitTime)] - splitTime[length(splitTime) -1]))
+                meanForce = c(meanForce, getMeanValue(sprintRawDynamics$time, sprintRawDynamics$rawForce,
+                                                      splitTime[length(splitTime) -1], 
splitTime[length(splitTime)]))
+                meanPower = c(meanPower, getMeanValue(sprintRawDynamics$time, sprintRawDynamics$rawPower,
+                                                      splitTime[length(splitTime) -1], 
splitTime[length(splitTime)]))
         }
         splitPosition = c(splitPosition, sprintRawDynamics$testLength)
         splitTime = c(splitTime, 
interpolateXAtY(sprintRawDynamics$time[sprintRawDynamics$startSample:sprintRawDynamics$endSample],
@@ -161,6 +180,15 @@ plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title
         meanSpeed = c(meanSpeed, (splitPosition[length(splitPosition)] - splitPosition[length(splitPosition) 
-1]) /
                               (splitTime[length(splitTime)] - splitTime[length(splitTime) -1]))
         
+        meanForce = c(meanForce, getMeanValue(sprintRawDynamics$time, sprintRawDynamics$rawForce,
+                                              splitTime[length(splitTime) -1], splitTime[length(splitTime)]))
+        meanPower = c(meanPower, getMeanValue(sprintRawDynamics$time, sprintRawDynamics$rawPower,
+                                              splitTime[length(splitTime) -1], splitTime[length(splitTime)]))
+        print("meanForce:")
+        print(meanForce)
+        print("meanPower:")
+        print(meanPower)
+        
         if(plotRawMeanSpeed)
         {
                 barplot(height = meanSpeed, width = diff(c(0,splitTime)), space = 0,