diff --git a/ShoulderCase/@Humerus/measureCenterAndRadius.m b/ShoulderCase/@Humerus/measureCenterAndRadius.m
new file mode 100644
index 0000000..2b4ca53
--- /dev/null
+++ b/ShoulderCase/@Humerus/measureCenterAndRadius.m
@@ -0,0 +1,6 @@
+function measureCenterAndRadius(obj)
+  % By fitting a sphere on humeral head landmarks
+  [center,radius,~,~] = fitSphere(obj.landmarks);
+  obj.center    = center';
+  obj.radius    = radius;
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Humerus/measureSubluxation.m b/ShoulderCase/@Humerus/measureSubluxation.m
new file mode 100644
index 0000000..81f1a4b
--- /dev/null
+++ b/ShoulderCase/@Humerus/measureSubluxation.m
@@ -0,0 +1,83 @@
+function measureSubluxation(obj)
+  % Scapulo-humeral subluxation (SHS)
+
+  % Vector from glenoid center to humeral head center
+  % glenoidToHH = obj.center - obj.shoulder.scapula.glenoid.center;
+  scapula = obj.shoulder.scapula;
+  glenoidToHH = Vector(scapula.glenoid.center, obj.center);
+
+  xAxis = Vector(...
+    scapula.coordSys.origin,...
+    scapula.coordSys.origin + scapula.coordSys.PA);
+  yAxis = Vector(...
+    scapula.coordSys.origin,...
+    scapula.coordSys.origin + scapula.coordSys.IS);
+  zAxis = Vector(...
+    scapula.coordSys.origin,...
+    scapula.coordSys.origin + scapula.coordSys.ML);
+
+  % Vector from humeral head center to scapular (z) axis,
+  % perpendicular to z-axis --> similar to projection in xy plane
+  glenoidToHHxy = zAxis.orthogonalComplementTo(glenoidToHH);
+
+
+
+
+  % SHS amplitude (ratio rather than amplitude)
+  obj.SHSAmpl = norm(glenoidToHHxy) / obj.radius;
+
+
+
+
+  % SHS orientation
+  obj.SHSOrient = rad2deg(angle(-xAxis, glenoidToHHxy));
+  % Superior is positive
+  obj.SHSOrient = sign(dot(yAxis, glenoidToHHxy)) * obj.SHSOrient;
+
+
+
+
+  % SHSAngle: angle between glenHumHead and zAxis
+  obj.SHSAngle = rad2deg(angle(glenoidToHH, zAxis));
+
+
+
+  % SHSPA is angle between zxProjection and zAxis
+  glenoidToHHzx = yAxis.orthogonalComplementTo(glenoidToHH);
+  obj.SHSPA = rad2deg(angle(glenoidToHHzx, zAxis));
+  % Convention that posterior is negative
+  obj.SHSPA = sign(dot(xAxis, glenoidToHHzx)) * obj.SHSPA;
+
+
+
+
+  % SHSIS is angle between yzProjection and zAxis
+  glenoidToHHyz = xAxis.orthogonalComplementTo(glenoidToHH);
+  obj.SHSIS = rad2deg(angle(glenoidToHHyz, zAxis));
+  % Convention that superior is negative
+  obj.SHSIS = sign(dot(yAxis, glenoidToHHyz)) * obj.SHSIS;
+
+
+
+
+  %% Gleno-humeral subluxation (GHS)
+  glenoidCenterLine = Vector(...
+    scapula.glenoid.center,...
+    scapula.glenoid.center + scapula.glenoid.centerLine);
+  % Vector from glenoid sphere centre to glenoid centerline (perpendicular)
+  GHS = glenoidCenterLine.orthogonalComplementTo(obj.center);
+  
+  % glenCenterlineNorm = glenCenterline/norm(glenCenterline);
+  % GHS = glenoidToHH - dot(glenoidToHH, glenCenterlineNorm)*glenCenterlineNorm;
+  % GHS  amplitude
+  obj.GHSAmpl = norm(GHS) / obj.radius;
+  % GHSAmpl = norm(GHS) / humHeadRadius;
+
+
+
+  % GHS orientation
+  GHSxy = zAxis.orthogonalComplementTo(GHS);
+  obj.GHSOrient = rad2deg(angle(-xAxis, GHSxy));
+  % Superior is positive
+  obj.GHSOrient = sign(dot(yAxis, GHSxy)) * obj.GHSOrient;
+end
\ No newline at end of file