% -- This is the Main Function that solves the ODE to get the motion --
% -- Developed to Solve the Flame motion --
% -- Authors: R. Anbarafshan, H. Azizi, R. M. Namin --
% -- Date: 12 / 30 / 2011 --

function freq = func (V)

x = 0;
y = 0.05;
q = 0;          % -- Initial Conditions

N = 5000;       % -- Number of Timesteps
T = 10;         % -- Total Sumilation Time
dt = T / N;

r = [x y q];

x0 = 0.02;

ii = 0;
for i = 1:N
    t = i * dt;
    [tt, rr] = ode45(@(tt, rr) derive(tt, rr, V), [t, t + dt], r);
    
            % -- This Solves the ODE using the differential function
            % -- attached as derive.m using Runge-Kutta 4-5.
    
    r = rr(end, :);
    R(i, :) = r;
    teta(i) = atan(r(1) / r(2));
    
    if(r(1) >= x0)      % -- Discharge Condition
        r(3) = 0;
        ii = ii + 1;
    end
end

freq = double(ii) / double(T);
t = dt:dt:T;
plot(t, R(:, 1:2));