%nozzle using the method of characteristics.  
%
%Bryce Roth     10-31-94

%get information from the user
%n=(input('Enter the number of expansion waves. ')+2);
%m=input('Enter the desired test section Mach number. ');
%gamma=input('Enter gamma. ');
n=12;
m=4;
gamma=1.4;

%initialize arrays (matricies, actually) for nu, theta, mu, x, y
nu=zeros(n,n);
theta=zeros(n,n);
mu=zeros(n,n);
x=zeros(n,n);
y=zeros(n,n);

%initialize wall points
nuts=pma(m,gamma);
for i=1:n-2
     theta(1,i+1)=i/(n-2)*nuts/2;
     y(1,i+1)=1;
     nu=theta;
     mu(1,i+1)=180/pi*asin(1/pmm(nu(1,i+1),gamma));
end

%calculate interior points
for i=2:n-1;
     for j=i:n-1;
          if (i==j)                   %isolate axis points
               theta(i,j)=0;
               nu(i,j)=theta(i-1,j)+nu(i-1,j);
               mu(i,j)=180/pi*asin(1/pmm(nu(i,j),gamma));
               y(i,j)=0;
               x(i,j)=-y(i-1,j)/tan((theta(i-1,j)*pi/180-mu(i-1,j)*pi/180+...
                    theta(i,j)*pi/180-mu(i,j)*pi/180)/2)+x(i-1,j);
          else
               theta(i,j)=(theta(i-1,j)+nu(i-1,j)+theta(i,j-1)-...
                    nu(i,j-1))/2;
               nu(i,j)=(theta(i-1,j)+nu(i-1,j)-theta(i,j-1)+...
                    nu(i,j-1))/2;
               mu(i,j)=180/pi*asin(1/pmm(nu(i,j),gamma));
               tan1=tan((theta(i-1,j)*pi/180-mu(i-1,j)*pi/180+...
                    theta(i,j)*pi/180-mu(i,j)*pi/180)/2);
               tan2=tan((theta(i,j-1)*pi/180+mu(i,j-1)*pi/180+...
                    theta(i,j)*pi/180+mu(i,j)*pi/180)/2);
               y(i,j)=(tan1*(-y(i,j-1)/tan2+x(i,j-1)-x(i-1,j))+...
                    y(i-1,j))/(1-tan1/tan2);
               x(i,j)=(y(i,j)-y(i,j-1))/tan2+x(i,j-1);
          end
     end
end

%find wall points
y(1,n)=1;
for i=1:n-1
     theta(i,n)=theta(i,n-1);
     nu(i,n)=nu(i,n-1);
     mu(i,n)=mu(i,n-1);
     if i > 1
          tans=tan((theta(i,n)*pi/180+theta(i-1,n)*pi/180)/2);
          tant=tan(theta(i,n)*pi/180+mu(i,n)*pi/180);
          p1=(y(i,n-1)-tant*x(i,n-1));
          p2=(y(i-1,n)-tans*x(i-1,n));
          x(i,n)=(p1-p2)/(tans-tant);
          y(i,n)=tant*x(i,n)+y(i,n-1)-tant*x(i,n-1);
     end
end

%plot the results
nozx=x(1:n-1,n);
nozy=y(1:n-1,n);
plot(nozx,nozy)


