diff --git a/ShoulderCase/CoordinateSystem.m b/ShoulderCase/CoordinateSystem.m
index b33cb56..90292de 100644
--- a/ShoulderCase/CoordinateSystem.m
+++ b/ShoulderCase/CoordinateSystem.m
@@ -1,209 +1,210 @@
 classdef CoordinateSystem < handle
 
 
 
   properties
     origin
     xAxis
     yAxis
     zAxis
   end
 
 
 
   properties (Hidden = true)
     rotationMatrix
 
     % enable retro-compatibility with codes that uses PA, IS, and ML axes
     PA
     IS
     ML
   end
 
 
 
   methods
 
     function obj = CoordinateSystem()
       obj.origin    = [0 0 0];
       obj.xAxis     = [1 0 0];
       obj.yAxis     = [0 1 0];
       obj.zAxis     = [0 0 1];
 
     end
 
     function set.xAxis(obj,value)
       % enable retro-compatibility with codes that uses a PA axis
       obj.xAxis               = value;
       obj.PA                  = value;
       obj.rotationMatrix(:,1) = value';
     end
 
     function set.yAxis(obj,value)
       % enable retro-compatibility with codes that uses an IS axis
       obj.yAxis               = value;
       obj.IS                  = value;
       obj.rotationMatrix(:,2) = value';
     end
 
     function set.zAxis(obj,value)
       % enable retro-compatibility with codes that uses a ML axis
       obj.zAxis               = value;
       obj.ML                  = value;
       obj.rotationMatrix(:,3) = value';
     end
 
     function points = express(obj,points)
       % Give the coordinates of points in present coordinate system
       if (size(points,2) ~= 3)
         error('Points have to be a Nx3 double array.')
         return;
       end
       points = (points-obj.origin) * obj.rotationMatrix;
     end
 
     function output = isOrthogonal(obj)
       % The axes comparison can't be done looking for strict equality due to
       % rounded values. Therefore the values must be evaluated with a given
       % tolerance
-      tolerance = 10^-9;
-      if dot(obj.xAxis,obj.yAxis) > tolerance
+      minAngle = (pi/2)-(10^-6); % Arbitrarily defined as a deviation of
+                                 % a millionth of a radian
+      if dot(obj.xAxis,obj.yAxis) > norm(obj.xAxis)*norm(obj.yAxis)*cos(minAngle)
         output = false;
         return
       end
-      if dot(obj.xAxis,obj.zAxis) > tolerance
+      if dot(obj.xAxis,obj.zAxis) > norm(obj.xAxis)*norm(obj.zAxis)*cos(minAngle)
         output = false;
         return
       end
-      if dot(obj.zAxis,obj.yAxis) > tolerance
+      if dot(obj.zAxis,obj.yAxis) > norm(obj.zAxis)*norm(obj.yAxis)*cos(minAngle)
         output = false;
         return
       end
       output = true;
     end
 
     function output = isRightHanded(obj)
       % The comparison between a theoretical right-handed axis and the actual
       % z Axis can't be done looking for strict vector equality due to rounded
       % values. Therefore the norm of the sum of the two vectors is evaluated.
       if not(obj.isOrthogonal)
         output = false;
         return
       end
       rightHandedAxis = cross(obj.xAxis,obj.yAxis);
       if norm(rightHandedAxis+obj.zAxis) < norm(rightHandedAxis)
         output = false;
         return
       end
       output = true;
     end
 
     function output = isLeftHanded(obj)
       if not(obj.isOrthogonal)
         output = false;
         return
       end
       if obj.isRightHanded
         output = false;
         return
       end
       output = true;
     end
 
     function setSystemWithScapulaLandmarks(obj,scapula)
       % Calculate the (EPFL) scapular coordinate system
       % The EPFL scapular coordinate system is defined for a right scapula see
       % paper (10.1302/0301-620X.96B4.32641). The X axis is
       % antero-posterior, the Y axis is infero-superior, the Z axis
       % is meddio-lateral. This system is right-handed.
       % This orientation corresponds approximatively to the ISB
       % recommendadtion (doi:10.1016/j.jbiomech.2004.05.042).
       % For left scapula we keep the x axis as postero-anterior and
       % get a left-handed coordinate system.
 
       if (length(scapula.groove) ~= 5)
         error(['Can''t set scapular coordinate system with actual',...
                'scapula.groove']);
         return;
       end
 
       % The four following methods have to be called in this specific order
       obj.setXAxisWithScapulaLandmarks(scapula);
       obj.setZAxisWithScapulaLandmarks(scapula);
       obj.setYAxisWithScapulaLandmarks(scapula);
       obj.setOriginWithScapulaLandmarks(scapula);
     end
 
     function setSystemWith(obj)
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       % New method to set a coordinate system can be written here
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
     end
 
   end
 
 
 
   methods (Access = private, Hidden = true)
 
     function setOriginWithScapulaLandmarks(obj,scapula)
       % Set the origin of the coordinate system to the spino-gemoid notch
       % projected on the scapular axis
       spinoGlenoid         = scapula.spinoGlenoidNotch;
       grooveMeanProjection = mean(scapula.plane.projectOnPlane(scapula.groove));
       obj.origin = grooveMeanProjection + ...
                    dot((spinoGlenoid - grooveMeanProjection),obj.zAxis)*obj.zAxis;
     end
 
     function setXAxisWithScapulaLandmarks(obj,scapula)
       % Set the postero-anterior axis to be the scapula plane normal
       obj.xAxis = scapula.plane.normal;
     end
 
     function setYAxisWithScapulaLandmarks(obj,scapula)
       % Set the infero-superior axis to be orthogonal with the two other axes
       superior    = scapula.spinoGlenoidNotch;
       inferior    = scapula.angulusInferior;
       obj.yAxis   = cross(obj.xAxis,obj.zAxis);
       obj.yAxis   = orientVectorToward(obj.yAxis,(superior-inferior));
     end
 
     function setZAxisWithScapulaLandmarks(obj,scapula)
       % Set the media-lateral axis to be the line fitted to the projection of
       % the groove points on the scapula plane
       lateral     = scapula.spinoGlenoidNotch;
       medial      = scapula.trigonumSpinae;
       groovePointsProjection = scapula.plane.projectOnPlane(scapula.groove);
       grooveAxis  = fitLine(groovePointsProjection);
       grooveAxis  = (grooveAxis/norm(grooveAxis))';
       obj.zAxis   = orientVectorToward(grooveAxis,(lateral-medial));
     end
 
   end
 
 
 
 end