1 #include <simgear/debug/logstream.hxx>
4 #include <Main/fg_props.hxx>
6 #include "Rotorpart.hpp"
14 SG_USING_STD(setprecision);
27 const float pi=3.14159;
29 static inline float sqr(float a) { return a * a; }
46 _base[0] = _base[1] = _base[2] = 0;
55 _forward[1]=_forward[2]=0;
56 _max_pitch=13./180*pi;
57 _maxcyclicail=10./180*pi;
58 _maxcyclicele=10./180*pi;
60 _mincyclicail=-10./180*pi;
61 _mincyclicele=-10./180*pi;
62 _min_pitch=-.5/180*pi;
66 _normal[0] = _normal[1] = 0;
68 _normal_with_yaw_roll[0]= _normal_with_yaw_roll[1]=0;
69 _normal_with_yaw_roll[2]=1;
71 _omega=_omegan=_omegarel=_omegarelneu=0;
81 _shared_flap_hinge=false;
82 _rellenteeterhinge=0.01;
89 _translift_maxfactor=1.3;
90 _ground_effect_constant=0.1;
91 _vortex_state_lift_factor=0.4;
104 _number_of_segments=1;
106 _rel_len_where_incidence_is_measured=0.7;
107 _torque_of_inertia=1;
111 _airfoil_incidence_no_lift=0;
112 _rel_len_blade_start=0;
113 _airfoil_lift_coefficient=0;
114 _airfoil_drag_coefficient0=0;
115 _airfoil_drag_coefficient1=0;
117 _global_ground[i] = _tilt_center[i] = 0;
118 _global_ground[2] = 1;
119 _global_ground[3] = -1e3;
120 _incidence_stall_zero_speed=18*pi/180.;
121 _incidence_stall_half_sonic_speed=14*pi/180.;
122 _lift_factor_stall=0.28;
123 _stall_change_over=2*pi/180.;
124 _drag_factor_stall=8;
129 for (int k=0;k<4;k++)
131 _groundeffectpos[k][i]=0;
132 _ground_effect_altitude=1;
134 _lift_factor=_f_ge=_f_vs=_f_tl=1;
135 _rotor_correction_factor=.65;
139 _num_ground_contact_pos=0;
140 _directions_and_postions_dirty=true;
159 for(i=0; i<_rotorparts.size(); i++) {
160 Rotorpart* r = (Rotorpart*)_rotorparts.get(i);
163 //untie the properties
166 SGPropertyNode * node = fgGetNode("/rotors", true)->getNode(_name,true);
167 node->untie("balance-ext");
168 node->untie("balance-int");
173 void Rotor::inititeration(float dt,float omegarel,float ddt_omegarel,float *rot)
178 _omega=_omegan*_omegarel;
179 _ddt_omega=_omegan*ddt_omegarel;
182 updateDirectionsAndPositions(drot);
183 Math::add3(rot,drot,drot);
184 for(i=0; i<_rotorparts.size(); i++) {
185 float s = Math::sin(float(2*pi*i/_number_of_parts+(_phi-pi/2.)*(_ccw?1:-1)));
186 float c = Math::cos(float(2*pi*i/_number_of_parts+(_phi-pi/2.)*(_ccw?1:-1)));
187 Rotorpart* r = (Rotorpart*)_rotorparts.get(i);
189 r->setDdtOmega(_ddt_omega);
190 r->inititeration(dt,drot);
191 r->setCyclic(_cyclicail*c+_cyclicele*s);
194 //calculate the normal of the rotor disc, for calcualtion of the downwash
195 float side[3],help[3];
196 Math::cross3(_normal,_forward,side);
197 Math::mul3(Math::cos(_yaw)*Math::cos(_roll),_normal,_normal_with_yaw_roll);
199 Math::mul3(Math::sin(_yaw),_forward,help);
200 Math::add3(_normal_with_yaw_roll,help,_normal_with_yaw_roll);
202 Math::mul3(Math::sin(_roll),side,help);
203 Math::add3(_normal_with_yaw_roll,help,_normal_with_yaw_roll);
206 if ((_balance1*_balance2 < 0.97) && (_balance1>-1))
208 _balance1-=(0.97-_balance1*_balance2)*(0.97-_balance1*_balance2)*0.005;
209 if (_balance1<-1) _balance1=-1;
213 float Rotor::calcStall(float incidence,float speed)
215 float stall_incidence=_incidence_stall_zero_speed
216 +(_incidence_stall_half_sonic_speed
217 -_incidence_stall_zero_speed)*speed/(343./2);
218 //missing: Temeperature dependency of sonic speed
219 incidence = Math::abs(incidence);
220 if (incidence > (90./180.*pi))
221 incidence = pi-incidence;
222 float stall = (incidence-stall_incidence)/_stall_change_over;
223 stall = Math::clamp(stall,0,1);
225 _stall_sum+=stall*speed*speed;
226 _stall_v2sum+=speed*speed;
231 float Rotor::getLiftCoef(float incidence,float speed)
233 float stall=calcStall(incidence,speed);
234 /* the next shold look like this, but this is the inner loop of
235 the rotor simulation. For small angles (and we hav only small
236 angles) the first order approximation works well
237 float c1= Math::sin(incidence-_airfoil_incidence_no_lift)*_liftcoef;
238 for c2 we would need higher order, because in stall the angle can be large
241 if (incidence > (pi/2))
243 else if (incidence <-(pi/2))
247 float c1= (i2-_airfoil_incidence_no_lift)*_liftcoef;
250 float c2= Math::sin(2*(incidence-_airfoil_incidence_no_lift))
251 *_liftcoef*_lift_factor_stall;
252 return (1-stall)*c1 + stall *c2;
258 float Rotor::getDragCoef(float incidence,float speed)
260 float stall=calcStall(incidence,speed);
261 float c1= (Math::abs(Math::sin(incidence-_airfoil_incidence_no_lift))
262 *_dragcoef1+_dragcoef0);
263 float c2= c1*_drag_factor_stall;
264 return (1-stall)*c1 + stall *c2;
267 int Rotor::getValueforFGSet(int j,char *text,float *f)
269 if (_name[0]==0) return 0;
270 if (4>numRotorparts()) return 0; //compile first!
273 sprintf(text,"/rotors/%s/cone-deg", _name);
274 *f=(_balance1>-1)?( ((Rotorpart*)getRotorpart(0))->getrealAlpha()
275 +((Rotorpart*)getRotorpart(1*(_number_of_parts>>2)))->getrealAlpha()
276 +((Rotorpart*)getRotorpart(2*(_number_of_parts>>2)))->getrealAlpha()
277 +((Rotorpart*)getRotorpart(3*(_number_of_parts>>2)))->getrealAlpha()
283 sprintf(text,"/rotors/%s/roll-deg", _name);
284 _roll = ( ((Rotorpart*)getRotorpart(0))->getrealAlpha()
285 -((Rotorpart*)getRotorpart(2*(_number_of_parts>>2)))->getrealAlpha()
287 *f=(_balance1>-1)?_roll *180/pi:0;
292 sprintf(text,"/rotors/%s/yaw-deg", _name);
293 _yaw=( ((Rotorpart*)getRotorpart(1*(_number_of_parts>>2)))->getrealAlpha()
294 -((Rotorpart*)getRotorpart(3*(_number_of_parts>>2)))->getrealAlpha()
296 *f=(_balance1>-1)?_yaw*180/pi:0;
301 sprintf(text,"/rotors/%s/rpm", _name);
302 *f=(_balance1>-1)?_omega/2/pi*60:0;
307 sprintf(text,"/rotors/%s/tilt/pitch-deg",_name);
308 *f=_tilt_pitch*180/pi;
312 sprintf(text,"/rotors/%s/tilt/roll-deg",_name);
313 *f=_tilt_roll*180/pi;
317 sprintf(text,"/rotors/%s/tilt/yaw-deg",_name);
322 sprintf(text,"/rotors/%s/balance", _name);
327 sprintf(text,"/rotors/%s/stall",_name);
328 *f=getOverallStall();
332 sprintf(text,"/rotors/%s/torque",_name);
338 if (b>=_number_of_blades)
343 sprintf(text,"/rotors/%s/blade[%i]/%s",
345 w==0?"position-deg":(w==1?"flap-deg":"incidence-deg"));
346 *f=((Rotorpart*)getRotorpart(0))->getPhi()*180/pi
347 +360*b/_number_of_blades*(_ccw?1:-1);
350 if (_balance1<=-1) *f=0;
357 rk=Math::clamp(rk,0,1);//Delete this
359 if(w==2) {k+=2;l+=2;}
361 if(w==1) {k+=1;l+=1;}
364 if (w==1) *f=rk*((Rotorpart*) getRotorpart(k*(_number_of_parts>>2)))->getrealAlpha()*180/pi
365 +rl*((Rotorpart*) getRotorpart(l*(_number_of_parts>>2)))->getrealAlpha()*180/pi;
366 else if(w==2) *f=rk*((Rotorpart*)getRotorpart(k*(_number_of_parts>>2)))->getIncidence()*180/pi
367 +rl*((Rotorpart*)getRotorpart(l*(_number_of_parts>>2)))->getIncidence()*180/pi;
372 void Rotorgear::setEngineOn(int value)
377 void Rotorgear::setRotorEngineMaxRelTorque(float lval)
379 _max_rel_torque=lval;
382 void Rotorgear::setRotorRelTarget(float lval)
384 _target_rel_rpm=lval;
387 void Rotor::setAlpha0(float f)
392 void Rotor::setAlphamin(float f)
397 void Rotor::setAlphamax(float f)
402 void Rotor::setAlpha0factor(float f)
407 int Rotor::numRotorparts()
409 return _rotorparts.size();
412 Rotorpart* Rotor::getRotorpart(int n)
414 return ((Rotorpart*)_rotorparts.get(n));
417 int Rotorgear::getEngineon()
422 float Rotor::getTorqueOfInertia()
424 return _torque_of_inertia;
427 void Rotor::setTorque(float f)
432 void Rotor::addTorque(float f)
437 void Rotor::strncpy(char *dest,const char *src,int maxlen)
440 while(src[n]&&n<(maxlen-1))
448 void Rotor::setNormal(float* normal)
451 float invsum,sqrsum=0;
452 for(i=0; i<3; i++) { sqrsum+=normal[i]*normal[i];}
454 invsum=1/Math::sqrt(sqrsum);
459 _normal_with_yaw_roll[i]=_normal[i] = normal[i]*invsum;
463 void Rotor::setForward(float* forward)
466 float invsum,sqrsum=0;
467 for(i=0; i<3; i++) { sqrsum+=forward[i]*forward[i];}
469 invsum=1/Math::sqrt(sqrsum);
472 for(i=0; i<3; i++) { _forward[i] = forward[i]*invsum; }
475 void Rotor::setForceAtPitchA(float force)
477 _force_at_pitch_a=force;
480 void Rotor::setPowerAtPitch0(float value)
482 _power_at_pitch_0=value;
485 void Rotor::setPowerAtPitchB(float value)
487 _power_at_pitch_b=value;
490 void Rotor::setPitchA(float value)
492 _pitch_a=value/180*pi;
495 void Rotor::setPitchB(float value)
497 _pitch_b=value/180*pi;
500 void Rotor::setBase(float* base)
503 for(i=0; i<3; i++) _base[i] = base[i];
506 void Rotor::setMaxCyclicail(float value)
508 _maxcyclicail=value/180*pi;
511 void Rotor::setMaxCyclicele(float value)
513 _maxcyclicele=value/180*pi;
516 void Rotor::setMinCyclicail(float value)
518 _mincyclicail=value/180*pi;
521 void Rotor::setMinCyclicele(float value)
523 _mincyclicele=value/180*pi;
526 void Rotor::setMinCollective(float value)
528 _min_pitch=value/180*pi;
531 void Rotor::setMinTiltYaw(float value)
533 _min_tilt_yaw=value/180*pi;
536 void Rotor::setMinTiltPitch(float value)
538 _min_tilt_pitch=value/180*pi;
541 void Rotor::setMinTiltRoll(float value)
543 _min_tilt_roll=value/180*pi;
546 void Rotor::setMaxTiltYaw(float value)
548 _max_tilt_yaw=value/180*pi;
551 void Rotor::setMaxTiltPitch(float value)
553 _max_tilt_pitch=value/180*pi;
556 void Rotor::setMaxTiltRoll(float value)
558 _max_tilt_roll=value/180*pi;
561 void Rotor::setTiltCenterX(float value)
563 _tilt_center[0] = value;
566 void Rotor::setTiltCenterY(float value)
568 _tilt_center[1] = value;
571 void Rotor::setTiltCenterZ(float value)
573 _tilt_center[2] = value;
576 void Rotor::setMaxCollective(float value)
578 _max_pitch=value/180*pi;
581 void Rotor::setDiameter(float value)
586 void Rotor::setWeightPerBlade(float value)
588 _weight_per_blade=value;
591 void Rotor::setNumberOfBlades(float value)
593 _number_of_blades=int(value+.5);
596 void Rotor::setRelBladeCenter(float value)
598 _rel_blade_center=value;
601 void Rotor::setDynamic(float value)
606 void Rotor::setDelta3(float value)
611 void Rotor::setDelta(float value)
616 void Rotor::setTranslift(float value)
621 void Rotor::setSharedFlapHinge(bool s)
623 _shared_flap_hinge=s;
626 void Rotor::setBalance(float b)
631 void Rotor::setC2(float value)
636 void Rotor::setStepspersecond(float steps)
638 _stepspersecond=steps;
641 void Rotor::setRPM(float value)
646 void Rotor::setPhiNull(float value)
651 void Rotor::setRelLenHinge(float value)
653 _rel_len_hinge=value;
656 void Rotor::setDownwashFactor(float value)
658 _downwash_factor=value;
661 void Rotor::setAlphaoutput(int i, const char *text)
663 strncpy(_alphaoutput[i],text,255);
666 void Rotor::setName(const char *text)
668 strncpy(_name,text,256);//256: some space needed for settings
671 void Rotor::setCcw(int ccw)
676 void Rotor::setNotorque(int value)
681 void Rotor::setRelLenTeeterHinge(float f)
683 _rellenteeterhinge=f;
686 void Rotor::setTeeterdamp(float f)
691 void Rotor::setMaxteeterdamp(float f)
696 void Rotor::setGlobalGround(double *global_ground, float* global_vel)
699 for(i=0; i<4; i++) _global_ground[i] = global_ground[i];
702 void Rotor::setParameter(char *parametername, float value)
704 #define p(a,b) if (strcmp(parametername,#a)==0) _##a = (value * (b)); else
706 p(translift_maxfactor,1)
707 p(ground_effect_constant,1)
708 p(vortex_state_lift_factor,1)
716 p(number_of_segments,1)
718 p(rel_len_where_incidence_is_measured,1)
721 p(airfoil_incidence_no_lift,pi/180.)
722 p(rel_len_blade_start,1)
723 p(incidence_stall_zero_speed,pi/180.)
724 p(incidence_stall_half_sonic_speed,pi/180.)
725 p(lift_factor_stall,1)
726 p(stall_change_over,pi/180.)
727 p(drag_factor_stall,1)
728 p(airfoil_lift_coefficient,1)
729 p(airfoil_drag_coefficient0,1)
730 p(airfoil_drag_coefficient1,1)
732 p(rotor_correction_factor,1)
733 SG_LOG(SG_INPUT, SG_ALERT,
734 "internal error in parameter set up for rotor: '" <<
735 parametername <<"'" << endl);
739 float Rotor::getLiftFactor()
744 void Rotorgear::setRotorBrake(float lval)
746 lval = Math::clamp(lval, 0, 1);
750 void Rotor::setTiltYaw(float lval)
752 lval = Math::clamp(lval, -1, 1);
753 _tilt_yaw = _min_tilt_yaw+(lval+1)/2*(_max_tilt_yaw-_min_tilt_yaw);
754 _directions_and_postions_dirty = true;
757 void Rotor::setTiltPitch(float lval)
759 lval = Math::clamp(lval, -1, 1);
760 _tilt_pitch = _min_tilt_pitch+(lval+1)/2*(_max_tilt_pitch-_min_tilt_pitch);
761 _directions_and_postions_dirty = true;
764 void Rotor::setTiltRoll(float lval)
766 lval = Math::clamp(lval, -1, 1);
767 _tilt_roll = _min_tilt_roll+(lval+1)/2*(_max_tilt_roll-_min_tilt_roll);
768 _directions_and_postions_dirty = true;
771 void Rotor::setCollective(float lval)
773 lval = Math::clamp(lval, -1, 1);
775 _collective=_min_pitch+(lval+1)/2*(_max_pitch-_min_pitch);
776 for(i=0; i<_rotorparts.size(); i++) {
777 ((Rotorpart*)_rotorparts.get(i))->setCollective(_collective);
781 void Rotor::setCyclicele(float lval,float rval)
783 lval = Math::clamp(lval, -1, 1);
784 _cyclicele=_mincyclicele+(lval+1)/2*(_maxcyclicele-_mincyclicele);
787 void Rotor::setCyclicail(float lval,float rval)
789 lval = Math::clamp(lval, -1, 1);
791 _cyclicail=-(_mincyclicail+(lval+1)/2*(_maxcyclicail-_mincyclicail));
794 void Rotor::getPosition(float* out)
797 for(i=0; i<3; i++) out[i] = _base[i];
800 void Rotor::calcLiftFactor(float* v, float rho, State *s)
802 //calculates _lift_factor, which is a foactor for the lift of the rotor
804 //- ground effect (_f_ge)
805 //- vortex state (_f_vs)
806 //- translational lift (_f_tl)
811 // Calculate ground effect
812 _f_ge=1+_diameter/_ground_effect_altitude*_ground_effect_constant;
814 // Now calculate translational lift
815 float v_vert = Math::dot3(v,_normal);
817 Math::cross3(v,_normal,help);
818 float v_horiz = Math::mag3(help);
819 _f_tl = ((1-Math::pow(2.7183,-v_horiz/_translift_ve))
820 *(_translift_maxfactor-1)+1)/_translift_maxfactor;
822 _lift_factor = _f_ge*_f_tl*_f_vs;
824 //store the gravity direction
825 Glue::geodUp(s->pos, _grav_direction);
828 void Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s)
830 _ground_effect_altitude=findGroundEffectAltitude(ground_cb,s,
831 _groundeffectpos[0],_groundeffectpos[1],
832 _groundeffectpos[2],_groundeffectpos[3]);
833 testForRotorGroundContact(ground_cb,s);
836 void Rotor::testForRotorGroundContact(Ground * ground_cb,State *s)
839 for (i=0;i<_num_ground_contact_pos;i++)
842 s->posLocalToGlobal(_ground_contact_pos[i], pt);
844 // Ask for the ground plane in the global coordinate system
845 double global_ground[4];
847 ground_cb->getGroundPlane(pt, global_ground, global_vel);
848 // find h, the distance to the ground
849 // The ground plane transformed to the local frame.
851 s->planeGlobalToLocal(global_ground, ground);
853 h = ground[3] - Math::dot3(_ground_contact_pos[i], ground);
854 // Now h is the distance from _ground_contact_pos[i] to ground
857 _balance1 -= (-h)/_diameter/_num_ground_contact_pos;
858 _balance1 = (_balance1<-1)?-1:_balance1;
862 float Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s,
863 float *pos0,float *pos1,float *pos2,float *pos3,
864 int iteration,float a0,float a1,float a2,float a3)
878 Math::add3(p[0],p[2],p[4]);
879 Math::mul3(0.5,p[4],p[4]);//the center
881 float mina=100*_diameter;
883 for (int i=0;i<5;i++)
885 if (a[i]==-1)//in the first iteration,(iteration==0) no height is
886 //passed to this function, these missing values are
890 s->posLocalToGlobal(p[i], pt);
892 // Ask for the ground plane in the global coordinate system
893 double global_ground[4];
895 ground_cb->getGroundPlane(pt, global_ground, global_vel);
896 // find h, the distance to the ground
897 // The ground plane transformed to the local frame.
899 s->planeGlobalToLocal(global_ground, ground);
901 a[i] = ground[3] - Math::dot3(p[i], ground);
902 // Now a[i] is the distance from p[i] to ground
908 if (mina>=10*_diameter)
909 return mina; //the ground effect will be zero
911 //check if further recursion is neccessary
912 //if the height does not differ more than 20%, than
913 //we can return then mean height, if not split
914 //zhe square to four parts and calcualte the height
916 //suma * 0.2 is the mean
917 //0.15 is the maximum allowed difference from the mean
918 //to the height at the center
920 ||(Math::abs(suma*0.2-a[4])<(0.15*0.2*suma*(1<<iteration))))
923 float pc[4][3],ac[4]; //pc[i]=center of pos[i] and pos[(i+1)&3]
924 for (int i=0;i<4;i++)
926 Math::add3(p[i],p[(i+1)&3],pc[i]);
927 Math::mul3(0.5,pc[i],pc[i]);
929 s->posLocalToGlobal(pc[i], pt);
931 // Ask for the ground plane in the global coordinate system
932 double global_ground[4];
934 ground_cb->getGroundPlane(pt, global_ground, global_vel);
935 // find h, the distance to the ground
936 // The ground plane transformed to the local frame.
938 s->planeGlobalToLocal(global_ground, ground);
940 ac[i] = ground[3] - Math::dot3(p[i], ground);
941 // Now ac[i] is the distance from pc[i] to ground
944 (findGroundEffectAltitude(ground_cb,s,p[0],pc[1],p[4],pc[3],
945 iteration+1,a[0],ac[0],a[4],ac[3])
946 +findGroundEffectAltitude(ground_cb,s,p[1],pc[0],p[4],pc[1],
947 iteration+1,a[1],ac[0],a[4],ac[1])
948 +findGroundEffectAltitude(ground_cb,s,p[2],pc[1],p[4],pc[2],
949 iteration+1,a[2],ac[1],a[4],ac[2])
950 +findGroundEffectAltitude(ground_cb,s,p[3],pc[2],p[4],pc[3],
951 iteration+1,a[3],ac[2],a[4],ac[3])
955 void Rotor::getDownWash(float *pos, float *v_heli, float *downwash)
957 float pos2rotor[3],tmp[3];
958 Math::sub3(_base,pos,pos2rotor);
959 float dist=Math::dot3(pos2rotor,_normal_with_yaw_roll);
960 //calculate incidence at 0.7r;
961 float inc = _collective+_twist *0.7
962 - _twist*_rel_len_where_incidence_is_measured;
965 if (dist<0) // we are not in the downwash region
967 downwash[0]=downwash[1]=downwash[2]=0.;
971 //calculate the mean downwash speed directly beneath the rotor disk
972 float v1bar = Math::sin(inc) *_omega * 0.35 * _diameter * 0.8;
973 //0.35 * d = 0.7 *r, a good position to calcualte the mean downwashd
974 //0.8 the slip of the rotor.
976 //calculate the time the wind needed from thr rotor to here
977 if (v1bar< 1) v1bar = 1;
978 float time=dist/v1bar;
980 //calculate the pos2rotor, where the rotor was, "time" ago
981 Math::mul3(time,v_heli,tmp);
982 Math::sub3(pos2rotor,tmp,pos2rotor);
984 //and again calculate dist
985 dist=Math::dot3(pos2rotor,_normal_with_yaw_roll);
986 //missing the normal is offen not pointing to the normal of the rotor
987 //disk. Rotate the normal by yaw and tilt angle
991 if (dist<0) // we are not in the downwash region
993 downwash[0]=downwash[1]=downwash[2]=0.;
996 //of course this could be done in a runge kutta integrator, but it's such
997 //a approximation that I beleave, it would'nt be more realistic
999 //calculate the dist to the rotor-axis
1000 Math::cross3(pos2rotor,_normal_with_yaw_roll,tmp);
1001 float r= Math::mag3(tmp);
1002 //calculate incidence at r;
1003 float rel_r = r *2 /_diameter;
1004 float inc_r = _collective+_twist * r /_diameter * 2
1005 - _twist*_rel_len_where_incidence_is_measured;
1007 //calculate the downwash speed directly beneath the rotor disk
1010 v1 = Math::sin(inc_r) *_omega * r * 0.8;
1012 //calcualte the downwash speed in a distance "dist" to the rotor disc,
1013 //for large dist. The speed is assumed do follow a gausian distribution
1014 //with sigma increasing with dist^2:
1015 //sigma is assumed to be half of the rotor diameter directly beneath the
1016 //disc and is assumed to the rotor diameter at dist = (diameter * sqrt(2))
1018 float sigma=_diameter/2 + dist * dist / _diameter /4.;
1019 float v2 = v1bar*_diameter/ (Math::sqrt(2 * pi) * sigma)
1020 * Math::pow(2.7183,-.5*r*r/(sigma*sigma))*_diameter/2/sigma;
1022 //calculate the weight of the two downwash velocities.
1023 //Directly beneath the disc it is v1, far away it is v2
1024 float g = Math::pow(2.7183,-2*dist/_diameter);
1025 //at dist = rotor radius it is assumed to be 1/e * v1 + (1-1/e)* v2
1027 float v = g * v1 + (1-g) * v2;
1028 Math::mul3(-v*_downwash_factor,_normal_with_yaw_roll,downwash);
1029 //the downwash is calculated in the opposite direction of the normal
1032 void Rotor::euler2orient(float roll, float pitch, float hdg, float* out)
1034 // the Glue::euler2orient, inverts y<z due to different bases
1035 // therefore the negation of all "y" and "z" coeffizients
1036 Glue::euler2orient(roll,pitch,hdg,out);
1037 for (int i=3;i<9;i++) out[i]*=-1.0;
1041 void Rotor::updateDirectionsAndPositions(float *rot)
1043 if (!_directions_and_postions_dirty)
1045 rot[0]=rot[1]=rot[2]=0;
1048 rot[0]=_old_tilt_roll-_tilt_roll;
1049 rot[1]=_old_tilt_pitch-_tilt_pitch;
1050 rot[2]=_old_tilt_yaw-_tilt_yaw;
1051 _old_tilt_roll=_tilt_roll;
1052 _old_tilt_pitch=_tilt_pitch;
1053 _old_tilt_yaw=_tilt_yaw;
1055 euler2orient(_tilt_roll, _tilt_pitch, _tilt_yaw, orient);
1059 Math::sub3(_base,_tilt_center,base);
1060 Math::vmul33(orient, base, base);
1061 Math::add3(base,_tilt_center,base);
1062 Math::vmul33(orient, _forward, forward);
1063 Math::vmul33(orient, _normal, normal);
1065 #define _forward forward
1066 #define _normal normal
1067 float directions[5][3];
1068 //pointing forward, right, ... the 5th is ony for calculation
1069 directions[0][0]=_forward[0];
1070 directions[0][1]=_forward[1];
1071 directions[0][2]=_forward[2];
1076 Math::cross3(directions[i-1],_normal,directions[i]);
1078 Math::cross3(_normal,directions[i-1],directions[i]);
1080 Math::set3(directions[4],directions[0]);
1081 // now directions[0] is perpendicular to the _normal.and has a length
1082 // of 1. if _forward is already normalized and perpendicular to the
1083 // normal, directions[0] will be the same
1084 //_num_ground_contact_pos=(_number_of_parts<16)?_number_of_parts:16;
1085 for (i=0;i<_num_ground_contact_pos;i++)
1088 float s = Math::sin(pi*2*i/_num_ground_contact_pos);
1089 float c = Math::cos(pi*2*i/_num_ground_contact_pos);
1090 Math::mul3(c*_diameter*0.5,directions[0],_ground_contact_pos[i]);
1091 Math::mul3(s*_diameter*0.5,directions[1],help);
1092 Math::add3(help,_ground_contact_pos[i],_ground_contact_pos[i]);
1093 Math::add3(_base,_ground_contact_pos[i],_ground_contact_pos[i]);
1097 Math::mul3(_diameter*0.7,directions[i],_groundeffectpos[i]);
1098 Math::add3(_base,_groundeffectpos[i],_groundeffectpos[i]);
1100 for (i=0;i<_number_of_parts;i++)
1102 Rotorpart* rp = getRotorpart(i);
1103 float lpos[3],lforceattac[3],lspeed[3],dirzentforce[3];
1104 float s = Math::sin(2*pi*i/_number_of_parts);
1105 float c = Math::cos(2*pi*i/_number_of_parts);
1106 float sp = Math::sin(float(2*pi*i/_number_of_parts-pi/2.+_phi));
1107 float cp = Math::cos(float(2*pi*i/_number_of_parts-pi/2.+_phi));
1108 float direction[3],nextdirection[3],help[3],direction90deg[3];
1109 float rotorpartmass = _weight_per_blade*_number_of_blades/_number_of_parts*.453;
1110 float speed=_rotor_rpm/60*_diameter*_rel_blade_center*pi;
1111 float lentocenter=_diameter*_rel_blade_center*0.5;
1112 float lentoforceattac=_diameter*_rel_len_hinge*0.5;
1113 float zentforce=rotorpartmass*speed*speed/lentocenter;
1115 Math::mul3(c ,directions[0],help);
1116 Math::mul3(s ,directions[1],direction);
1117 Math::add3(help,direction,direction);
1119 Math::mul3(c ,directions[1],help);
1120 Math::mul3(s ,directions[2],direction90deg);
1121 Math::add3(help,direction90deg,direction90deg);
1123 Math::mul3(cp ,directions[1],help);
1124 Math::mul3(sp ,directions[2],nextdirection);
1125 Math::add3(help,nextdirection,nextdirection);
1127 Math::mul3(lentocenter,direction,lpos);
1128 Math::add3(lpos,_base,lpos);
1129 Math::mul3(lentoforceattac,nextdirection,lforceattac);
1130 //nextdirection: +90deg (gyro)!!!
1132 Math::add3(lforceattac,_base,lforceattac);
1133 Math::mul3(speed,direction90deg,lspeed);
1134 Math::mul3(1,nextdirection,dirzentforce);
1135 rp->setPosition(lpos);
1136 rp->setNormal(_normal);
1137 rp->setZentipetalForce(zentforce);
1138 rp->setPositionForceAttac(lforceattac);
1139 rp->setSpeed(lspeed);
1140 rp->setDirectionofZentipetalforce(dirzentforce);
1141 rp->setDirectionofRotorPart(direction);
1146 _directions_and_postions_dirty=false;
1149 void Rotor::compile()
1151 // Have we already been compiled?
1152 if(_rotorparts.size() != 0) return;
1154 //rotor is divided into _number_of_parts parts
1155 //each part is calcualted at _number_of_segments points
1158 //and make it a factor of 4
1159 _number_of_parts=(int(Math::clamp(_number_of_parts,4,256))>>2)<<2;
1161 _dynamic=_dynamic*(1/ //inverse of the time
1162 ( (60/_rotor_rpm)/4 //for rotating 90 deg
1163 +(60/_rotor_rpm)/(2*_number_of_blades) //+ meantime a rotorblade
1164 //will pass a given point
1166 //normalize the directions
1167 Math::unit3(_forward,_forward);
1168 Math::unit3(_normal,_normal);
1169 _num_ground_contact_pos=(_number_of_parts<16)?_number_of_parts:16;
1170 float rotorpartmass = _weight_per_blade*_number_of_blades/_number_of_parts*.453;
1171 //was pounds -> now kg
1173 _torque_of_inertia = 1/12. * ( _number_of_parts * rotorpartmass) * _diameter
1174 * _diameter * _rel_blade_center * _rel_blade_center /(0.5*0.5);
1175 float speed=_rotor_rpm/60*_diameter*_rel_blade_center*pi;
1176 float lentocenter=_diameter*_rel_blade_center*0.5;
1177 float lentoforceattac=_diameter*_rel_len_hinge*0.5;
1178 float zentforce=rotorpartmass*speed*speed/lentocenter;
1179 float pitchaforce=_force_at_pitch_a/_number_of_parts*.453*9.81;
1180 // was pounds of force, now N, devided by _number_of_parts
1181 //(so its now per rotorpart)
1183 float torque0=0,torquemax=0,torqueb=0;
1184 float omega=_rotor_rpm/60*2*pi;
1186 float omega0=omega*Math::sqrt(1/(1-_rel_len_hinge));
1187 float delta_theoretical=pitchaforce/(_pitch_a*omega*lentocenter*2*rotorpartmass);
1188 _delta*=delta_theoretical;
1190 float relamp=(omega*omega/(2*_delta*Math::sqrt(sqr(omega0*omega0-omega*omega)
1191 +4*_delta*_delta*omega*omega)))*_cyclic_factor;
1192 float relamp_theoretical=(omega*omega/(2*delta_theoretical*Math::sqrt(sqr(omega0*omega0-omega*omega)
1193 +4*delta_theoretical*delta_theoretical*omega*omega)))*_cyclic_factor;
1194 _phi=Math::acos(_rel_len_hinge);
1195 _phi-=Math::atan(_delta3);
1198 torque0=_power_at_pitch_0/_number_of_parts*1000/omega;
1199 // f*r=p/w ; p=f*s/t; r=s/t/w ; r*w*t = s
1200 torqueb=_power_at_pitch_b/_number_of_parts*1000/omega;
1201 torquemax=_power_at_pitch_b/_number_of_parts*1000/omega/_pitch_b*_max_pitch;
1211 Rotorpart* rps[256];
1213 for (i=0;i<_number_of_parts;i++)
1215 Rotorpart* rp=rps[i]=newRotorpart(zentforce,pitchaforce,_delta3,rotorpartmass,
1216 _translift,_rel_len_hinge,lentocenter);
1217 int k = i*4/_number_of_parts;
1218 rp->setAlphaoutput(_alphaoutput[k&1?k:(_ccw?k^2:k)],0);
1219 rp->setAlphaoutput(_alphaoutput[4+(k&1?k:(_ccw?k^2:k))],1+(k>1));
1220 _rotorparts.add(rp);
1221 rp->setTorque(torquemax,torque0);
1222 rp->setRelamp(relamp);
1223 rp->setTorqueOfInertia(_torque_of_inertia/_number_of_parts);
1224 rp->setDirection(2*pi*i/_number_of_parts);
1226 for (i=0;i<_number_of_parts;i++)
1228 rps[i]->setlastnextrp(rps[(i-1+_number_of_parts)%_number_of_parts],
1229 rps[(i+1)%_number_of_parts],
1230 rps[(i+_number_of_parts/2)%_number_of_parts],
1231 rps[(i-_number_of_parts/4+_number_of_parts)%_number_of_parts],
1232 rps[(i+_number_of_parts/4)%_number_of_parts]);
1235 updateDirectionsAndPositions(drot);
1236 for (i=0;i<_number_of_parts;i++)
1238 rps[i]->setCompiled();
1240 float lift[4],torque[4], v_wind[3];
1241 v_wind[0]=v_wind[1]=v_wind[2]=0;
1242 rps[0]->setOmega(_omegan);
1244 if (_airfoil_lift_coefficient==0)
1246 //calculate the lift and drag coefficients now
1250 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_a,0,0,
1251 &(torque[0]),&(lift[0])); //max_pitch a
1252 _liftcoef = pitchaforce/lift[0];
1255 rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,&(torque[0]),&(lift[0]));
1260 rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,&(torque[1]),&(lift[1]));
1265 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,&(torque[2]),&(lift[2]));
1270 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,&(torque[3]),&(lift[3]));
1275 _dragcoef1=torque0/torque[1];
1276 _dragcoef0=(torqueb-_dragcoef1*torque[3])/torque[2];
1280 _dragcoef1=(torque0/torque[0]-torqueb/torque[2])
1281 /(torque[1]/torque[0]-torque[3]/torque[2]);
1282 _dragcoef0=(torqueb-_dragcoef1*torque[3])/torque[2];
1287 _liftcoef=_airfoil_lift_coefficient/_number_of_parts*_number_of_blades;
1288 _dragcoef0=_airfoil_drag_coefficient0/_number_of_parts*_number_of_blades*_c2;
1289 _dragcoef1=_airfoil_drag_coefficient1/_number_of_parts*_number_of_blades*_c2;
1293 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_a,0,0,
1294 &(torque[0]),&(lift[0])); //pitch a
1295 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,
1296 &(torque[1]),&(lift[1])); //pitch b
1297 rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,
1298 &(torque[3]),&(lift[3])); //pitch 0
1299 SG_LOG(SG_GENERAL, SG_INFO,
1300 "Rotor: coefficients for airfoil:" << endl << setprecision(6)
1301 << " drag0: " << _dragcoef0*_number_of_parts/_number_of_blades/_c2
1302 << " drag1: " << _dragcoef1*_number_of_parts/_number_of_blades/_c2
1303 << " lift: " << _liftcoef*_number_of_parts/_number_of_blades
1305 << "at 10 deg:" << endl
1306 << "drag: " << (Math::sin(10./180*pi)*_dragcoef1+_dragcoef0)
1307 *_number_of_parts/_number_of_blades/_c2
1308 << " lift: " << Math::sin(10./180*pi)*_liftcoef*_number_of_parts/_number_of_blades
1310 << "Some results (Pitch [degree], Power [kW], Lift [N])" << endl
1311 << 0.0f << "deg " << Math::abs(torque[3]*_number_of_parts*_omegan/1000) << "kW "
1312 << lift[3]*_number_of_parts << endl
1313 << _pitch_a*180/pi << "deg " << Math::abs(torque[0]*_number_of_parts*_omegan/1000)
1314 << "kW " << lift[0]*_number_of_parts << endl
1315 << _pitch_b*180/pi << "deg " << Math::abs(torque[1]*_number_of_parts*_omegan/1000)
1316 << "kW " << lift[1]*_number_of_parts << endl << endl );
1318 //first calculation of relamp is wrong
1319 //it used pitchaforce, but this was unknown and
1320 //on the default value
1321 _delta*=lift[0]/pitchaforce;
1322 relamp=(omega*omega/(2*_delta*Math::sqrt(sqr(omega0*omega0-omega*omega)
1323 +4*_delta*_delta*omega*omega)))*_cyclic_factor;
1324 for (i=0;i<_number_of_parts;i++)
1326 rps[i]->setRelamp(relamp);
1328 rps[0]->setOmega(0);
1335 //tie the properties
1336 /* After reset these values are totally wrong. I have to find out why
1337 SGPropertyNode * node = fgGetNode("/rotors", true)->getNode(_name,true);
1338 node->tie("balance_ext",SGRawValuePointer<float>(&_balance2),false);
1339 node->tie("balance_int",SGRawValuePointer<float>(&_balance1));
1343 std::ostream & operator<<(std::ostream & out, Rotor& r)
1345 #define i(x) << #x << ":" << r.x << endl
1346 #define iv(x) << #x << ":" << r.x[0] << ";" << r.x[1] << ";" <<r.x[2] << ";" << endl
1347 out << "Writing Info on Rotor "
1350 i(_omega) i(_omegan) i(_omegarel) i(_ddt_omega) i(_omegarelneu)
1353 i( _airfoil_incidence_no_lift)
1355 i( _airfoil_lift_coefficient)
1356 i( _airfoil_drag_coefficient0)
1357 i( _airfoil_drag_coefficient1)
1359 i( _number_of_segments)
1360 i( _number_of_parts)
1362 iv( _groundeffectpos[0])iv( _groundeffectpos[1])iv( _groundeffectpos[2])iv( _groundeffectpos[3])
1363 i( _ground_effect_altitude)
1365 iv( _normal_with_yaw_roll)
1368 i( _number_of_blades)
1369 i( _weight_per_blade)
1370 i( _rel_blade_center)
1373 i( _force_at_pitch_a)
1375 i( _power_at_pitch_0)
1376 i( _power_at_pitch_b)
1393 i( _alphamin) i(_alphamax) i(_alpha0) i(_alpha0factor)
1394 i( _teeterdamp) i(_maxteeterdamp)
1395 i( _rellenteeterhinge)
1397 i( _translift_maxfactor)
1398 i( _ground_effect_constant)
1399 i( _vortex_state_lift_factor)
1400 i( _vortex_state_c1)
1401 i( _vortex_state_c2)
1402 i( _vortex_state_c3)
1403 i( _vortex_state_e1)
1404 i( _vortex_state_e2)
1405 i( _vortex_state_e3)
1406 i( _lift_factor) i(_f_ge) i(_f_vs) i(_f_tl)
1411 i( _twist) //outer incidence = inner inner incidence + _twist
1412 i( _rel_len_where_incidence_is_measured)
1413 i( _torque_of_inertia)
1414 i( _rel_len_blade_start)
1415 i( _incidence_stall_zero_speed)
1416 i( _incidence_stall_half_sonic_speed)
1417 i( _lift_factor_stall)
1418 i( _stall_change_over)
1419 i( _drag_factor_stall)
1426 i( _cyclic_factor) <<endl;
1428 for(j=0; j<r._rotorparts.size(); j++) {
1429 out << *((Rotorpart*)r._rotorparts.get(j));
1436 void Rotor:: writeInfo()
1439 std::ostringstream buffer;
1441 FILE*f=fopen("c:\\fgmsvc\\bat\\log.txt","at");
1442 if (!f) f=fopen("c:\\fgmsvc\\bat\\log.txt","wt");
1445 fprintf(f,"%s",(const char *)buffer.str().c_str());
1450 Rotorpart* Rotor::newRotorpart(float zentforce,float maxpitchforce,
1451 float delta3,float mass,float translift,float rellenhinge,float len)
1453 Rotorpart *r = new Rotorpart();
1454 r->setDelta3(delta3);
1455 r->setDynamic(_dynamic);
1456 r->setTranslift(_translift);
1459 r->setRelLenHinge(rellenhinge);
1460 r->setSharedFlapHinge(_shared_flap_hinge);
1461 r->setOmegaN(_omegan);
1462 r->setPhi(_phi_null);
1463 r->setAlpha0(_alpha0);
1464 r->setAlphamin(_alphamin);
1465 r->setAlphamax(_alphamax);
1466 r->setAlpha0factor(_alpha0factor);
1468 r->setDiameter(_diameter);
1470 #define p(a) r->setParameter(#a,_##a);
1472 p(number_of_segments)
1473 p(rel_len_where_incidence_is_measured)
1474 p(rel_len_blade_start)
1475 p(rotor_correction_factor)
1480 void Rotor::interp(float* v1, float* v2, float frac, float* out)
1482 out[0] = v1[0] + frac*(v2[0]-v1[0]);
1483 out[1] = v1[1] + frac*(v2[1]-v1[1]);
1484 out[2] = v1[2] + frac*(v2[2]-v1[2]);
1487 void Rotorgear::initRotorIteration(float *lrot,float dt)
1491 if (!_rotors.size()) return;
1492 Rotor* r0 = (Rotor*)_rotors.get(0);
1493 omegarel= r0->getOmegaRelNeu();
1494 for(i=0; i<_rotors.size(); i++) {
1495 Rotor* r = (Rotor*)_rotors.get(i);
1496 r->inititeration(dt,omegarel,0,lrot);
1500 void Rotorgear::calcForces(float* torqueOut)
1503 torqueOut[0]=torqueOut[1]=torqueOut[2]=0;
1504 // check,<if the engine can handle the torque of the rotors.
1505 // If not reduce the torque to the fueselage and change rotational
1506 // speed of the rotors instead
1509 float omegarel,omegan;
1510 Rotor* r0 = (Rotor*)_rotors.get(0);
1511 omegarel= r0->getOmegaRel();
1513 float total_torque_of_inertia=0;
1514 float total_torque=0;
1515 for(i=0; i<_rotors.size(); i++) {
1516 Rotor* r = (Rotor*)_rotors.get(i);
1517 omegan=r->getOmegan();
1518 total_torque_of_inertia+=r->getTorqueOfInertia()*omegan*omegan;
1519 //FIXME: this is constant, so this can be done in compile
1521 total_torque+=r->getTorque()*omegan;
1523 float max_torque_of_engine=0;
1524 SGPropertyNode * node=fgGetNode("/rotors/gear", true);
1527 max_torque_of_engine=_max_power_engine*_max_rel_torque;
1528 float df=_target_rel_rpm-omegarel;
1529 df/=_engine_prop_factor;
1530 df = Math::clamp(df, 0, 1);
1531 max_torque_of_engine = df * _max_power_engine*_max_rel_torque;
1535 float rel_torque_engine=1;
1536 if (total_torque<=0)
1537 rel_torque_engine=0;
1539 if (max_torque_of_engine>0)
1540 rel_torque_engine=1/max_torque_of_engine*total_torque;
1542 rel_torque_engine=0;
1544 //add the rotor brake and the gear fritcion
1546 if (r0->_rotorparts.size()) dt=((Rotorpart*)r0->_rotorparts.get(0))->getDt();
1548 float rotor_brake_torque;
1549 rotor_brake_torque=_rotorbrake*_max_power_rotor_brake+_rotorgear_friction;
1550 //clamp it to the value you need to stop the rotor
1551 //to avod accelerate the rotor to neagtive rpm:
1552 rotor_brake_torque=Math::clamp(rotor_brake_torque,0,
1553 total_torque_of_inertia/dt*omegarel);
1554 max_torque_of_engine-=rotor_brake_torque;
1556 //change the rotation of the rotors
1557 if ((max_torque_of_engine<total_torque) //decreasing rotation
1558 ||((max_torque_of_engine>total_torque)&&(omegarel<_target_rel_rpm))
1559 //increasing rotation due to engine
1560 ||(total_torque<0) ) //increasing rotation due to autorotation
1562 _ddt_omegarel=(max_torque_of_engine-total_torque)/total_torque_of_inertia;
1563 if(max_torque_of_engine>total_torque)
1565 //check if the acceleration is due to the engine. If yes,
1566 //the engine self limits the accel.
1567 float lim1=-total_torque/total_torque_of_inertia;
1568 //accel. by autorotation
1570 if (lim1<_engine_accel_limit) lim1=_engine_accel_limit;
1571 //if the accel by autorotation greater than the max. engine
1572 //accel, then this is the limit, if not: the engine is the limit
1573 if (_ddt_omegarel>lim1) _ddt_omegarel=lim1;
1575 if (_ddt_omegarel>5.5)_ddt_omegarel=5.5;
1576 //clamp it to avoid overflow. Should never be reached
1577 if (_ddt_omegarel<-5.5)_ddt_omegarel=-5.5;
1579 if (_max_power_engine<0.001) {omegarel=1;_ddt_omegarel=0;}
1580 //for debug: negative or no maxpower will result
1581 //in permanet 100% rotation
1583 omegarel+=dt*_ddt_omegarel;
1585 if (omegarel>2.5) omegarel=2.5;
1586 //clamp it to avoid overflow. Should never be reached
1587 if (omegarel<-.5) omegarel=-.5;
1589 r0->setOmegaRelNeu(omegarel);
1590 //calculate the torque, which is needed to accelerate the rotors.
1591 //Add this additional torque to the body
1592 for(j=0; j<_rotors.size(); j++) {
1593 Rotor* r = (Rotor*)_rotors.get(j);
1594 for(i=0; i<r->_rotorparts.size(); i++) {
1595 float torque_scalar=0;
1596 Rotorpart* rp = (Rotorpart*)r->_rotorparts.get(i);
1598 rp->getAccelTorque(_ddt_omegarel,torque);
1599 Math::add3(torque,torqueOut,torqueOut);
1603 _total_torque_on_engine=total_torque+_ddt_omegarel*total_torque_of_inertia;
1607 void Rotorgear::addRotor(Rotor* rotor)
1613 void Rotorgear::compile()
1616 for(int j=0; j<_rotors.size(); j++) {
1617 Rotor* r = (Rotor*)_rotors.get(j);
1622 void Rotorgear::getDownWash(float *pos, float * v_heli, float *downwash)
1625 downwash[0]=downwash[1]=downwash[2]=0;
1626 for(int i=0; i<_rotors.size(); i++) {
1627 Rotor* ro = (Rotor*)_rotors.get(i);
1628 ro->getDownWash(pos,v_heli,tmp);
1629 Math::add3(downwash,tmp,downwash); // + downwash
1633 void Rotorgear::setParameter(char *parametername, float value)
1635 #define p(a,b) if (strcmp(parametername,#a)==0) _##a = (value * (b)); else
1636 p(max_power_engine,1000)
1637 p(engine_prop_factor,1)
1638 p(yasimdragfactor,1)
1639 p(yasimliftfactor,1)
1640 p(max_power_rotor_brake,1000)
1641 p(rotorgear_friction,1000)
1642 p(engine_accel_limit,0.01)
1643 SG_LOG(SG_INPUT, SG_ALERT,
1644 "internal error in parameter set up for rotorgear: '"
1645 << parametername <<"'" << endl);
1648 int Rotorgear::getValueforFGSet(int j,char *text,float *f)
1652 sprintf(text,"/rotors/gear/total-torque");
1653 *f=_total_torque_on_engine;
1657 Rotorgear::Rotorgear()
1662 _max_power_rotor_brake=1;
1663 _rotorgear_friction=1;
1664 _max_power_engine=1000*450;
1665 _engine_prop_factor=0.05f;
1669 _engine_accel_limit=0.05f;
1670 _total_torque_on_engine=0;
1675 Rotorgear::~Rotorgear()
1677 for(int i=0; i<_rotors.size(); i++)
1678 delete (Rotor*)_rotors.get(i);
1681 }; // namespace yasim