%% Calculating max copterspeed from max thrust and aerodynamic forces
%by Dr. William Thielicke aka Willa aka Shrediquette
clear all;
clc;

rpm=24500; %revs per minute, estimated from electricrcaircraftguy's equations
prop_pitch = 4; %inches;
prop_diameter = 5; %inches;
cd_arm = 1.5; % drag coefficient of the arms (source: Hoerner)
weight = 0.6 * 9.81;% N, weight of the copter

%Three different copter types.
%% Standard racer:
%%{
cd_frame = 1; %Realistic for a bluff body
area_arm = 0.0037 ; %m^2, area of the arms below the propeller          
angle_tilted_body = 0/180*pi; %rad, how much is the nose of the frame tilted up in the design
area_frame_frontal = 0.0025;% m^2 area of the main body of the frame as viewed from the front 
area_frame_dorsal = 0.0107; %m^2 area of the main body of the frame as viewed from the top 
frame_width = 0.05; %m
frame_length = 0.17; %m
%%}

%% QRC5
%{
cd_frame = 0.47;
area_arm = 0.0024 ;           
angle_tilted_body = 30/180*pi; 
area_frame_frontal = 0.004166;
area_frame_dorsal = 0.009377; 
frame_width = 0.069; 
frame_length = 0.154; 
%}

%% DERBE
%{
cd_frame = 1;
area_arm = 0.0015 ;         
angle_tilted_body = 40/180*pi; 
area_frame_frontal = 0.0017;
area_frame_dorsal = 0.0085;
frame_width = 0.05; 
frame_length = 0.18; 
%}

%Initialize variables
angle = 0; %rad, pitch angle of the copter
vertical_force = weight;
while vertical_force >= weight   
    %Initialize variables
    vertical_force = weight;
    copter_drag = 0;
    horizontal_force = 0;
    v=0; % m/s, horizontal speed
    angle = angle+1/180*pi;
    while horizontal_force >= copter_drag
        v=v+0.1/3.6; %slowly increase speed as long as drag is smaller than horizontal force
        v_in_prop = sin(angle)*v; %the relevant component of the propeller disk inflow velocity
        thrust = 4.392399*10^-8*rpm*((prop_diameter^3.5)/sqrt(prop_pitch))*(4.233333*10^-4*rpm*prop_pitch-v_in_prop) * 4; %takern from electricrcaircraftguy
        v_out_prop = 147/3.6; %m/s jet speed of the propeller. Calculated from electricrcaircraftguy
        area_frame_horizontal = abs(sin(angle-angle_tilted_body)*area_frame_dorsal) + abs(cos(angle-angle_tilted_body)*area_frame_frontal); %m^2
        downforce=((frame_width/frame_length)/(frame_width/frame_length+2)) * sin(angle-angle_tilted_body)*area_frame_dorsal*0.5*1.293*v^2; % negative lift is created at positive pitch angles
        armdrag = cd_arm * area_arm * 0.5 * 1.293 * v_out_prop^2; %Drag created by the arms in the propeller jet
        remaining_thrust = thrust - armdrag;
        horizontal_force = sin(angle)*remaining_thrust; %force in the horizontal direction ("thrust")
        vertical_force = cos(angle) * remaining_thrust - downforce; %"lift", must offset weight....
        copter_drag = cd_frame * area_frame_horizontal * 0.5 * 1.293 * v^2; %Dragcoefficient * Area * 0.5 * rho(air) * velocity^2
    end
end
%Output result
disp(['Maximum velocity:  ' num2str(v*3.6) '  km/h'])
disp(['Pitch angle:       ' num2str(angle/pi*180) '   degrees'])