diff --git a/ShoulderCase/@Shoulder/Shoulder.m b/ShoulderCase/@Shoulder/Shoulder.m
index 4852a0f..21a14b9 100644
--- a/ShoulderCase/@Shoulder/Shoulder.m
+++ b/ShoulderCase/@Shoulder/Shoulder.m
@@ -1,39 +1,230 @@
 classdef (Abstract) Shoulder < handle
 
   properties
     side % R or L
     scapula
     humerus
     muscles
     CTscan
     comment
   end
 
 
 
   properties (Hidden = true)
     SCase
   end
 
 
 
   methods (Abstract, Access = protected)
     createScapula(obj);
   end
 
 
 
   methods
     function obj = Shoulder(SCase)
       obj.side    = '';
       obj.humerus = Humerus(obj);
       obj.muscles = Muscles(obj);
       obj.CTscan  = '';
       obj.comment = '';
 
       obj.SCase   = SCase;
 
       obj.createScapula;
     end
+
+    function output = plot(obj)
+      % Create a figure and plot inside:
+      % the scapula surface, coordSys, groove points and landmarks;
+      % the glenoid surface, sphere and centerLine;
+      % the humerus sphere and centerLine.
+
+      figure;
+      hold on;
+
+
+      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      %% Scapula
+      % Plot the segmented surface
+      if not(isempty(obj.scapula.surface))
+          points = obj.scapula.surface.points;
+          faces  = obj.scapula.surface.faces;
+          x = points(:,1);
+          y = points(:,2);
+          z = points(:,3);
+          trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',0.5,'EdgeColor','none','FaceLighting','gouraud');
+      end
+
+      % Plot scapula groove
+      if not(isempty(obj.scapula.groove))
+        points = obj.scapula.groove;
+        x = points(:,1);
+        y = points(:,2);
+        z = points(:,3);
+        scatter3(x,y,z,100,'blue');
+      end
+
+      % Plot scapula markers
+      points = [...
+          obj.scapula.angulusInferior;
+          obj.scapula.trigonumSpinae;
+          obj.scapula.processusCoracoideus;
+          obj.scapula.acromioClavicular;
+          obj.scapula.angulusAcromialis;
+          obj.scapula.spinoGlenoidNotch...
+          ];
+      x = points(:,1);
+      y = points(:,2);
+      z = points(:,3);
+      scatter3(x,y,z,100,'blue','LineWidth',2);
+
+
+      % Plot scapular plane
+      axisLength = 50;
+      pt0 = obj.scapula.coordSys.origin;
+      pt1 = pt0 + (obj.scapula.coordSys.ML + obj.scapula.coordSys.IS) * axisLength ;
+      pt2 = pt1 - obj.scapula.coordSys.IS * pdist([pt1 ; obj.scapula.angulusInferior]);
+      pt3 = pt2 - obj.scapula.coordSys.ML * pdist([pt1 ; obj.scapula.trigonumSpinae] );
+      pt4 = pt3 + obj.scapula.coordSys.IS * pdist([pt1 ; obj.scapula.angulusInferior]);
+      X = [pt1(1) pt2(1) pt3(1) pt4(1)];
+      Y = [pt1(2) pt2(2) pt3(2) pt4(2)];
+      Z = [pt1(3) pt2(3) pt3(3) pt4(3)];
+      patch(X,Y,Z,'black','FaceAlpha',0.3); % Transparent plane
+
+      % Plot scapular axis
+      pt1 = obj.scapula.coordSys.origin;
+      axisLength = 10 + pdist([pt1 ; obj.scapula.trigonumSpinae]);
+      pt2 = pt1 - obj.scapula.coordSys.ML * axisLength;
+      x = [pt1(1), pt2(1)];
+      y = [pt1(2), pt2(2)];
+      z = [pt1(3), pt2(3)];
+      plot3(x,y,z,'Color','blue','LineWidth', 5);
+
+      % Plot scapular coordinate system
+      axisLength = 50; % Length of the coord. sys. axes
+      pt0 = obj.scapula.coordSys.origin;
+      pt1 = pt0 + obj.scapula.coordSys.PA*axisLength; % X
+      pt2 = pt0 + obj.scapula.coordSys.IS*axisLength; % Y
+      pt3 = pt0 + obj.scapula.coordSys.ML*axisLength; % Z
+      x = [pt0(1), pt1(1)];
+      y = [pt0(2), pt1(2)];
+      z = [pt0(3), pt1(3)];
+      plot3(x,y,z,'Color','red','LineWidth', 5);
+      x = [pt0(1), pt2(1)];
+      y = [pt0(2), pt2(2)];
+      z = [pt0(3), pt2(3)];
+      plot3(x,y,z,'Color','green','LineWidth', 5);
+      x = [pt0(1), pt3(1)];
+      y = [pt0(2), pt3(2)];
+      z = [pt0(3), pt3(3)];
+      plot3(x,y,z,'Color','blue','LineWidth', 5);
+
+      % plot scapular axis on glenoid center
+      axisLength = 50;
+      pt1 = obj.scapula.glenoid.center;
+      pt2 = pt1 + obj.scapula.coordSys.ML*axisLength;
+      x = [pt1(1), pt2(1)];
+      y = [pt1(2), pt2(2)];
+      z = [pt1(3), pt2(3)];
+      plot3(x,y,z,'Color','magenta','LineWidth', 5);
+
+
+      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      % Glenoid
+
+      % Plot the glenoid surface
+      % We might move them in the direction of th glenoid center
+      % instead of the ML axis
+
+      points = obj.scapula.glenoid.surface.points +...
+          obj.scapula.coordSys.ML * 0.2;
+      faces  = obj.scapula.glenoid.surface.faces;
+      x = points(:,1);
+      y = points(:,2);
+      z = points(:,3);
+      trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','magenta','FaceLighting','none');
+
+
+      % plot glenoid centerline
+      pt1 = obj.scapula.glenoid.center;
+      pt2 = pt1 + obj.scapula.glenoid.centerLine;
+      x = [pt1(1), pt2(1)];
+      y = [pt1(2), pt2(2)];
+      z = [pt1(3), pt2(3)];
+      plot3(x,y,z,'Color','magenta','LineWidth', 5);
+
+      % We might add spherical cap
+
+      % plot glenoid sphere
+      n = 20; % number of lines
+      [X,Y,Z] = sphere(n);
+      X = X*obj.scapula.glenoid.radius;
+      Y = Y*obj.scapula.glenoid.radius;
+      Z = Z*obj.scapula.glenoid.radius;
+      Xc = pt2(1);
+      Yc = pt2(2);
+      Zc = pt2(3);
+      mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
+
+
+      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      % Humerus
+
+      % plot humeral head sphere
+      n = 20; % number of lines
+      [X,Y,Z] = sphere(n);
+      X = X*obj.humerus.radius;
+      Y = Y*obj.humerus.radius;
+      Z = Z*obj.humerus.radius;
+      Xc = obj.humerus.center(1);
+      Yc = obj.humerus.center(2);
+      Zc = obj.humerus.center(3);
+      mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
+
+
+      % plot humeral head centerline
+      pt1 = obj.scapula.glenoid.center;
+      pt2 = obj.humerus.center;
+      x = [pt1(1), pt2(1)];
+      y = [pt1(2), pt2(2)];
+      z = [pt1(3), pt2(3)];
+      plot3(x,y,z,'Color','magenta','LineWidth', 5);
+
+
+      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      % Graph properties
+
+      axis off; % remove axis
+      lightPosition = obj.scapula.glenoid.center + ...
+          (...
+          obj.scapula.coordSys.ML + ...
+          obj.scapula.coordSys.IS + ...
+          obj.scapula.coordSys.PA ...
+          ) * 100;
+      light('Position',lightPosition,'Style','local');
+      grid on;
+      xlabel('X (CT)');
+      ylabel('Y (CT)');
+      zlabel('Z (CT)');
+      axis square;
+      axis vis3d; % fixed limits and scaling
+      viewPoint = obj.scapula.coordSys.ML + ...
+                  obj.scapula.coordSys.IS + ...
+                  obj.scapula.coordSys.PA;
+      viewPoint = obj.scapula.coordSys.PA;
+      view(viewPoint); % view from lateral side
+      % Align IS axis with window vertical axis
+      camup(obj.scapula.coordSys.IS);
+
+      figureHandle = gca;
+      figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis
+      figureHandle.Box = 'On';
+      set(gcf,'color','w'); % Set background color to white
+
+      output = 1;
+    end
   end
 end