filename = "/Users/andy2/Documents/work/Rear Right.stl"
ps = 1 % conversion from mm RD-GL coords
thickness = 0.302 % thickness of layers
pagewidth = 285; % paper size
pageheight = 195;

[F V C] = stl_read(filename);

X = [V(F(:,1),1) V(F(:,2),1) V(F(:,3),1)]; % X coords for each face
Y = [V(F(:,1),2) V(F(:,2),2) V(F(:,3),2)]; %Y
Z = [V(F(:,1),3) V(F(:,2),3) V(F(:,3),3)]; %Z

T = Z; %quick and dirty rotation
Z = X;
X = T;

clear T

% find extents

zmin= min(Z(:))
zmax = max(Z(:))

page_x = 0; page_y = 0;

%for z = zmin:.005:zmax

z = mean([zmax zmin])% test porpoises

I1 = find(sum((Z>z),2) == 1); % Index of faces with 1 vertex greater than z
I2 = find(sum((Z>z),2) == 2); % index of faces with 2 > z

C = [];

for a = 1:length(I1)
	[zs i] = sort(Z(I1(a),:)) ;% put biggest Z last, and store order
	xs = X(I1(a),i);
	ys = Y(I1(a),i);

	x1 = xs(1)+(xs(3)-xs(1))*(z-zs(1))/(zs(3)-zs(1));
	x2 = xs(2)+(xs(3)-xs(2))*(z-zs(2))/(zs(3)-zs(2));

	y1 = ys(1)+(ys(3)-ys(1))*(z-zs(1))/(zs(3)-zs(1));
	y2 = ys(2)+(ys(3)-ys(2))*(z-zs(2))/(zs(3)-zs(2));

	C = [C; x1 y1 x2 y2];

end

for a = 1:length(I2)

	[zs i] = sort(Z(I2(a),:)) ;% put biggest Z last, and store order
	xs = X(I2(a),i);
	ys = Y(I2(a),i);

	x1 = xs(2)+(xs(1)-xs(2))*(z-zs(2))/(zs(1)-zs(2));
	x2 = xs(3)+(xs(1)-xs(3))*(z-zs(3))/(zs(1)-zs(3));

	y1 = ys(2)+(ys(1)-ys(2))*(z-zs(2))/(zs(1)-zs(2));
	y2 = ys(3)+(ys(1)-ys(3))*(z-zs(3))/(zs(1)-zs(3));

	C = [C; x1 y1 x2 y2];


end

xmax = max(max(C(:,[1 3])));
xmin = min(min(C(:,[1 3])));
ymax = max(max(C(:,[2 4])));
ymin = min(min(C(:,[2 4])));

clg
hold on
plot (xmin, ymin, 'b');
legend("off")

%plot([0;0;xmax-xmin;xmax-xmin;0] , [ymax-ymin;0;0;ymax-ymin;ymax-ymin])
%plot([C(:,1) C(:,3)]' , [C(:,2) C(:,4)]')
%hold off

if (xmax-xmin > pagewidth | ymax-ymin > pageheight)
	disp("Page not large enough for object")
	stop
endif

C = C * ps; % Convert to points/mm/whatever

if (page_y + ymax - ymin > pageheight) % top of page reached

	if (f ~= 0) % file already open

		fprintf(f, '%%%%EOF\n')
		fclose(f)
		fputs(p, "SP 0;\n")

		% do something with plotter here
	endif
endif


f = fopen(['/Users/andy2/Documents/gash' num2str(z) '.eps'], 'wt');

fprintf(f, '%%!PS-Adobe-3.0 EPSF-3.0\n')
fprintf(f, '%%%%BoundingBox: %0.2f %0.2f %0.2f %0.2f\n', 0, 0, (xmax-xmin)*ps*1.02, (ymax-ymin)*ps*1.02)
fprintf(f, '0.75 setlinewidth\n')

%Try to string line segments together into loops.

g = 1:length(C); % use this vector to store indices to the currently unassigned segments

lastx = 0 ; lasty = 0;

while (any(g))
	% find nearest (x,y) to (lastx, lasty)

		gi = 1; % starting position

		fprintf(f, '%0.2f %0.2f moveto\n', C(g(gi), 1)-xmin*ps, C(g(gi), 2)-ymin*ps);
		fprintf(f, '%0.2f %0.2f lineto\n', C(g(gi), 3)-xmin*ps, C(g(gi), 4)-ymin*ps);
		plot([C(g(gi), 1) C(g(gi), 3)] , [C(g(gi), 2) C(g(gi), 4)], 'b');

		lastx = C(g(gi),3) ; lasty = C(g(gi), 4);

	do

		g = [g(1:gi-1) g(gi+1:end)]; % remove used index from list

		 % find x and y match looking in both start and end
		gi = [find((C(g,1) == lastx) & (C(g,2) == lasty)) ; ...
		     find((C(g,3) == lastx) & (C(g,4) == lasty))];

		if any(gi)
			gi = gi(1); % choose the first match
			xi = [3 1](find(C(g(gi),[1 3]) == lastx)); % find the new positions
			yi = [4 2](find(C(g(gi),[2 4]) == lasty));
			if (length(xi) == 1 | length(yi) == 1) % needed for zero-length segments
				fprintf(f, '%0.2f %0.2f lineto\n', ...
				 C(g(gi), xi(1))-xmin*ps, C(g(gi), yi(1))-ymin*ps);

				plot([lastx C(g(gi), xi(1))], [lasty C(g(gi), yi(1))], 'b')

				lastx = C(g(gi),xi(1)) ; lasty = C(g(gi), yi(1));
			endif
		endif

		if !any(gi)
			fprintf(f, 'stroke\n')
		endif

	until !any(gi) % ie no matches for last end point


endwhile



fprintf(f, '%%%%EOF\n')

fclose(f)

%end