%Author Terry Cornall with advice and feedback from Dr. Andrew Price.
%Copyright 2006 Terry Cornall.
%This code is for use with a hot-wire foam cutter. It does a calculation
%that allows a wing to be cut that would otherwise not fit on the table, by rotating the wing
%so that the trailing edge is parallel to the face of the block of foam.
%This requires that the wing profile as used by the foam cutter be
%re-calculated so that when the real root and tip chords are cut at an angle into the
%foam wing then the correct profile is restored.
%Reads the named airfoil (profile) file and creates a wing surface
%description that is then rotated to make the trailing edge horizontal
%to a foam block for cutting.
%Produces a new airfoil file to be used by the foam cutter.
%The new airfoil is the projection of the original airfoil (after it is
%rotated) onto the z=0 plane
%The input is a file called 'mh61-1.cor'
%(Change that to suit your own needs)
%The output is a file called 'mh61_rotated.cor'. It contains the new profile
%and also the new root and tip chords and the span and the angle that the wing is rotated by.
%Also produces some plots of rotated and unrotated wing and compares the
%new and original profiles so you can do sanity checks.
%assumes that the first point in the chord file is the trailing edge and
%that the leading edge of the root chord is at 0,0. Assumes that the
%maximum (in x dimension) point of the root chord is the (normalised)trailing edge
%Assumes there is no twist (i.e. washout) in the wing.
%assumes a swept-back wing.
%assumes tip profile is same as root profile except scaled down.
%Does calculate an x,y,z point for tip and chord leading edges but doesn't
%insert those points in the rotated profile, rather it does a transform on the points
%in the original profile.
%Although the results of this code have been checked agains a CAD package,
%there is no guarantee that there are no mistakes or miscalculations
%herein,
%You are welcome to use and change the code at your own risk, for non-profit purposes.
filename='mh61-1.cor'
outputfilename='mh61_rotated.cor'
span=1000; %original wing span in mm
root_chord=500; %length of root
tip_chord=250; %length of tip
sweep_deg=27.5; %amount of sweep on the leading edge of the wing, in degrees (+ve is sweep back)
sweep=pi/180*sweep_deg; %sweep in radians becaue Matlab uses radians in sin, cos and tan etc
%span=R*cos(sweep); %length of original leading edge
R=span/cos(sweep);
dx_chords=R*sin(sweep); %distance back of the original tip LE from the original root LE
tip_position=[dx_chords;0]; %position of original tip leading edge
%The tan of the angle between the trailing edge and horizontal is the ratio of the
%(the difference between the trailing edge of root and tip )over the span.
%tan(rot_angle)=(dx_chords+tip_chord-root_chord)/span;
rot_angle=-atan(((tip_position(1)+tip_chord)-(root_chord))/span) ;
rot_angle_deg=rot_angle*180/pi;
%this is the angle to rotate the wing by to get the leading edge horizontal and aligned to the foam block
infile=fopen(filename,'r'); %open input airfoil file
outfile=fopen(outputfilename,'w'); %open output airfoil file
if(infile==-1)
['Cant open ',filename]
return;
end
if(outfile==-1)
['Cant open ',outputfilename]
return;
end
n=0;
s=fgets(infile) %read first line whcih we will treat as a title
n=1;
while(~feof(infile)) %until finished
s=fgets(infile); %read a line
[parsed,count]=sscanf(s,'%f %f'); %get x and y points
points(1,n)=parsed(1);
points(2,n)=parsed(2);
n=n+1; %keep track of number of (data) lines read
end
fclose(infile); %finished with input
root_profile=points*root_chord;
tip_profile=points*tip_chord; %scale to mm
%make 3d. x axis is from leading edge to trailing edge.
%y axis is perpendicular to wing.
%z axis is along the length of the wing.
root_3d(1,:)=root_profile(1,:);
root_3d(2,:)=root_profile(2,:);
root_3d(3,:)=0; %make an x,y,z version of the root profile and set it at z=0
tip_3d(1,:)=tip_profile(1,:)+tip_position(1);
tip_3d(2,:)=tip_profile(2,:)+tip_position(2);
tip_3d(3,:)=span; %and put the tip at z=span
plot3(root_3d(1,:),root_3d(2,:),root_3d(3,:),'b'); %plot the root and tip in a 3d plot
hold on;
plot3(tip_3d(1,:),tip_3d(2,:),tip_3d(3,:),'b');
% line([0,tip_position(1)],[0,tip_position(2)],[0,span]); %leading edge
% line([0+root_chord,tip_position(1)+tip_chord],[0,tip_position(2)],[0,span]); %trailing edge
%Now I assume that the points in the profiles are arranged
%Trailing edge first point, leading edge middle of the set, trailing edge last point and that's the way this
%.cor file seems to be
l=line([root_3d(1,1),tip_3d(1,1)],[root_3d(2,1),tip_3d(2,1)],[root_3d(3,1),tip_3d(3,1)]); %trailingedge
set(l,'Color',[0,1,1]); %make it greeny-blue to show that we know our trailing edge from our leading edge
%draw the leading edge of the original unrotated wing
%The position of the leading edge of the root is 0,0,0.
%The position of the leading edge of the tip is the trailing edge of the tip-the tip chord for x, 0 for y, plus the span for z
l=line([0,tip_3d(1,1)-tip_chord],[0,0],[0,span]); %leading edge of tip
set(l,'Color',[0,0,1]); %make it blue
%axis([0,root_chord,-0.5*root_chord,0.5*root_chord,0,span]); %set sizes of axes
view(-6.5, 52); %give us a nice slightly rotated view of the 3d data
%OK, so far the code has been for verification and setup. Now for the work.
%Rotate the 3d datapoints for the root and tip profiles so that the
%trailing edge is aligned with the z axis. The rotation is in the x,z plane,
%i.e. around the y axis, so y values stay the same
rot=[cos(rot_angle),0,sin(rot_angle);
0 ,1 ,0; %rotation operator for rotation in x,z plane
-sin(rot_angle),0,cos(rot_angle)];
rot_root_3d=rot*root_3d; %apply the rotation opearator
rot_tip_3d=rot*tip_3d; % to the points on root and tip profiles
plot3(rot_root_3d(1,:),rot_root_3d(2,:),rot_root_3d(3,:),'g'); %plot the rotated root profile
plot3(rot_tip_3d(1,:),rot_tip_3d(2,:),rot_tip_3d(3,:),'g'); %plot the rotated tip profile
l=line([rot_root_3d(1,1),rot_tip_3d(1,1)],[rot_root_3d(2,1),rot_tip_3d(2,1)],[rot_root_3d(3,1),rot_tip_3d(3,1)]); %rotated trailing edge
set(l,'Color',[1,0,1]); %make it purple to show that we know our trailing edge from our leading edge
%Now draw the leading edge of the rotated wing
%Leading edge of root is still 0,0,0 and leading edge of tip is now rot*[tip_3d(1,1)-tip_chord;0;span]
rot_tip_3d_le=rot*[tip_3d(1,1)-tip_chord;0;span];
l=line([0,rot_tip_3d_le(1)],[0,rot_tip_3d_le(2)],[0,rot_tip_3d_le(3)]); %rotated leading edge
set(l,'Color',[0,1,0]); %make it green
axis ('equal'); %and make aspect ratio 1:1:1
title('Original (blue) and rotated (green) wing, with projected profiles (red)');
%now work out the new span. It'll be from the rotated tip TE to the origin,
%in the z dimension.
new_span=abs(rot_tip_3d(3,1)); %TE is at point 1 in array
%Now for the real work. Work out what the projections along the lines from the tip profile of the rotated root
%profile onto the z=0 plane are. This is sort of like ray-tracing and I take
%a similar approach from Heckbert's 'Graphic Gems IV' published by AP Professional in 1994.
%Knowing 2 points of the ray
%i.e. a rotated root point rr and corresponding rotated tip profile point rt,
%then we can write a vector equation for a ray on the rotated wing surface from tip to root as
% p=rr+t*(rt-rr) where p is an x,y,z point on the ray
%and t is a parametric scalar
%and then , writing the the equation of the
%plane Z=0 as k dot (p-p0)=0 where k is the unit vector in z direction, [0,0,1] and dot is
%inner product (dot product) operator,
%p is a point [x,y,0] on the z=0 plane and p0 is any other point in the z=0 plane
%Make p0= [0,0,0] so equation of plane becomes
% k dot p=0
%, then the equation for p, the point of intersection of ray with z=0 plane is given by
% k dot (rr+t*(rt-rr)=0; where all the values are known except t which can then be found
% using
% k.rr= -t*(k.rt-rr) => t= -(k.rr)/(k.(rt-rr))
%and plugged back into p=rr+t*(rt-rr) to find the desired point p.
rr=rot_root_3d; rt=rot_tip_3d; %shorten names of profile points arrays
%t= -(k dot rr)/(k dot (rt-rr))
%k dot rr is just the z component of rr, i.e. rr(3) but rr is array of points, so need rr(3,:)
%could use -([0,0,1] * rr )./ ([0,0,1]* (rt-rr)) but it's the same thing.
g=rt-rr; %g is an array of vectors representing the lines between the root and tip
t= -rr(3,:) ./ g(3,:); % ./ is needed rather than / because we have an array
%of results and we just want element by element division rather
%than vector division
%t represents the scale of extension of each ray needed to hit the z=0 plane
for(n=1:length(t)) %have to do the * in a loop because t is nx[1x1] but g is nx[3x1]
p(:,n)=rr(:,n)+t(n)*g(:,n); %and matlab has no way to do it concurrently in that case
end %that I know of anyway. Anybody?
plot3(p(1,:),p(2,:),p(3,:),'r'); %plot the projection of the rotated root profile onto the z=0 plane
s=['Angle of rotation is ',num2str(rot_angle_deg,3)];
text(100,100,s);
%now have to do the projection of tip profile onto z=new_span plane
%k dot (p-p0)=0 becomes k dot (pt-pns) and choose
%pns=[0,0,new_span] and pt is the projection of rotated tip
%profile. Again, k=[0,0,1] so just get the z component of rt-pns
%t= -(k dot (rt-pns))/(k dot (rt-rr)) becomes
%t= -(rt-pns)(3,:) ./ (rt-rr)(3,:) (except matlab won't let you write
%it that way.
m=rt; %want to get m=rt-pns;
m(3,:)=rt(3,:)-new_span; % already have g=rt-rr from above
t= -m(3,:) ./ g(3,:); % amount to stretch ray to hit z=new_span plane
for(n=1:length(t)) %have to do the * in a loop because t is nx[1x1] but g is nx[3x1]
pt(:,n)=rt(:,n)+t(n)*g(:,n); %and matlab has no way to do it concurrently in that case
end %that I know of anyway. Anybody?
plot3(pt(1,:),pt(2,:),pt(3,:),'r'); %plot the projection of the rotated tip profile onto the z=new_span plane
%similarly, need to do the projection of the tip leading edge (not
%necessarily in the profile)
m=rot_tip_3d_le;
m(3)=rot_tip_3d_le(3)-new_span;
g=rot_tip_3d_le-0;
t=-m(3) / g(3);
pt_le=rot_tip_3d_le+t*g;
new_tip_chord=abs(pt_le(1)-pt(1,1)); %ta-da, the end of required calculations...
%now normalise the root profile
new_root_chord=abs(p(1,1)); %difference between trailing edge and leading edge of the root profile
p=p/new_root_chord; %normalise in the x direction and shift to zero if needed
figure;
plot(p(1,:),p(2,:),'r'); %plot new chord
hold;
plot(points(1,:),points(2,:),'b'); %compare with original chord
axis('equal');
legend('normalised projection of rotated profile','original profile');
%write to file, with explanatory header and other parameters
fprintf(outfile,'%s with sweep of %0.2f rotated by %0.2f and new cutting profile calculated and normalised. New span=%0.2f, new root chord=%0.2f new tip chord=%0.2f\n',filename,sweep_deg,rot_angle_deg,new_span,new_root_chord, new_tip_chord);
for(n=1:length(p))
fprintf(outfile,'%0.4f %0.4f\n',p(1,n),p(2,n)); %write out the new profile to the output file
end
fclose(outfile);