%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))])
