filename = "/Users/andy2/Documents/work/Rear Right.stl"
ps = 72/25.4 % conversion from mm to points
thickness = 0.302 % thickness of layers
[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

xmin= min(X(:))
xmax = max(X(:))
ymin= min(Y(:))
ymax = max(Y(:))
zmin= min(Z(:))
zmax = max(Z(:))

%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

clg
plot([0;0;xmax-xmin;xmax-xmin;0] , [ymax-ymin;0;0;ymax-ymin;ymax-ymin])
hold on
legend("off")

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))-xmin;
	x2 = xs(2)+(xs(3)-xs(2))*(z-zs(2))/(zs(3)-zs(2))-xmin;

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

	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))-xmin;
	x2 = xs(3)+(xs(1)-xs(3))*(z-zs(3))/(zs(1)-zs(3))-xmin;

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

	C = [C; x1 y1 x2 y2];


end

plot([C(:,1) C(:,3)]' , [C(:,2) C(:,4)]')
hold off

C = C * ps; % Convert to points

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

while (any(g))


		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);
		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);
				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