function rockvis3(ad,bd,hem,cle)

% General properties.
if (nargin==2)
   str='Usage: rockvis3(ad, bd, l or u, c). Default parameter l is used.';
   warning(str);
   hem='l';
elseif (cle=='c')
   close; close; close;
else 
   cle='h';
end

%-------------------------------------------------------------------------------------------
%%                               SPHERE SCRIPT
% This script constructs the ('l')ower or ('u')pper hemisphere of a sphere of radius R.

n=100; % n is the number of points in the drawing of the sphere.
R=1;   % OLD VERSION: R=100;
theta=pi*(0:5:2*n)/n;
if hem=='l', % phi defines the hemisphere (L,U).
   phi=(pi/2)*(-n:5:0)'/n; 
   h=-1;
elseif hem=='u', 
   phi=(pi/2)*(0:5:n)'/n;
   h=1;
end 
r=R*cos(phi);
% Parametric Equation of a sphere X(r,theta,phi),Y(r,theta,phi),Z(r,theta,phi):
xsphere=r*cos(theta);
ysphere=r*sin(theta);
zsphere=R*sin(phi)*ones(size(theta));
%% END SPHERE SCRIPT

%-------------------------------------------------------------------------------------------
%%                               PLANE SCRIPT
% an/bn are the dip direction/dip amount of the normal to the inclined plane n.
an=ad+180;	if an>360 an=an-360; end
bn=-h*(90-bd); 

% Reverse calculation:
%   ad=an-2*pi; 		
%   bd=pi/2-bn;

% Convert degrees to radians:
   ad=ad*pi/180;
   bd=bd*pi/180;
   an=an*pi/180;
   bn=bn*pi/180;
   
   %--------
   m=10; % m is the number of points into which the Range (-R,R) is  divided.
   % The RANGE of the function Z=Z(X,Y) on the xy plane is:
   [xplane,yplane]=meshgrid(-R:R,-R:R); 
   
   % Direction cosines of a vector n in xyz-space:
   ax=sin(an)*cos(bn);
   ay=cos(an)*cos(bn);
   az=h*sin(bn);
   
   % Definition of the PLANE that passes through O(0,0,0) and is perpendicular to vector n
   % in parametric form:
   zplane=-ax/az*xplane-ay/az*yplane;
%% END PLANE SCRIPT

%===========================================================================================
% The CIRCLE OF PROJECTION (PROJECTION PLANE sigma) [of radius R]
theta=pi*(0:2*n)/n;

xproj=R*cos(theta);
yproj=R*sin(theta);
%=======================================

% The normal to the plane
if hem=='l', aN=an; bN=bn; else aN=an+pi; bN=-bn; end
xN = R*sin(aN)*tan(pi/4-bN/2);
yN = R*cos(aN)*tan(pi/4-bN/2);


% The Lowest Point of the Line of Maximum Dip (ML)
xM=-R*sin(an)*tan(bn/2);
yM=-R*cos(an)*tan(bn/2);

%-----
% The physical center of the great circle (GC or g)
rg=R*tan(pi/2-bn);
xC=rg*sin(an);
yC=rg*cos(an);

Rg=R/cos(bd); % the radius of the GC. Check: Rg-sqrt(R^2+rg^2)==0 

% A,B is the projections of the outermost points of the GC
xA=R*cos(an);
yA=-R*sin(an);

xB=-R*cos(an);
yB=R*sin(an);    %% Check: sqrt(xB^2+yB^2)==R
% Translate to x1y1 reference system:
x1B=xB-xC;       
y1B=yB-yC;

a1B=atan(y1B/x1B); 
if x1B<0 
   a1B=a1B+pi; 
elseif y1B<0 
   a1B=a1B+2*pi;
end

azGC=linspace(a1B,a1B+2*bn,n);
x1GCproj=Rg*cos(azGC); 
y1GCproj=Rg*sin(azGC);

%-----------------------------------------------------------------------
% THE PROJECTION OF THE GREAT CIRCLE GC ON THE PROJECTION PLANE sigma=0.
%
xGCproj = x1GCproj + xC*ones(size(x1GCproj));
yGCproj = y1GCproj + yC*ones(size(y1GCproj));

%===============================================================================================
%                                INVERSE TRANSFORMATION
% The 3-DIMENSIONAL INTERSECTION of the SPHERE and the PLANE is:
%if hem == 'l', 
   k = 2./(xGCproj.^2+yGCproj.^2+1); % Analogy factor.
   Z = (ones(size(k))-k);            % Z=1-k=(k-1)*(-1) ; z=-1 for the sigma plane.
%elseif hem == 'u',
%   k=-4./(xGCproj.^2+yGCproj.^2+4);  % Z=2*k+1=(2*k+1)*1 ; z=1 for the sigma plane.
%   Z = 2*k+ones(size(k));
%end   
X = k.*xGCproj;
Y = k.*yGCproj;

%plot3(X,Y,Z);

%===============================================================================================
% Compute the angles:
% Y./X
% Z./R
a=atan(Y./X) + (1-sign(X))/2*pi + (1-sign(Y))/2*2*pi;

b=asin(-Z./R);

%---------------------------------------------------
% Project:
r =1.4142*R*cos(pi/4+b/2); % 1.4142 = sqrt(2)

xEA=r.*sin(a);
yEA=r.*cos(a);

%=============================================================================================
%                                       PLOTTING
%=============================================================================================
% function plotHemProj
% close; close; close; close;

% --------------------------------------------------------------------------------------------
%                                    PLOT THE 3D PART
figure(1);
mesh(xsphere, ysphere, zsphere, ones(size(zsphere)));
axis equal;

hold on;    % Holds previous screen so that 
%--
surf(xplane, yplane, zplane, ones(size(zplane)));

%% view(-20,50);
title('3D VIEW OF THE HEMISPHERICAL PROJECTION');
% colormap hot;
%-- OCTAVE/GNUPLOT code ---------
%   gset parametric;
%   gsplot(xplane,yplane,zplane); 
%--------------------------------
rotate3d on;
colormap lines;
% ----------------------------------------------------------------------------------------
%                                     PLOT THE 2D PART
%-----------------------------------------------------------------------------------------

%---------------------------------- EQUAL-ANGLE PROJECTION -------------------------------
figure(2);
hold on; 
plot(0,0,'+');               
plot(xproj,yproj);          
plot(xM,yM,'o');
text(xM,yM,' M');
plot(xN,yN,'*');
text(xN,yN,' N'); 

if bd<=1.05 % 1.05 rad ~= 60 deg.
   plot(xC,yC,'x');
   text(xC,yC,' C');
end
plot(xB,yB,'o');
text(xB,yB,' B');  
plot(xA,yA,'o');
text(xA,yA,' A'); 
plot(xGCproj,yGCproj,'m');    % The great circle
axis equal;

%str1=strcat('EQUAL-ANGLE PROJECTION ON A PLANE: ad/bd=',num2str(ad*180/pi),'/',num2str(bd*180/pi));
%h1=text(-R,R+0.145,str1);
title(['EQUAL-ANGLE PROJECTION ON A PLANE: ad/bd=',num2str(ad*180/pi),'/',num2str(bd*180/pi)]);
%set(h1,'EraseMode','normal');
text(0,R+0.025,'|N(0)');
text(R,0,'--E(90)');
text(0-0.005,-R-0.025,'| S(180)');
text(-R-0.28,0,'W(270)--');
%-------------------------------------------------------------------------------------
%                            PLOT THE EQUAL AREA PROJECTION
figure(1);
rotate3d on;
%-----------
figure(3);
[XEA YEA]=pol2cart(a,r);
plot(XEA,YEA,'r');
hold on;
plot(0,0,'+');               
plot(xproj,yproj); 
axis equal;

%--- Test: Plot (XA,YA,ZA), (XB,YB,ZB), (XM, YM, ZM)
aA=ad-pi/2;
rA=1.4142*R*cos(pi/4+0/2);
plot(rA*sin(aA), rA*cos(aA), 'ro');
text(rA*sin(aA), rA*cos(aA), ' A');

aB=ad+pi/2;
rB=1.4142*R*cos(pi/4+0/2);
plot(rB*sin(aB), rB*cos(aB), 'ro');
text(rB*sin(aB), rB*cos(aB), ' B');

%aM==ad.
if hem=='l', aM=ad; else aM=ad+pi; end
rM=1.4142*R*cos(pi/4+bd/2);
plot(rM*sin(aM), rM*cos(aM), 'ro');
text(rM*sin(aM), rM*cos(aM), ' M');

%text(-R,R+0.145,'');
%str1=strcat('EQUAL-AREA PROJECTION ON A PLANE: ad/bd=',num2str(ad*180/pi),'/',num2str(bd*180/pi));
%h2=text(-R,R+0.145,str1);
title(['EQUAL-AREA PROJECTION ON A PLANE: ad/bd=',num2str(ad*180/pi),'/',num2str(bd*180/pi)]);
%set(h2,'EraseMode','normal');
text(0,R+0.01,'|N(0)');
text(R,0,'--E(90)');
text(0-0.005,-R-0.025,'| S(180)');
text(-R-0.28,0,'W(270)--');
%==================================================================================================
