55

FirePath

A code base that generates thrust data for a hybrid rocket engine using mathematical models and the Chemical Equilibrium with Applications (CEA) tool, HRETS-CEA

The entire codebase in Matlab - Codebase for calculating the area and perimeter of a variable non-convex polygons and animates burn Working Progress - CEA Implimentation

clf
clear 
Radius = input('Radius of Fueslage:');
t = linspace(0,2*pi,1000);
plot(Radius*cos(t),Radius*sin(t),'b')
axis equal
 
% point input
loop = 0;
while loop < 2
    loop = input('How many points :');
end
disp('Pick points inside the circle to make a non - convex polgon')
temp = zeros(2,loop);
hold on
cord = ones(2,(19*loop));
for i = 1:loop
     temp(:,i) =[ 1000000000; 1000000000];
     while (temp(1,i)^2 + temp(2,i)^2 > Radius^2)
         [temp(1,i),temp(2,i)] = ginput(1);
         if (temp(1,i)^2 + temp(2,i)^2 > Radius^2)
             disp('Point shold be inside circle')
         end
     end
     plot(temp(1,i),temp(2,i),'ro')
     text(temp(1,i),temp(2,i),int2str(i),'FontSize',14)
end
 
%ordering the points 
if loop > 3
    x_m = mean(temp(1,:));
    y_m = mean(temp(2,:));
    d = zeros(loop,loop);
 
    for i = 1:loop
        for j = 1:loop
            if j ~= i
            d(i,j) = mid_dist(temp(1,i),temp(2,i),temp(1,j),temp(2,j),x_m,y_m);
            end
        end
    end
 
    %creating point map based on farthest distance of line from centroid
    pointmap = 1:1:loop;
    i = 1;
    [m ,pointmap(i)] = max(d(:,i));
    d(i,:) =  zeros(1,loop);
    i = pointmap(i);
    while i ~= 1
        [m ,pointmap(i)] = max(d(:,i));
        d(i,:) =  zeros(1,loop);
        i = pointmap(i);
    end
 
    temp1 = temp;
 
    % creating lines from the points
    i = pointmap(1);
    for k = 2:1:loop
        temp(1,k) = temp1(1,i); 
        temp(2,k) = temp1(2,i); 
        i = pointmap(i);
    end
 
end
 
 
% creating lines from the points 
for i = 1:loop
     if i > 1 
        segment_x = linspace(temp(1,i-1),temp(1,i),20);
        segment_y = linspace(temp(2,i-1),temp(2,i),20);
     else
        segment_x = linspace(temp(1,loop),temp(1,1),20);
        segment_y = linspace(temp(2,loop),temp(2,1),20);
     end 
        segment_x = segment_x(2:end);
        segment_y = segment_y(2:end);
        cord(1,((19*(i-1)+1):19*i)) = segment_x;
        cord(2,((19*(i-1)+1):19*i)) = segment_y;
end
 
temp = [temp temp(:,1)];
hold on
% plot(temp(1,:),temp(2,:),'r')
plot(cord(1,:),cord(2,:),'go')
 
%Burn code
[t_b,O] = burn_time(cord,0.1,Radius);
plot(O(1),O(2),'k*');
r = 0.1 * ones(1,t_b);
[x,y] = burn(cord,r,Radius,t_b);
 
%plot cord
h = plot(x(1,:),y(1,:),'go');
axis equal
 
for j = 1:t_b
    set(h,'XData',x(j,:));
    set(h,'YData',y(j,:));
    drawnow;
    pause(0.075)
end
 
function[d] = mid_dist(x_1,y_1,x_2,y_2,x_0,y_0)
   x_m = (x_1 + x_2)/2;
   y_m = (y_1 + y_2)/2;
   d =  sqrt((x_0 - x_m)^2 + (y_0 - y_m)^2);
end
 
function[d] = perp_dist(x_1,y_1,x_2,y_2,x_0,y_0)
   d =  abs(((x_2 - x_1)*(y_1 - y_0)  - (x_1 - x_0)*(y_2 - y_1))/(sqrt((x_2 - x_1)^2 + (y_2 - y_1)^2))) ;
end
 
function [P] = intersection(a_1,b_1,a_d,b_d,r)
 syms y_1 x_1 x_d y_d m y x 
if a_1 ==  a_d
    P(1) = a_1;
    temp_y = sqrt(r^2 - a_1^2 );
    if temp_y > b_d
        temp_y = -temp_y;
    end
    P(2) = temp_y;
else 
 
    
    m = (y_d - y_1)/(x_d - x_1);
    eqn = (y_d - y)/(x_d - x) == m;
    S = solve(eqn,y);
    eq2 = x^2 + S^2 == r^2;
    K = solve(eq2,x);
    x_1 = a_1;
    y_1 = b_1;
    x_d = a_d;
    y_d = b_d;
    L =  subs(K);
    M = subs(S);
    
    
    if L(1) > a_1 &&  L(1) < a_d
        P(1) = L(1);
    else
        P(1) = L(2);
    end
    
    P(2) = subs(M,x,P(1));
end
end
 
function[ t,c ] = burn_time(cord,r_b,R)
x_cord = [cord(1,:)];
y_cord = [cord(2,:)];
x_m = mean(x_cord);
y_m = mean(y_cord);
 
x_cord_new = x_cord;
y_cord_new = y_cord;
 
k = size(x_cord);
n = (x_cord.^2 + y_cord.^2).^(0.5);
 
[M,index] = min(n);
 
i = index;
 
if (i > 1) && (i < k(2)-1)
    tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(i+1) - x_cord(i)))/2....
        , -((y_cord(i) - y_cord(i-1)) + (y_cord(i+1) - y_cord(i)))/2  ] ;
elseif i == 1
    tangent_dir = [((x_cord(i) - x_cord(k(2))) + (x_cord(i+1) - x_cord(i)))/2....
        , -((y_cord(i) - y_cord(k(2))) + (y_cord(i+1) - y_cord(i)))/2  ] ;
elseif i == k(2)
    tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(1) - x_cord(i)))/2....
        , -((y_cord(i) - y_cord(i-1)) + (y_cord(1) - y_cord(i)))/2  ] ;
end 
 
 
temp_burn_dir = flip(tangent_dir) / norm(tangent_dir);
 
point = [x_cord(i) - x_m , y_cord(i) - y_m] ;
new_point  = point + temp_burn_dir ; 
 
if norm(new_point) < norm(point)
   burn_dir = -temp_burn_dir;
else
   burn_dir = temp_burn_dir;
end
 
x_cord_new(i) = r_b*burn_dir(1) + x_cord(i);
y_cord_new(i) = r_b*burn_dir(2) + y_cord(i);
 
t = 1;
 
while x_cord_new(i)^2 + y_cord_new(i)^2 < R^2
    x_cord_new(i) = r_b*burn_dir(1) + x_cord_new(i);
    y_cord_new(i) = r_b*burn_dir(2) + y_cord_new(i);
    t = t+1;
end
t = t+10;
c = [x_cord(i) ; y_cord(i)];
end
 
function[x,y] = burn(cord,r_b,R,t_b)
 
x_cord = [cord(1,:)];
y_cord = [cord(2,:)];
plot(x_cord,y_cord,'ro')
 
x_cord_new = x_cord;
y_cord_new = y_cord;
k = size(x_cord);
x = ones(t_b,k(2));
y = ones(t_b,k(2));
x(1,:) = x_cord;
y(1,:) = y_cord;
end_flag = zeros(1,k(2));
 
 
for j = 1:t_b
    
    for i = 1:k(2)
        if end_flag(i) == 0 
            if (i > 1) && (i < k(2)-1)
                tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(i+1) - x_cord(i)))/2....
                    , -((y_cord(i) - y_cord(i-1)) + (y_cord(i+1) - y_cord(i)))/2  ] ;
            elseif i == 1
                tangent_dir = [((x_cord(i) - x_cord(k(2))) + (x_cord(i+1) - x_cord(i)))/2....
                    , -((y_cord(i) - y_cord(k(2))) + (y_cord(i+1) - y_cord(i)))/2  ] ;
            elseif i == k(2)
                tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(1) - x_cord(i)))/2....
                    , -((y_cord(i) - y_cord(i-1)) + (y_cord(1) - y_cord(i)))/2  ] ;
            end 
            temp_burn_dir = r_b(j)*flip(tangent_dir) / norm(tangent_dir);
            
            point = [x_cord(i) - mean(x_cord), y_cord(i)-mean(y_cord)] ;
            new_point  = point + temp_burn_dir ; 
            if norm(new_point) < norm(point)
               burn_dir = -temp_burn_dir;
            else
               burn_dir = temp_burn_dir;
            end
            
            x_cord_new(i) = burn_dir(1)+x_cord(i);
            y_cord_new(i) = burn_dir(2)+y_cord(i);
    
             if sqrt(x_cord_new(i)^2 + y_cord_new(i)^2) > R
                 K = intersection(x_cord(i),y_cord(i),x_cord_new(i),y_cord_new(i),R);
                 x_cord_new(i) = K(1);
                 y_cord_new(i) = K(2);
                 end_flag(i) = 1;
             end
        else 
            x_cord_new(i) = x_cord(i);
            y_cord_new(i) = y_cord(i);
        end
    end
    x_cord = x_cord_new;
    y_cord = y_cord_new;
    x(j,:) = x_cord_new;
    y(j,:) = y_cord_new;
 
end
end