LCOV - code coverage report
Current view: top level - SITL - SIM_Motor.cpp (source / functions) Hit Total Coverage
Test: lcov.info Lines: 96 103 93.2 %
Date: 2026-09-09 11:02:08 Functions: 6 6 100.0 %

          Line data    Source code
       1             : /*
       2             :    This program is free software: you can redistribute it and/or modify
       3             :    it under the terms of the GNU General Public License as published by
       4             :    the Free Software Foundation, either version 3 of the License, or
       5             :    (at your option) any later version.
       6             : 
       7             :    This program is distributed in the hope that it will be useful,
       8             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
       9             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      10             :    GNU General Public License for more details.
      11             : 
      12             :    You should have received a copy of the GNU General Public License
      13             :    along with this program.  If not, see <http://www.gnu.org/licenses/>.
      14             :  */
      15             : /*
      16             :   simple electric motor simulator class
      17             : */
      18             : 
      19             : #include "SIM_Motor.h"
      20             : #include <AP_Motors/AP_Motors.h>
      21             : 
      22             : using namespace SITL;
      23             : 
      24             : // calculate rotational accel and thrust for a motor
      25   452958054 : void Motor::calculate_forces(const struct sitl_input &input,
      26             :                              uint8_t motor_offset,
      27             :                              Vector3f &torque,
      28             :                              Vector3f &thrust,
      29             :                              const Vector3f &velocity_air_bf,
      30             :                              const Vector3f &gyro,
      31             :                              float air_density,
      32             :                              float voltage,
      33             :                              bool use_drag)
      34             : {
      35             : 
      36   452958054 :     const float pwm = input.servos[motor_offset+servo];
      37             :     float command;
      38   452958054 :     float thrust_sign = 1;
      39   452958054 :     if (rsc_servo >= 0) {
      40             :         // variable-pitch rotor; the motor's servo channel commands
      41             :         // blade pitch about mid-PWM and rotor speed comes from the
      42             :         // shared rotor-speed-control servo. The governor is assumed
      43             :         // to hold the rotor at nominal speed from half RSC output;
      44             :         // thrust scales with blade pitch and the square of rotor
      45             :         // speed
      46     2214400 :         const float rotor_speed = constrain_float((input.servos[motor_offset+rsc_servo]-1000)*0.002, 0, 1);
      47     2214400 :         const float collective = constrain_float((pwm-1500)*0.002, -1, 1);
      48     2214400 :         if (is_negative(collective)) {
      49     1244803 :             thrust_sign = -1;
      50             :         }
      51     2214400 :         command = fabsf(collective) * sq(rotor_speed);
      52             :     } else {
      53   450743654 :         command = pwm_to_command(pwm);
      54             :     }
      55   452958054 :     float voltage_scale = voltage / voltage_max;
      56             : 
      57   452958054 :     if (voltage_scale < 0.1) {
      58             :         // battery is dead
      59        4764 :         torque.zero();
      60        4764 :         thrust.zero();
      61        4764 :         current = 0;
      62        4764 :         return;
      63             :     }
      64             : 
      65             :     // apply slew limiter to command
      66   452953290 :     uint64_t now_us = AP_HAL::micros64();
      67   452953290 :     if (last_calc_us != 0 && slew_max > 0) {
      68   452419854 :         float dt = (now_us - last_calc_us)*1.0e-6;
      69   452419854 :         float slew_max_change = slew_max * dt;
      70   452419854 :         command = constrain_float(command, last_command-slew_max_change, last_command+slew_max_change);
      71             :     }
      72   452953290 :     last_calc_us = now_us;
      73   452953290 :     last_command = command;
      74             : 
      75             :     // velocity of motor through air
      76   452953290 :     Vector3f motor_vel = velocity_air_bf;
      77             : 
      78             :     // add velocity of motor about center due to vehicle rotation
      79   452953290 :     motor_vel += -(position % gyro);
      80             : 
      81             :     // calculate velocity into prop, clipping at zero
      82   452953290 :     float velocity_in = MAX(0, -motor_vel.projected(thrust_vector).z);
      83             : 
      84             :     // get thrust for untilted motor
      85   452953290 :     float motor_thrust = thrust_sign * calc_thrust(command, air_density, velocity_in, voltage_scale);
      86             : 
      87             :     // the yaw torque of the motor. A variable-pitch rotor spins in
      88             :     // the same direction regardless of blade pitch sign, so the drag
      89             :     // torque reaction on the frame does not reverse with negative
      90             :     // thrust
      91   452953290 :     const float yaw_scale = 0.05 * diagonal_size * motor_thrust * thrust_sign;
      92   452953290 :     Vector3f rotor_torque = thrust_vector * yaw_factor * command * yaw_scale * -1.0;
      93             : 
      94             :     // thrust in bodyframe NED
      95   452953290 :     thrust = thrust_vector * motor_thrust;
      96             : 
      97             :     // work out roll and pitch of motor relative to it pointing straight up
      98   452953290 :     float roll = 0, pitch = 0;
      99             : 
     100   452953290 :     uint64_t now = AP_HAL::micros64();
     101             :     
     102             :     // possibly roll and/or pitch the motor
     103   452953290 :     if (roll_servo >= 0) {
     104     1277918 :         uint16_t servoval = update_servo(input.servos[roll_servo+motor_offset], now, last_roll_value);
     105     1277918 :         if (roll_min < roll_max) {
     106           0 :             roll = constrain_float(roll_min + (servoval-1000)*0.001*(roll_max-roll_min), roll_min, roll_max);
     107             :         } else {
     108     1277918 :             roll = constrain_float(roll_max + (2000-servoval)*0.001*(roll_min-roll_max), roll_max, roll_min);
     109             :         }
     110             :     }
     111   452953290 :     if (pitch_servo >= 0) {
     112     6622906 :         uint16_t servoval = update_servo(input.servos[pitch_servo+motor_offset], now, last_pitch_value);
     113     6622906 :         if (pitch_min < pitch_max) {
     114           0 :             pitch = constrain_float(pitch_min + (servoval-1000)*0.001*(pitch_max-pitch_min), pitch_min, pitch_max);
     115             :         } else {
     116     6622906 :             pitch = constrain_float(pitch_max + (2000-servoval)*0.001*(pitch_min-pitch_max), pitch_max, pitch_min);
     117             :         }
     118             :     }
     119   452953290 :     last_change_usec = now;
     120             : 
     121             :     // possibly rotate the thrust vector and the rotor torque
     122   452953290 :     if (!is_zero(roll) || !is_zero(pitch)) {
     123     4845874 :         Matrix3f rotation;
     124     4845874 :         rotation.from_euler(radians(roll), radians(pitch), 0);
     125     4845874 :         thrust = rotation * thrust;
     126     4845874 :         rotor_torque = rotation * rotor_torque;
     127             :     }
     128             : 
     129   452953290 :     if (use_drag) {
     130             :         // Momentum drag is modelled as isotropic, scaled by the total
     131             :         // thrust, which handles tilted motors correctly. The component
     132             :         // along the rotor axis is really a 'hack' to compensate for
     133             :         // incorrect modelling of the variation of thrust with inflow
     134             :         // velocity. If not applied, the vehicle will climb at an
     135             :         // unrealistic rate during operation in STABILIZE. TODO replace
     136             :         // prop and motor model with one based on DC motor, momentum
     137             :         // disc and blade element theory.
     138   452953290 :         const float momentum_drag_factor = momentum_drag_coefficient * sqrtf(air_density * true_prop_area * thrust.length());
     139   452953290 :         thrust -= motor_vel * momentum_drag_factor;
     140             :     }
     141             : 
     142             :     // calculate total torque in newton-meters
     143   452953290 :     torque = (position % thrust) + rotor_torque;
     144             : 
     145             :     // calculate current
     146   452953290 :     float power = power_factor * fabsf(motor_thrust);
     147   452953290 :     current = power / MAX(voltage, 0.1);
     148             : }
     149             : 
     150             : /*
     151             :   update and return current value of a servo. Calculated as 1000..2000
     152             :  */
     153     7900824 : uint16_t Motor::update_servo(uint16_t demand, uint64_t time_usec, float &last_value) const
     154             : {
     155     7900824 :     if (servo_rate <= 0) {
     156           0 :         return demand;
     157             :     }
     158     7900824 :     if (servo_type == SERVO_RETRACT) {
     159             :         // handle retract servos
     160      731082 :         if (demand > 1700) {
     161      333956 :             demand = 2000;
     162      397126 :         } else if (demand < 1300) {
     163      397126 :             demand = 1000;
     164             :         } else {
     165           0 :             demand = last_value;
     166             :         }
     167             :     }
     168     7900824 :     demand = constrain_int16(demand, 1000, 2000);
     169     7900824 :     float dt = (time_usec - last_change_usec) * 1.0e-6f;
     170             :     // assume servo moves through 90 degrees over 1000 to 2000
     171     7900824 :     float max_change = 1000 * (dt / servo_rate) * 60.0f / 90.0f;
     172     7900824 :     last_value = constrain_float(demand, last_value-max_change, last_value+max_change);
     173     7900824 :     return uint16_t(last_value+0.5);
     174             : }
     175             : 
     176             : 
     177             : // calculate current and voltage
     178   452958054 : float Motor::get_current(void) const
     179             : {
     180   452958054 :     return current;
     181             : }
     182             : 
     183             : // setup PWM ranges for this motor
     184   452962818 : void Motor::setup_params(uint16_t _pwm_min, uint16_t _pwm_max, float _spin_min, float _spin_max, float _expo, float _slew_max,
     185             :                          float _diagonal_size, float _power_factor, float _voltage_max, float _effective_prop_area,
     186             :                          float _velocity_max, Vector3f _position, Vector3f _thrust_vector, float _yaw_factor, 
     187             :                          float _true_prop_area, float _momentum_drag_coefficient)
     188             : {
     189   452962818 :     mot_pwm_min = _pwm_min;
     190   452962818 :     mot_pwm_max = _pwm_max;
     191   452962818 :     mot_spin_min = _spin_min;
     192   452962818 :     mot_spin_max = _spin_max;
     193   452962818 :     mot_expo = _expo;
     194   452962818 :     slew_max = _slew_max;
     195   452962818 :     power_factor = _power_factor;
     196   452962818 :     voltage_max = _voltage_max;
     197   452962818 :     effective_prop_area = _effective_prop_area;
     198   452962818 :     max_outflow_velocity = _velocity_max;
     199   452962818 :     true_prop_area = _true_prop_area;
     200   452962818 :     momentum_drag_coefficient = _momentum_drag_coefficient;
     201   452962818 :     diagonal_size = _diagonal_size;
     202             : 
     203   452962818 :     if (!_position.is_zero()) {
     204           0 :         position = _position;
     205             :     } else {
     206   452962818 :         position.x = cosf(radians(angle)) * _diagonal_size;
     207   452962818 :         position.y =  sinf(radians(angle)) * _diagonal_size;
     208   452962818 :         position.z = 0;
     209             :     }
     210             : 
     211   452962818 :     if (!_thrust_vector.is_zero()) {
     212           0 :         thrust_vector = _thrust_vector;
     213             :     }
     214   452962818 :     if (!is_zero(_yaw_factor)) {
     215           0 :         yaw_factor = _yaw_factor;
     216             :     }
     217   452962818 : }
     218             : 
     219             : /*
     220             :   convert a PWM value to a command value from 0 to 1
     221             : */
     222   450743654 : float Motor::pwm_to_command(float pwm) const
     223             : {
     224   450743654 :     const float pwm_thrust_max = mot_pwm_min + mot_spin_max * (mot_pwm_max - mot_pwm_min);
     225   450743654 :     const float pwm_thrust_min = mot_pwm_min + mot_spin_min * (mot_pwm_max - mot_pwm_min);
     226   450743654 :     const float pwm_thrust_range = pwm_thrust_max - pwm_thrust_min;
     227   450743654 :     return constrain_float((pwm-pwm_thrust_min)/pwm_thrust_range, 0, 1);
     228             : }
     229             : 
     230             : /*
     231             :   calculate thrust given a command value
     232             : */
     233   452953290 : float Motor::calc_thrust(float command, float air_density, float velocity_in, float voltage_scale) const
     234             : {
     235   452953290 :     float velocity_out = voltage_scale * max_outflow_velocity * sqrtf((1-mot_expo)*command + mot_expo*sq(command));
     236   452953290 :     float ret = 0.5 * air_density * effective_prop_area * (sq(velocity_out) - sq(velocity_in));
     237             : #if 0
     238             :     if (command > 0) {
     239             :         ::printf("air_density=%f effective_prop_area=%f velocity_in=%f velocity_max=%f\n",
     240             :                  air_density, effective_prop_area, velocity_in, voltage_scale * max_outflow_velocity);
     241             :         ::printf("calc_thrust %.3f %.3f\n", command, ret);
     242             :     }
     243             : #endif
     244   452953290 :     return ret;
     245             : }

Generated by: LCOV version 1.14