%test
clear all

% make inner combustion chamber boundary
Rc=3;
sides = 36;
theta = 0;
for i=1:sides+1
   [chamber(i,1),chamber(i,2)] = pol2cart(theta, Rc);
   theta = theta + 360/(sides)*(pi/180);	% radians
end


%====================
Rg_o = .5;
sides = 6;
ports = 4;
dr = .2;
%====================


% make the grain port surface

if ports == 1					% -- single port grain --
	theta = 0;
   
   % make points
   for i=1:sides+1
      [grain_o(i,1,ports),grain_o(i,2,ports)] = pol2cart(theta, Rg_o);
   	theta = theta + 360/(sides)*(pi/180);	% radians
   end

   
elseif ports == 2					% -- double port grain --
   
   	% make top port
      y = (Rc-Rg_o);
      x = sqrt(Rg_o^2-y^2);
      [theta,r] = cart2pol(x,y);
      dtheta = (pi-2*theta)/(sides-1);
      for i=1:sides
         [grain_o(i,1,1),grain_o(i,2,1)] = pol2cart(theta, Rg_o);
         theta = theta + dtheta;
      end
      grain_o(sides+1,1,1) = grain_o(1,1,1);
      grain_o(sides+1,2,1) = grain_o(1,2,1);
      
      % make bottom port
      y = -(Rc-Rg_o);
      x = -sqrt(Rg_o^2-y^2);
      [theta,r] = cart2pol(x,y);
      %dtheta = (pi-2*theta)/(sides-1);
      for i=1:sides
         [grain_o(i,1,2),grain_o(i,2,2)] = pol2cart(theta, Rg_o);
         theta = theta + dtheta;
      end
      grain_o(sides+1,1,2) = grain_o(1,1,2);
      grain_o(sides+1,2,2) = grain_o(1,2,2);
      
elseif ports == 3
   grain_o(1,1,1) = 2;
   grain_o(1,2,1) = 0;
   grain_o(2,1,1) = 1;
   grain_o(2,2,1) = 1;
   grain_o(3,1,1) = 1;
   grain_o(3,2,1) = -1;
   grain_o(4,1,1) = 2;
   grain_o(4,2,1) = 0;
   
   grain_o(1,1,2) = -1;
   grain_o(1,2,2) = 2;
   grain_o(2,1,2) = -2;
   grain_o(2,2,2) = 1;
   grain_o(3,1,2) = -.5;
   grain_o(3,2,2) = 1;
   grain_o(4,1,2) = -1;
   grain_o(4,2,2) = 2;
   
   grain_o(1,1,3) = -1;
   grain_o(1,2,3) = -1;
   grain_o(2,1,3) = -1.5;
   grain_o(2,2,3) = -1.5;
   grain_o(3,1,3) = -.5;
   grain_o(3,2,3) = -2;
   grain_o(4,1,3) = -1;
   grain_o(4,2,3) = -1;
   
elseif ports == 4
	theta = 0;
   
   % make points
   for i=1:sides+1
      [grain_o(i,1),grain_o(i,2)] = pol2cart(theta, Rg_o);
   	theta = theta + 360/(sides)*(pi/180);	% radians
   end
   
   grain_o(:,1,1) = grain_o(:,1) + Rc/2;
   grain_o(:,2,2) = grain_o(:,2) + Rc/2;
   grain_o(:,1,3) = grain_o(:,1) - Rc/2;
   grain_o(:,2,4) = grain_o(:,2) - Rc/2;
   
end %if 
   
   
% plot the unburned grain
figure (1)
plot(chamber(:,1),chamber(:,2),'r')
hold on
for i=1:ports
	plot(grain_o(:,1,i),grain_o(:,2,i),'b')
	hold on
end


% simulate the regression

for k = 1:floor((Rc-Rg_o)/dr)
  
   %make vectors (from the prior point set)
   for h=1:ports
   for i=1:sides
      if i < sides
      	gvect(i,:,h) = [grain_o(i+1,1,h)-grain_o(i,1,h) ...
            grain_o(i+1,2,h)-grain_o(i,2,h) ...
            0];
      	gnorm(i,:,h) = cross(gvect(i,:,h),[0 0 1])/norm(gvect(i,:,h));
   	else
      	gvect(i,:,h) = [grain_o(1,1,h)-grain_o(i,1,h) ...
            grain_o(1,2,h)-grain_o(i,2,h) ...
            0];
      	gnorm(i,:,h) = cross(gvect(i,:,h),[0 0 1])/norm(gvect(i,:,h));         
      end
   end %i
   end %h
  
  
   % make the new points
   for h=1:ports
   j = 1;
   for i=1:sides
      
      grain(j,1,h) = grain_o(i,1,h) + dr*gnorm(i,1,h);
		grain(j,2,h) = grain_o(i,2,h) + dr*gnorm(i,2,h);
      
		if i<sides
   		grain(j+1,1,h) = grain_o(i+1,1,h) + dr*gnorm(i,1,h);
   	   grain(j+1,2,h) = grain_o(i+1,2,h) + dr*gnorm(i,2,h);
   	else
   		grain(j+1,1,h) = grain_o(1,1,h) + dr*gnorm(i,1,h);
   	   grain(j+1,2,h) = grain_o(1,2,h) + dr*gnorm(i,2,h);
      end
      
   	j = j+2;
   end %i
	end %h

   [len_g a b] = size(grain);
   [len_g_o a b] = size(grain_o);
   
   for h=1:ports
   	for i=len_g_o+1:len_g
      	grain_o(i,:,h) = [1 1];
      end   
   end
   
   for h = 1:ports
   	plot(grain(:,1,h),grain(:,2,h))
      hold on;
      grain_o(:,:,h) = grain(:,:,h);
   end
   
	sides = sides*2;

end %k

   
   
hold off