ardupilot/libraries/SITL/SIM_Plane.cpp

292 lines
9.0 KiB
C++

/*
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
/*
very simple plane simulator class. Not aerodynamically accurate,
just enough to be able to debug control logic for new frame types
*/
#include "SIM_Plane.h"
#include <stdio.h>
using namespace SITL;
Plane::Plane(const char *home_str, const char *frame_str) :
Aircraft(home_str, frame_str)
{
mass = 2.0f;
/*
scaling from motor power to Newtons. Allows the plane to hold
vertically against gravity when the motor is at hover_throttle
*/
thrust_scale = (mass * GRAVITY_MSS) / hover_throttle;
frame_height = 0.1f;
if (strstr(frame_str, "-heavy")) {
mass = 8;
}
if (strstr(frame_str, "-revthrust")) {
reverse_thrust = true;
}
if (strstr(frame_str, "-elevon")) {
elevons = true;
} else if (strstr(frame_str, "-vtail")) {
vtail = true;
}
if (strstr(frame_str, "-elevrev")) {
reverse_elevator_rudder = true;
}
ground_behavior = GROUND_BEHAVIOR_FWD_ONLY;
if (strstr(frame_str, "-ice")) {
ice_engine = true;
}
}
/*
the following functions are from last_letter
https://github.com/Georacer/last_letter/blob/master/last_letter/src/aerodynamicsLib.cpp
many thanks to Georacer!
*/
float Plane::liftCoeff(float alpha) const
{
const float alpha0 = coefficient.alpha_stall;
const float M = coefficient.mcoeff;
const float c_lift_0 = coefficient.c_lift_0;
const float c_lift_a0 = coefficient.c_lift_a;
double sigmoid = ( 1+exp(-M*(alpha-alpha0))+exp(M*(alpha+alpha0)) ) / (1+exp(-M*(alpha-alpha0))) / (1+exp(M*(alpha+alpha0)));
double linear = (1.0-sigmoid) * (c_lift_0 + c_lift_a0*alpha); //Lift at small AoA
double flatPlate = sigmoid*(2*copysign(1,alpha)*pow(sin(alpha),2)*cos(alpha)); //Lift beyond stall
float result = linear+flatPlate;
return result;
}
float Plane::dragCoeff(float alpha) const
{
const float b = coefficient.b;
const float s = coefficient.s;
const float c_drag_p = coefficient.c_drag_p;
const float c_lift_0 = coefficient.c_lift_0;
const float c_lift_a0 = coefficient.c_lift_a;
const float oswald = coefficient.oswald;
double AR = pow(b,2)/s;
double c_drag_a = c_drag_p + pow(c_lift_0+c_lift_a0*alpha,2)/(M_PI*oswald*AR);
return c_drag_a;
}
// Torque calculation function
Vector3f Plane::getTorque(float inputAileron, float inputElevator, float inputRudder, const Vector3f &force) const
{
const float alpha = angle_of_attack;
const float s = coefficient.s;
const float c = coefficient.c;
const float b = coefficient.b;
const float c_l_0 = coefficient.c_l_0;
const float c_l_b = coefficient.c_l_b;
const float c_l_p = coefficient.c_l_p;
const float c_l_r = coefficient.c_l_r;
const float c_l_deltaa = coefficient.c_l_deltaa;
const float c_l_deltar = coefficient.c_l_deltar;
const float c_m_0 = coefficient.c_m_0;
const float c_m_a = coefficient.c_m_a;
const float c_m_q = coefficient.c_m_q;
const float c_m_deltae = coefficient.c_m_deltae;
const float c_n_0 = coefficient.c_n_0;
const float c_n_b = coefficient.c_n_b;
const float c_n_p = coefficient.c_n_p;
const float c_n_r = coefficient.c_n_r;
const float c_n_deltaa = coefficient.c_n_deltaa;
const float c_n_deltar = coefficient.c_n_deltar;
const Vector3f &CGOffset = coefficient.CGOffset;
float rho = air_density;
//read angular rates
double p = gyro.x;
double q = gyro.y;
double r = gyro.z;
//calculate aerodynamic torque
double qbar = 1.0/2.0*rho*pow(airspeed,2)*s; //Calculate dynamic pressure
double la, na, ma;
if (is_zero(airspeed))
{
la = 0;
ma = 0;
na = 0;
}
else
{
la = qbar*b*(c_l_0 + c_l_b*beta + c_l_p*b*p/(2*airspeed) + c_l_r*b*r/(2*airspeed) + c_l_deltaa*inputAileron + c_l_deltar*inputRudder);
ma = qbar*c*(c_m_0 + c_m_a*alpha + c_m_q*c*q/(2*airspeed) + c_m_deltae*inputElevator);
na = qbar*b*(c_n_0 + c_n_b*beta + c_n_p*b*p/(2*airspeed) + c_n_r*b*r/(2*airspeed) + c_n_deltaa*inputAileron + c_n_deltar*inputRudder);
}
// Add torque to to force misalignment with CG
// r x F, where r is the distance from CoG to CoL
la += CGOffset.y * force.z - CGOffset.z * force.y;
ma += -CGOffset.x * force.z + CGOffset.z * force.x;
na += -CGOffset.y * force.x + CGOffset.x * force.y;
return Vector3f(la, ma, na);
}
// Force calculation function from last_letter
Vector3f Plane::getForce(float inputAileron, float inputElevator, float inputRudder) const
{
const float alpha = angle_of_attack;
const float c_drag_q = coefficient.c_drag_q;
const float c_lift_q = coefficient.c_lift_q;
const float s = coefficient.s;
const float c = coefficient.c;
const float b = coefficient.b;
const float c_drag_deltae = coefficient.c_drag_deltae;
const float c_lift_deltae = coefficient.c_lift_deltae;
const float c_y_0 = coefficient.c_y_0;
const float c_y_b = coefficient.c_y_b;
const float c_y_p = coefficient.c_y_p;
const float c_y_r = coefficient.c_y_r;
const float c_y_deltaa = coefficient.c_y_deltaa;
const float c_y_deltar = coefficient.c_y_deltar;
float rho = air_density;
//request lift and drag alpha-coefficients from the corresponding functions
double c_lift_a = liftCoeff(alpha);
double c_drag_a = dragCoeff(alpha);
//convert coefficients to the body frame
double c_x_a = -c_drag_a*cos(alpha)+c_lift_a*sin(alpha);
double c_x_q = -c_drag_q*cos(alpha)+c_lift_q*sin(alpha);
double c_z_a = -c_drag_a*sin(alpha)-c_lift_a*cos(alpha);
double c_z_q = -c_drag_q*sin(alpha)-c_lift_q*cos(alpha);
//read angular rates
double p = gyro.x;
double q = gyro.y;
double r = gyro.z;
//calculate aerodynamic force
double qbar = 1.0/2.0*rho*pow(airspeed,2)*s; //Calculate dynamic pressure
double ax, ay, az;
if (is_zero(airspeed))
{
ax = 0;
ay = 0;
az = 0;
}
else
{
ax = qbar*(c_x_a + c_x_q*c*q/(2*airspeed) - c_drag_deltae*cos(alpha)*fabs(inputElevator) + c_lift_deltae*sin(alpha)*inputElevator);
// split c_x_deltae to include "abs" term
ay = qbar*(c_y_0 + c_y_b*beta + c_y_p*b*p/(2*airspeed) + c_y_r*b*r/(2*airspeed) + c_y_deltaa*inputAileron + c_y_deltar*inputRudder);
az = qbar*(c_z_a + c_z_q*c*q/(2*airspeed) - c_drag_deltae*sin(alpha)*fabs(inputElevator) - c_lift_deltae*cos(alpha)*inputElevator);
// split c_z_deltae to include "abs" term
}
return Vector3f(ax, ay, az);
}
void Plane::calculate_forces(const struct sitl_input &input, Vector3f &rot_accel, Vector3f &body_accel)
{
float aileron = filtered_servo_angle(input, 0);
float elevator = filtered_servo_angle(input, 1);
float rudder = filtered_servo_angle(input, 3);
float throttle;
if (reverse_elevator_rudder) {
elevator = -elevator;
rudder = -rudder;
}
if (elevons) {
// fake an elevon plane
float ch1 = aileron;
float ch2 = elevator;
aileron = (ch2-ch1)/2.0f;
// the minus does away with the need for RC2_REV=-1
elevator = -(ch2+ch1)/2.0f;
} else if (vtail) {
// fake a vtail plane
float ch1 = elevator;
float ch2 = rudder;
// this matches VTAIL_OUTPUT==2
elevator = (ch2-ch1)/2.0f;
rudder = (ch2+ch1)/2.0f;
}
if (reverse_thrust) {
throttle = filtered_servo_angle(input, 2);
} else {
throttle = filtered_servo_range(input, 2);
}
float thrust = throttle;
if (ice_engine) {
thrust = icengine.update(input);
}
// calculate angle of attack
angle_of_attack = atan2f(velocity_air_bf.z, velocity_air_bf.x);
beta = atan2f(velocity_air_bf.y,velocity_air_bf.x);
Vector3f force = getForce(aileron, elevator, rudder);
rot_accel = getTorque(aileron, elevator, rudder, force);
// simulate engine RPM
rpm1 = thrust * 7000;
// scale thrust to newtons
thrust *= thrust_scale;
accel_body = Vector3f(thrust, 0, 0) + force;
accel_body /= mass;
// add some noise
if (thrust_scale > 0) {
add_noise(fabsf(thrust) / thrust_scale);
}
if (on_ground(position)) {
// add some ground friction
Vector3f vel_body = dcm.transposed() * velocity_ef;
accel_body.x -= vel_body.x * 0.3f;
}
}
/*
update the plane simulation by one time step
*/
void Plane::update(const struct sitl_input &input)
{
Vector3f rot_accel;
update_wind(input);
calculate_forces(input, rot_accel, accel_body);
update_dynamics(rot_accel);
// update lat/lon/altitude
update_position();
// update magnetic field
update_mag_field_bf();
}