%Horn design 3
%Gordon C. Kirkwood December 2013
%We desire a horn whose cross sectional area changes, over the length of
%the horn, with an exponential function of the general form:
%area = A0 * e^(m*x)
%radius = r0*e^(m*x/2)
% x is axis of horn, position
%intercept of initial slope of horn in 2/m from A0 coordinate.
%what factor m to use? solve for M based on desired height, initial
%radius, final radius.
%then solve for radius as a function of height.
%then figure coordinates of 2 edges of a gore panel, a panel of the faceted cone.
%then figure the distance ALONG THE curved central line of each gore / panel (as opposed to previous
%parametrization, which was according to height), of each coordinate;
%then 'flatten' by recording new gore panel DXF:
close all
clear all
addpath('/Users/gordonkirkwood/Documents/MATLAB/DXFlib'); %Change to point to the location of DXF writing functions.
%DXFlib can be downloaded at: http://www.mathworks.com/matlabcentral/fileexchange/33884-dxflib
%all length units in inches
numfacets = 5 %the number of facets of the faceted exponential horn.
%
height = 43 %inches length of centra axis of horn
dia1=34 %diameter of horn at widest end, measured as the diameter of the circle containing the faceted cross section (not the inscribed circle).
dia0= 3 %diameter of the horn at narrowest end
%make the filename encode the parameters of the horn, so you don't
%accidentally overwrite or mislabel a file.
fname = ['dia_',num2str(dia1),'to',num2str(dia0),'_H',num2str(height),'_facets_',num2str(numfacets)];
filename = ['exponential_horn_',fname,'.dxf'];
%%%%%%%%%%%%%%%%%%%%%%%%
r1 = dia1/2
r0 = dia0/2
expbase = 2
m = log(r1/r0)*2/height
samples = 30;
h = linspace(0,height,samples)';
r = r0*exp(m*h/2);
figure(2)
subplot(2,1,1);
plot(h,r,h,-1*r); axis equal;
subplot(2,1,2);
plot(h,pi*r.^2)
title('horn shape, +/- radius vs. axial position')
%draw horn
theta = linspace(0,2*pi,60);
figure(1)
for ih =1:size(h,1)
plot3(r(ih).*cos(theta),r(ih).*sin(theta),ones(size(theta))*h(ih))
hold on;
end
axis equal
%now figure the coordinates of each point on perimeter, rotated around
%central axis by (2*pi/numfacets angle.
angle = [0:numfacets-1].*2*pi/numfacets %radians
uv= exp(i*angle); %unit vector in new angle direction, one for each entry in angle vector
uvs = repmat(uv,samples,1);
rs = repmat(r,1,numfacets);
coords = uvs.*rs;
%ok, coords now has (complex number xy coordinates) for each angle position
%of cone.
%now figure the length accumulated between each sample position. this
%reparametrizes curve from height to 'length along curve' necessary for
%flattening.
%
c_x= real(coords);
c_y = imag(coords);
dx= diff(c_x')'; %distance between points on the same layer of horn, on different radial subdivisions. basically, the radius times the included chord of the number of facets
dy = diff(c_y')';
dist = sqrt(dx.^2 + dy.^2); %<--- the width of a panel as a function of sample
centerline = (coords(:,1) + coords(:,2))/2;
center_x = real(centerline);
center_y = imag(centerline);
dxc = diff(center_x);
dyc = diff(center_y);
dh = diff(h);
dl = sqrt(dxc.^2 + dyc.^2 + dh.^2);
plot3(c_x,c_y,h);%coords)
hold on
plot3(center_x,center_y,h,'color','c'); %looks good.
%now I've got the panel widths (dist) as a function of height. but I need
%them as a function of panel length. since the panels are not parallel to
%the axis, and indeed curve (it is an exponential horn with exponent m!=0)
%height is not the same as curve position. measure height by measuring
%distance between each successive point, along a given XY curve.
l = cumsum([0; dl],1);
figure
plot(l,dist/2, l, -1*dist/2,'r-')
title('facet shape');
hold on
%new coordinates are 1...samples
%+dist(:)/2,l(:)
%-dist(:)/2,l(:)
Y = [dist(2:end,1)';dist(1:end-1,1)']./2;
X = [l(2:end)'; l(1:end-1)';];
Xright = [X [l(end); l(end)]]
Yright = [Y [dist(end,1)/2; -1*dist(end,1)/2]]
line(Xright,Yright,'color','g');
hold on
Xleft = [X(:,end:-1:1),[0; 0]];
Yleft = [-1*Y(:,end:-1:1), [-1*dist(1)/2; dist(1)/2]];
line(Xleft,Yleft,'color','g');
axis equal
title('line plots')
totalX = [Xleft,Xright];
totalY = [Yleft,Yright];
figure
line(totalX,totalY)
title('total X and Y line matricies drawn. ');
axis equal
%OK, looks good, now to
%now test- debug
figure(1)
FID = dxf_open(filename);
FID = dxf_set(FID,'Color',[1 1 0],'Layer',20);
Z = zeros(size(totalY));
dxf_polyline(FID,totalX,totalY,Z);
dxf_close(FID);
disp([filename,' created'])
disp(['total length =',num2str(max(l))])