[chronojump] raceAnalyzer. Start testing with a new getArea function

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

commit 910eccda953cafc11f79aeb2ed8656a244decf3f
Author: Xavier Padullés <x padulles gmail com>
Date:   Fri Mar 22 12:36:33 2019 +0100

    raceAnalyzer. Start testing with a new getArea function

 r-scripts/scripts-util.R  | 38 +++++++++++++++++++++++++++++++-------
 r-scripts/sprintEncoder.R | 47 +++++++++++++++++++++++++----------------------
 2 files changed, 56 insertions(+), 29 deletions(-)
---
diff --git a/r-scripts/scripts-util.R b/r-scripts/scripts-util.R
index c5fe6b52..1485b1e4 100644
--- a/r-scripts/scripts-util.R
+++ b/r-scripts/scripts-util.R
@@ -1,14 +1,14 @@
 #Function to get the interpolated x at a given y
 interpolateXAtY <- function(X, Y, desiredY){
         if(max(Y) < desiredY){
-                print("desiredY is greater than max(Y)")
+                # print("desiredY is greater than max(Y)")
                 return(max(Y))
         }
                 
         #find the closest sample
         nextSample = 1
-        print("Calculating the interpolation")
-        print(paste("desiredY:", desiredY))
+        # print("Calculating the interpolation")
+        # print(paste("desiredY:", desiredY))
         while (Y[nextSample] < desiredY){
                 nextSample = nextSample +1
         }
@@ -20,7 +20,7 @@ interpolateXAtY <- function(X, Y, desiredY){
         } else {
                 correspondingX = X[previousSample] + (desiredY  - Y[previousSample]) * (X[nextSample] - 
X[previousSample]) / (Y[nextSample] - Y[previousSample])
         }
-        print(paste("CorrespondingX:", correspondingX))
+        # print(paste("CorrespondingX:", correspondingX))
         return(correspondingX)
 }
 
@@ -50,18 +50,41 @@ getAreaUnderCurve <- function(x, y)
         return(totalArea/2) #The area of the parallelograms are twice the triangles areas
 }
 
+getAreaUnderCurve2 <- function(x, y)
+{
+        print("Calculating Area")
+        # print("X:")
+        # print(x)
+        # print("Y:")
+        # print(y)
+
+        totalArea = 0
+        # print(paste("V",1," = ", "(",x[1 + 1] - x[1],",", y[1 + 1] - y[1], ")", sep = ""))
+        for(i in 1:length(x))
+        {
+                barArea = y[i]*(x[i+1] - x[i])
+                
+                # print(paste("V",i," = ", "(",x[i + 1] - x[1],",", y[i + 1] - y[1], ")", sep = ""))
+                # print(parallelogramArea/2)
+                
+                totalArea = totalArea + barArea
+        }
+        # 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, "]"))
+        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")
+        #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")
+        #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)]
@@ -72,6 +95,7 @@ getMeanValue <- function(X, Y, startX, endX)
         
         #calculating the area under the curve (integral)
         area = getAreaUnderCurve(X , Y)
+
         return(area / (X[length(X)] - X[1]))
 }
 
diff --git a/r-scripts/sprintEncoder.R b/r-scripts/sprintEncoder.R
index de371942..537783bb 100644
--- a/r-scripts/sprintEncoder.R
+++ b/r-scripts/sprintEncoder.R
@@ -60,23 +60,23 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
         Af = 0.2025*(Height^0.725)*(Mass^0.425)*0.266 # Model of the frontal area
         Ka = 0.5*ro*Af*Cd
         
-        encoderCarrera = read.csv2(file = filename, sep = ";")
-        colnames(encoderCarrera) = c("displacement", "time", "force")
-        encoderCarrera$force = encoderCarrera$force * 0.140142 /2  #TODO: Implement the calibration factor 
comming from the arduino
-        totalTime = encoderCarrera$time/1E6     #Converting microseconds to seconds
+        raceAnalyzer = read.csv2(file = filename, sep = ";")
+        colnames(raceAnalyzer) = c("displacement", "time", "force")
+        raceAnalyzer$force = raceAnalyzer$force * 0.140142 /2  #TODO: Implement the calibration factor 
comming from the arduino
+        totalTime = raceAnalyzer$time/1E6     #Converting microseconds to seconds
         elapsedTime = diff(c(0,totalTime))      #The elapsed time between each sample
         
         #The encoder can be used with both loose ends of the rope. This means that the encoder can rotate
         #in both directions. We considre always the speed as positive.
-        if(mean(encoderCarrera$displacement) < 0)
-                encoderCarrera$displacement = -encoderCarrera$displacement
+        if(mean(raceAnalyzer$displacement) < 0)
+                raceAnalyzer$displacement = -raceAnalyzer$displacement
         diameter = 0.16         #Diameter of the pulley
-        #encoderCarrera$displacement = encoderCarrera$displacement * 2 * pi * diameter / 200
+        #raceAnalyzer$displacement = raceAnalyzer$displacement * 2 * pi * diameter / 200
         #In 30m there are 12267 pulses.
         #TODO: measure this several times to have an accurate value
-        encoderCarrera$displacement = encoderCarrera$displacement * 30 / 12267
-        position = cumsum(encoderCarrera$displacement)
-        speed = encoderCarrera$displacement / elapsedTime
+        raceAnalyzer$displacement = raceAnalyzer$displacement * 30 / 12267
+        position = cumsum(raceAnalyzer$displacement)
+        speed = raceAnalyzer$displacement / elapsedTime
         accel = (speed[3] - speed[1]) / (totalTime[3] - totalTime[1])
         for(i in 3:(length(speed) -1))
         {
@@ -84,7 +84,7 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
         }
         accel = c(accel[1],accel, accel[length(accel)])
         forceBody = accel * Mass + Ka*(speed - Vw)^2
-        totalForce = forceBody + encoderCarrera$force
+        totalForce = forceBody + raceAnalyzer$force
         power = totalForce * speed
         # print("position:")
         # print(position)
@@ -92,10 +92,10 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
         # print(speed)
         # print("accel:")
         # print(accel)
-        # print("forceBody:")
-        # print(forceBody)
-        # print("forceRope:")
-        # print(encoderCarrera$force)
+        print("forceBody:")
+        print(forceBody)
+        print("forceRope:")
+        print(raceAnalyzer$force)
         # print("power:")
         # print(power)
         
@@ -135,7 +135,7 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
         #Zeroing position to the initial acceleration sample
         position = position - position[trimmingSamples$start]
         data = data.frame(time = time[trimmingSamples$start:trimmingSamples$end], speed = 
speed[trimmingSamples$start:trimmingSamples$end])
-        print(data)
+        #print(data)
         
         print("Trying nls")
         regression = tryNLS(data)
@@ -174,7 +174,7 @@ getSprintFromEncoder <- function(filename, testLength, Mass, Temperature = 25, H
 
 plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title = "Test graph",
                                   plotRawMeanSpeed = TRUE, plotRawSpeed = TRUE, plotRawAccel = FALSE, 
plotRawForce = FALSE, plotMeanRawForce = TRUE, plotRawPower = FALSE, plotMeanRawPower = TRUE,
-                                  plotFittedSpeed = TRUE, plotFittedAccel = FALSE, plotFittedForce = FALSE, 
plotFittedPower = FALSE)
+                                  plotFittedSpeed = TRUE, plotFittedAccel = FALSE, plotFittedForce = FALSE, 
plotFittedPower = TRUE)
 {
         #Plotting position
         # plot(sprintRawDynamics$time[sprintRawDynamics$startSample:sprintRawDynamics$endSample], 
sprintRawDynamics$rawPosition[sprintRawDynamics$startSample:sprintRawDynamics$endSample],
@@ -315,6 +315,7 @@ plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title
         if(plotRawForce|| plotFittedForce)
         {
                 ylimits = 
c(min(sprintRawDynamics$rawForce[sprintRawDynamics$startSample:sprintRawDynamics$endSample])*0.95, 
max(c(sprintRawDynamics$rawFmax, sprintFittedDynamics$fmax.fitted)*1.05))
+                #ylimits = c(0,sprintFittedDynamics$fmax.fitted)
                 if (plotRawForce)
                 {
                         #Plotting rawForce
@@ -339,6 +340,8 @@ plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title
                         legendColor = c(legendColor, "blue")
                 }
                 axis(side = 4, col = "blue", line = 2)
+                print("Mean force from the model")
+                print(getMeanValue(sprintFittedDynamics$t.fitted, sprintRawDynamics$force.fitted, 0, 1.004))
         }     
 
         if(plotMeanRawForce)
@@ -357,7 +360,7 @@ plotSprintFromEncoder <- function(sprintRawDynamics, sprintFittedDynamics, title
                         text(splitTime[n] + (splitTime[n+1] - splitTime[n])*0.2, meanForce[n+1], 
paste(round(meanForce[n+1], digits = 2), "N"), col = "blue")
                         
                 }
-                axis(side = 4, col = "blue", line = 2)
+                #axis(side = 4, col = "blue", line = 2)
         }
         
         if(plotRawPower|| plotFittedPower)
@@ -479,8 +482,8 @@ tryNLS <- function(data){
                 {
                         model = nls(speed ~ Vmax*(1-exp(-K*time)), data,
                                     start = list(Vmax = max(data[,"speed"]), K = 1), 
control=nls.control(warnOnly=TRUE))
-                        print("model:")
-                        print(model)
+                        # print("model:")
+                        # print(model)
                         return(list(regressionDone = TRUE, model = model))
                 }, 
                 error=function(cond)
@@ -494,8 +497,8 @@ tryNLS <- function(data){
 testEncoderCJ <- function(filename, testLength, mass, personHeight, tempC)
 {
         sprintRawDynamics = getSprintFromEncoder(filename, testLength, op$mass, op$tempC, op$personHeight)
-        print("sprintRawDynamics:")
-        print(sprintRawDynamics)
+        # print("sprintRawDynamics:")
+        # print(sprintRawDynamics)
         if (sprintRawDynamics$longEnough & sprintRawDynamics$regressionDone)
         {
                 sprintFittedDynamics = getDynamicsFromSprint(K = sprintRawDynamics$K, Vmax = 
sprintRawDynamics$Vmax, mass, tempC, personHeight)