% Matlab script wullf4 to generate Wulff nets
% I need to eliminate great circle lines within the 10-degree small circles
% I need to get the size right to match my keys

% Definition of variables
% x,y: center of arc
% r: radius of arc
% thetaa: lower limt of arc range
% thetab: upper limit of arc range

% Clear screen
clf;
figure(1)
clf
%	Set radius of Wulff net primitive circle
bigr = 1.2;
phid = [2:2:88];	% Angular range for 
phir = phid*pi/180;
omegad = 90 - phid; 
omegar = pi/2-phir;

% Set up for plotting great circles with centers along
% positive x-axis
x1 = bigr.*tan(phir);
y1 = zeros(size(x1));
r1 =  bigr./cos(phir);
theta1ad = (180-80)*ones(size(x1));
theta1ar = theta1ad*pi/180;
theta1bd = (180+80)*ones(size(x1));
theta1br = theta1bd*pi/180;

% Set up for plotting great circles 
% with centers along the negative x-axis
x2 = -1*x1;
y2 = y1;
r2 = r1;
theta2ad = -80*ones(size(x2));
theta2ar = theta2ad*pi/180;
theta2bd =  80*ones(size(x2));
theta2br = theta2bd*pi/180;


% Set up for plotting small circles
% with centers along the positive y-axis
y3 = bigr./sin(omegar);
x3 = zeros(size(y3));
r3 =  bigr./tan(omegar);
theta3ad = 3*90-omegad;
theta3ar = 3*pi/2-omegar;
theta3bd = 3*90+omegad;
theta3br = 3*pi/2+omegar;

% Set up for plotting small circles
% with centers along the negative y-axis
y4 = -1*y3;
x4 = x3;
r4 =  r3;
theta4ad = 90-omegad;
theta4ar = pi/2-omegar;
theta4bd = 90+omegad;
theta4br = pi/2+omegar;

% Group all x, y, r, and theta information for great cricles 
phi = [phid, phid];
x = [x1, x2];
y = [y1, y2];
r = [r1, r2];
thetaad = [theta1ad, theta2ad];
thetaar = [theta1ar, theta2ar];
thetabd = [theta1bd, theta2bd];
thetabr = [theta1br, theta2br];
	
% Plot portions of all great circles that lie inside the
% primitive circle, with thick lines (1 pt.) at 10 degree increments
for i=1:length(x)
 	thd = thetaad(i):1:thetabd(i);
 	thr = thetaar(i):pi/180:thetabr(i);
	xunit = x(i) + r(i).*cos(thr);
	yunit = y(i) + r(i).*sin(thr);
	p = plot(xunit,yunit,'LineWidth',0.5);
	hold on
end

% Now "blank out" the portions of the great circle cyclographic traces 
% within 10 degrees of the poles of the primitive circle.
rr =  bigr./tan(80*pi/180);
ang1 = 0:pi/180:pi;
xx = zeros(size(ang1)) + rr.*cos(ang1);
yy = bigr./cos(10*pi/180).*ones(size(ang1)) - rr.*sin(ang1);
p = fill(xx,yy,'w')
yy = -bigr./cos(10*pi/180).*ones(size(ang1)) + rr.*sin(ang1);
p = fill(xx,yy,'w')

for i=1:length(x)
 	thd = thetaad(i):1:thetabd(i);
 	thr = thetaar(i):pi/180:thetabr(i);
	xunit = x(i) + r(i).*cos(thr);
	yunit = y(i) + r(i).*sin(thr);
	if mod(phi(i),10) == 0
		p = plot(xunit,yunit,'LineWidth',1);
		angg = thetaad(i)
	end
	hold on
end

% Now "blank out" the portions of the great circle cyclographic traces 
% within 2 degrees of the poles of the primitive circle.
rr =  bigr./tan(88*pi/180);
ang1 = 0:pi/180:pi;
xx = zeros(size(ang1)) + rr.*cos(ang1);
yy = bigr./cos(2*pi/180).*ones(size(ang1)) - rr.*sin(ang1);
p = fill(xx,yy,'w')
yy = -bigr./cos(2*pi/180).*ones(size(ang1)) + rr.*sin(ang1);
p = fill(xx,yy,'w')


% Group all x, y, r, and theta information for small circles
phi = [phid, phid];
x = [x3, x4];
y = [y3, y4];
r = [r3, r4];
thetaad = [theta3ad, theta4ad];
thetaar = [theta3ar, theta4ar];
thetabd = [theta3bd, theta4bd];
thetabr = [theta3br, theta4br];

% Plot primitive circle
thd = 0:1:360;
thr = 0:pi/180:2*pi;
xunit = bigr.*cos(thr);
yunit = bigr.*sin(thr);
p = plot(xunit,yunit);
hold on
	
% Plot portions of all small circles that lie inside the
% primitive circle, with thick lines (1 pt.) at 10 degree increments
for i=1:length(x)
 	thd = thetaad(i):1:thetabd(i);
 	thr = thetaar(i):pi/180:thetabr(i);
	xunit = x(i) + r(i).*cos(thr);
	yunit = y(i) + r(i).*sin(thr);
%	blug = mod(thetaad(i),10)
	if mod(phi(i),10) == 0
		p = plot(xunit,yunit,'LineWidth',1);
		angg = thetaad(i)
	else
		p = plot(xunit,yunit,'LineWidth',0.5);
	end
	hold on
end

% Draw thick north-south and east-west diameters
xunit = [-bigr,bigr];
yunit = [0,0];
p = plot(xunit,yunit,'LineWidth',1);
hold on
xunit = [0,0];
yunit = [-bigr,bigr];
p = plot(xunit,yunit,'LineWidth',1);
hold on

% Parameters to control appearance of plot
% THESE COME AFTER THE PLOT COMMANDS!!!
axis([-bigr bigr -bigr bigr])
% axis ('square').  BAD way to get aspect ratio of plot.  It
% also considers titles and axis labels when scaling the figure!
set(gca,'DataAspectRatio',[bigr,bigr,bigr])
%axes('Position',[0,0,1,1]);
%axes('AspectRatio',[1,1]);
set(gca,'Visible','off');		% This turns off the visibility of the axes
% figure('PaperPosition',[1,3,6,6]);
print -dill wulffnet.ill
print -deps wulffnet.eps
% end