%
% Tveggja NACAmptt prófíla-reiknir sem táknbreytir ferlum yfir
%
% í G-kóða fyrir frauðplastskera sem staðsettur er í vindgangarými
%
% Háskólans í Reykjavík.
%
%
%
% Innlagsbreytur: m, p, t, c, b, TEo, TEv, inc, n, filename
%
% Hver innlagsbreyta fyrir sig er útskýrð þar sem hún kemur fyrir í
%
% kóðanum
%
%
%
% Úttaksbreyta: G-kóði sem skrifaður er í skránna 'filename.gcode'
%
%
%
% Höfundur: Kristján Orri Magnússon, Háskólinn í Reykjavík, Haust 2011
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear all
% Nafn úttaksskráar. Endingin skal vera .gcode. 
filename=['NACA0015test.gcode']
%One digit max camber as % of chord, m2 er prófíllinn sem að A-B plan
%frauðplastskerans teiknar og m er fyrir X-Y planið.
m=0 
m2=0
%One digit describing the distance of maximum camber from the airfoil
%leading edge in tens of percents of the chord. p->X-Y plan, p2->A-B plan 
p=0
p2=0
%Two digits describing maximum thickness of the airfoil as percent of the
%chord. t-> X-Y plan, t2-> A-B plan 
t=15
t2=15
%chord length (mm). c ->X-Y plan c2->A-B plan 
c=130
c2=130
% Span length(mm). Vænghaf vængs sem skera á út.
b=500
%T/E horizontal offset (mm). Lárétt hliðrun trailing edge milli prófíla í
%XY og AB vængsniðum.
TEo=10
%T/E vertical offset (mm). Lóðrétt hliðrun trailing edge milli prófíla í
%XY og AB vængsniðum.
TEv=30
%Incline(rad). Vindingur, horn milli vænglína sitt hvors prófíls.
inc=20*pi/180
%number of points (80 is minimum for L/E precision)
n=80

%Total span (mm). Hornrétt vegalengd milli XY og AB plana, fasti. 
tsp=935
%Scale up for shorter than tsp wings 
inc=inc*tsp/b;
c2=c-(c-c2)*tsp/b;
TEo=TEo*tsp/b; 
TEv=TEv*tsp/b;
%scaling for machine units
c=c*8.889;
c2=c2*8.889;
TEo=TEo*8.889; TEv=TEv*8.889;
%dx þarf að ganga upp í c til að loka prófílnum 
dx=c/n;
dx2=c2/n;
%Scaling for equations m=m/100;
m2=m2/100;
p=p/10; 
p2=p2/10; 
t=t/100; 
t2=t2/100;
x=0:dx:c; %Lárétt staðsetning miðað við vænglínu í XY plani 
x2=0:dx2:c2;%Lóðrétt staðsetning miðað við vænglínu í XY plani
%symmetrical airfoil generator
yt=(t*c/0.2).*(0.2969.*(x./c).^0.5-0.1281.*(x./c)- 0.3516.*(x./c).^2+0.2843.*(x./c).^3-0.1015.*(x./c).^4);
yt2=(t2*c2/0.2).*(0.2969.*(x2./c2).^0.5-0.1281.*(x2./c2)- 0.3516.*(x2./c2).^2+0.2843.*(x2./c2).^3-0.1015.*(x2./c2).^4);
%camber line 1 (XY plan)
for i=1:length(x) 
    if x(i)<p*c 
    yc(i)=m.*(x(i)./p^2).*(2*p-x(i)./c) 
    else
    yc(i)=m.*(c-x(i))./(1-p)^2.*(1+x(i)./c-2*p);
    end
end
th(1)=1;
for i=2:length(x)
th(i)=atan((yc(i)-yc(i-1))/dx); %Snertilhorn camber line (XY)
end
%camber line 2 (AB plan)
for i=1:length(x2) 
    if x2(i)<p2*c2
     yc2(i)=m2.*(x2(i)./p2^2).*(2*p2-x2(i)./c2); 
    else
    yc2(i)=m2.*(c2-x2(i))./(1-p2)^2.*(1+x2(i)./c2-2*p2);
    end
end
th2(1)=1;
for i=2:length(x2)
th2(i)=atan((yc2(i)-yc2(i-1))/dx); %Snertilhorn camber line (AB)
end
%points for profile 1 (XY)
xu=x-yt.*sin(th); % X upper
xl=x+yt.*sin(th); % X lower 
yu=yc+yt.*cos(th); % Y upper 
yl=yc-yt.*cos(th); % Y lower
xu=fliplr(xu); %Snúið við til að fá samfella hreyfingu yfir prófílinn
yu=fliplr(yu);

%points for profile 2 (AB)
xu2=x2-yt2.*sin(th2)+100000; %A upper (án vindings)
xl2=x2+yt2.*sin(th2)+100000; %A lower (án vindings)
yu2=yc2+yt2.*cos(th2); %B upper
yl2=yc2-yt2.*cos(th2); %B lower
fu=atan(yu2./xu2); %Vænglínu AB prófíls snúið um hornið inc 
fl=atan(yl2./xl2);
yu2=xu2.*sin(fu+inc)-100000*sin(inc);
yl2=xl2.*sin(fl+inc)-100000*sin(inc); 
xu2=xu2.*cos(fu+inc)-100000*cos(inc);
xl2=xl2.*cos(fl+inc)-100000*cos(inc);
yu2=yu2-yu2(end)+TEv; %Hliðrun A upper í lárétta stefnu
yl2=yl2-yl2(end)+TEv; %Hliðrun A lower í lárétta stefnu 
xu2=xu2+c-xu2(end)-TEo;%Hliðrun B upper í lóðrétta stefnu 
xl2=xl2+c-xl2(end)-TEo;
%Hliðrn B lower í lóðrétta stefnu
xu2=fliplr(xu2); %Snúið við til að fá samfella hreyfingu yfir prófílinn 
yu2=fliplr(yu2);

%Wing plotted in 3D for visualization 
z11=zeros(length(xu))
z12=zeros(length(xl)) 
z21=zeros(length(xu2)) 
z22=zeros(length(xl2)) 
z21=z21+tsp*8.889
z22=z22+tsp*8.889
figure(1) 
hold on 
plot3(xu,yu,z11,'r',xl,yl,z12,'r') 
plot3(xu2,yu2,z21,'r',xl2,yl2,z22,'r') 
for i=1:10:length(xu)
plot3([xu(i) xu2(i)],[yu(i) yu2(i)],[z11(i) z21(i)])
plot3([xl(i) xl2(i)],[yl(i) yl2(i)],[z11(i) z21(i)])
end
axis equal
%grid on
hold off



% GCODE táknbreytir -> breytir línuvigrum í færslulínur G-kóða.
% Fyrst er vigrum breytt úr staðsetningarvigrum yfir í færsuvigra
for i=1:n
xx(i)=xu(i+1)-xu(i);
yy(i)=yu(i+1)-yu(i);

xx(i+length(x))=xl(i+1)-xl(i);
yy(i+length(x))=yl(i+1)-yl(i);
end

for i=1:n
aa(i)=xu2(i+1)-xu2(i);
bb(i)=-1*(yu2(i+1)-yu2(i));
aa(i+length(x2))=xl2(i+1)-xl2(i);
bb(i+length(x2))=-1*(yl2(i+1)-yl2(i));
end
%Svo er skráin opnuð
fid=fopen(filename,'wt');
%Og g-kóðinn skrifaður í hana
fprintf(fid,['G1 A' num2str(-TEo) ' B' num2str(-TEv) ' F3500\n'])
%Hliðrun
for i=1: 2*n+1
fprintf(fid,['G1 X' num2str(xx(i)) ' Y' num2str(yy(i)) ' F4000\n']);
fprintf(fid,['G1 A' num2str(aa(i)) ' B' num2str(bb(i)) ' F4000\n']);

end