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