function h=beachball(mag, strike, dip);
%BEACHBALL--produces beachball plot of a simple earthquake
%
% h=beachball(mag,strike,dip)
%
% mag=magnitude of the earthquake (radius of ball).  This should be
%     a normalized value (e.g. between 0 and 1)
%
% strike=strike angle in degrees.  This is the direction (compass
%     heading) of the fault along the surface of the earth.  The 
%     values should be between -90 and 90 degrees (-90=west, 0=north
%     90=east--values outside this range are symmetric)
%
% dip=dip angle in degrees.  This is the angle at which the fault
%     plane intersects the earth's surface.  The values are between
%     0 and 180
%
% Returns two handles: h(1)=handle to circle, h(2)=handle to sector

strikerange=[-90 90];
diprange=[0 180];
magrange=[0 1];
if(strike<strikerange(1) | strike> strikerange(2))
	error(['strike angle should be between -90 and 90 degrees']);
end
if(dip<diprange(1) | dip> diprange(2))
	error(['dip angle should be between 0 and 180 degrees']);
end
if(mag<magrange(1) | mag> magrange(2))
	error(['mag should be between 0 and 1 (normalized units)']);
end

deg2rad=pi/180;

strike=deg2rad*strike;
cols=['w','k'];%colors for ball, cols(1)=background, cols(2)=between planes
if(dip>90);
	cols=cols([2,1]);  %we need to flip colors
   dip=dip-90;
end
dip=deg2rad*dip;

m=100;
z=linspace(-1,1,m)';

%1) Compute equation for first plane using dip
%   The equation for the intersection of this plane
%   with the surface is r=1,longitude=dip
%
%	 In cartesian coordinates, x=cos(dip)*cos(lat)
%                             z=sin(lat)
%
%	 Solving for lat in terms of z and subsituting into x equation
%   yields:
%                             x=cos(dip)*sqrt(1-z^2)
%
coslat=sqrt(1-z.*z);%cosine of latitude
x1=zeros(m,2);
x1(:,1)=cos(dip)*coslat;
x1(:,2)=z;

%2) Compute equation for second plane which is at lon=dip+90
%
x2=zeros(m,2);
x2(:,1)=cos(dip+0.5*pi)*coslat;
x2(:,2)=z;


%Now, we need to rotate both planes by strike so that they intersect
%at lon=0 and lat=strike
%
cs=cos(strike);ss=sin(strike);
rotmat=[cs,-ss;ss,cs];%rotation matrix

x1=x1*rotmat;
x2=x2*rotmat;

%stack x1 and x2 so we can patch
x1=[x1;flipud(x2(2:m-1,:))];

%Produce the graphics
holdstate=get(gca,'nextplot');%save hold status
theta=linspace(0,2*pi,200);
h(1)=patch(mag*cos(theta),mag*sin(theta),cols(1));hold on

h(2)=patch(mag*x1(:,1),mag*x1(:,2),cols(2));
set(gca,'nextplot',holdstate);
