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;
158 for(i=0; i<_rotorparts.size(); i++) {
159 Rotorpart* r = (Rotorpart*)_rotorparts.get(i);
162 //untie the properties
165 SGPropertyNode * node = fgGetNode("/rotors", true)->getNode(_name,true);
166 node->untie("balance-ext");
167 node->untie("balance-int");
172 void Rotor::inititeration(float dt,float omegarel,float ddt_omegarel,float *rot)
177 _omega=_omegan*_omegarel;
178 _ddt_omega=_omegan*ddt_omegarel;
181 updateDirectionsAndPositions(drot);
182 Math::add3(rot,drot,drot);
183 for(i=0; i<_rotorparts.size(); i++) {
184 float s = Math::sin(float(2*pi*i/_number_of_parts+(_phi-pi/2.)*(_ccw?1:-1)));
185 float c = Math::cos(float(2*pi*i/_number_of_parts+(_phi-pi/2.)*(_ccw?1:-1)));
186 Rotorpart* r = (Rotorpart*)_rotorparts.get(i);
188 r->setDdtOmega(_ddt_omega);
189 r->inititeration(dt,drot);
190 r->setCyclic(_cyclicail*c+_cyclicele*s);
193 //calculate the normal of the rotor disc, for calcualtion of the downwash
194 float side[3],help[3];
195 Math::cross3(_normal,_forward,side);
196 Math::mul3(Math::cos(_yaw)*Math::cos(_roll),_normal,_normal_with_yaw_roll);
198 Math::mul3(Math::sin(_yaw),_forward,help);
199 Math::add3(_normal_with_yaw_roll,help,_normal_with_yaw_roll);
201 Math::mul3(Math::sin(_roll),side,help);
202 Math::add3(_normal_with_yaw_roll,help,_normal_with_yaw_roll);
205 if ((_balance1*_balance2 < 0.97) && (_balance1>-1))
207 _balance1-=(0.97-_balance1*_balance2)*(0.97-_balance1*_balance2)*0.005;
208 if (_balance1<-1) _balance1=-1;
212 float Rotor::calcStall(float incidence,float speed)
214 float stall_incidence=_incidence_stall_zero_speed
215 +(_incidence_stall_half_sonic_speed
216 -_incidence_stall_zero_speed)*speed/(343./2);
217 //missing: Temeperature dependency of sonic speed
218 incidence = Math::abs(incidence);
219 if (incidence > (90./180.*pi))
220 incidence = pi-incidence;
221 float stall = (incidence-stall_incidence)/_stall_change_over;
222 stall = Math::clamp(stall,0,1);
224 _stall_sum+=stall*speed*speed;
225 _stall_v2sum+=speed*speed;
230 float Rotor::getLiftCoef(float incidence,float speed)
232 float stall=calcStall(incidence,speed);
233 /* the next shold look like this, but this is the inner loop of
234 the rotor simulation. For small angles (and we hav only small
235 angles) the first order approximation works well
236 float c1= Math::sin(incidence-_airfoil_incidence_no_lift)*_liftcoef;
237 for c2 we would need higher order, because in stall the angle can be large
240 if (incidence > (pi/2))
242 else if (incidence <-(pi/2))
246 float c1= (i2-_airfoil_incidence_no_lift)*_liftcoef;
249 float c2= Math::sin(2*(incidence-_airfoil_incidence_no_lift))
250 *_liftcoef*_lift_factor_stall;
251 return (1-stall)*c1 + stall *c2;
257 float Rotor::getDragCoef(float incidence,float speed)
259 float stall=calcStall(incidence,speed);
260 float c1= (Math::abs(Math::sin(incidence-_airfoil_incidence_no_lift))
261 *_dragcoef1+_dragcoef0);
262 float c2= c1*_drag_factor_stall;
263 return (1-stall)*c1 + stall *c2;
266 int Rotor::getValueforFGSet(int j,char *text,float *f)
268 if (_name[0]==0) return 0;
269 if (4>numRotorparts()) return 0; //compile first!
272 sprintf(text,"/rotors/%s/cone-deg", _name);
273 *f=(_balance1>-1)?( ((Rotorpart*)getRotorpart(0))->getrealAlpha()
274 +((Rotorpart*)getRotorpart(1*(_number_of_parts>>2)))->getrealAlpha()
275 +((Rotorpart*)getRotorpart(2*(_number_of_parts>>2)))->getrealAlpha()
276 +((Rotorpart*)getRotorpart(3*(_number_of_parts>>2)))->getrealAlpha()
282 sprintf(text,"/rotors/%s/roll-deg", _name);
283 _roll = ( ((Rotorpart*)getRotorpart(0))->getrealAlpha()
284 -((Rotorpart*)getRotorpart(2*(_number_of_parts>>2)))->getrealAlpha()
286 *f=(_balance1>-1)?_roll *180/pi:0;
291 sprintf(text,"/rotors/%s/yaw-deg", _name);
292 _yaw=( ((Rotorpart*)getRotorpart(1*(_number_of_parts>>2)))->getrealAlpha()
293 -((Rotorpart*)getRotorpart(3*(_number_of_parts>>2)))->getrealAlpha()
295 *f=(_balance1>-1)?_yaw*180/pi:0;
300 sprintf(text,"/rotors/%s/rpm", _name);
301 *f=(_balance1>-1)?_omega/2/pi*60:0;
306 sprintf(text,"/rotors/%s/tilt/pitch-deg",_name);
307 *f=_tilt_pitch*180/pi;
311 sprintf(text,"/rotors/%s/tilt/roll-deg",_name);
312 *f=_tilt_roll*180/pi;
316 sprintf(text,"/rotors/%s/tilt/yaw-deg",_name);
321 sprintf(text,"/rotors/%s/balance", _name);
326 sprintf(text,"/rotors/%s/stall",_name);
327 *f=getOverallStall();
331 sprintf(text,"/rotors/%s/torque",_name);
337 if (b>=_number_of_blades)
342 sprintf(text,"/rotors/%s/blade[%i]/%s",
344 w==0?"position-deg":(w==1?"flap-deg":"incidence-deg"));
345 *f=((Rotorpart*)getRotorpart(0))->getPhi()*180/pi
346 +360*b/_number_of_blades*(_ccw?1:-1);
349 if (_balance1<=-1) *f=0;
356 rk=Math::clamp(rk,0,1);//Delete this
358 if(w==2) {k+=2;l+=2;}
360 if(w==1) {k+=1;l+=1;}
363 if (w==1) *f=rk*((Rotorpart*) getRotorpart(k*(_number_of_parts>>2)))->getrealAlpha()*180/pi
364 +rl*((Rotorpart*) getRotorpart(l*(_number_of_parts>>2)))->getrealAlpha()*180/pi;
365 else if(w==2) *f=rk*((Rotorpart*)getRotorpart(k*(_number_of_parts>>2)))->getIncidence()*180/pi
366 +rl*((Rotorpart*)getRotorpart(l*(_number_of_parts>>2)))->getIncidence()*180/pi;
371 void Rotorgear::setEngineOn(int value)
376 void Rotorgear::setRotorEngineMaxRelTorque(float lval)
378 _max_rel_torque=lval;
381 void Rotorgear::setRotorRelTarget(float lval)
383 _target_rel_rpm=lval;
386 void Rotor::setAlpha0(float f)
391 void Rotor::setAlphamin(float f)
396 void Rotor::setAlphamax(float f)
401 void Rotor::setAlpha0factor(float f)
406 int Rotor::numRotorparts()
408 return _rotorparts.size();
411 Rotorpart* Rotor::getRotorpart(int n)
413 return ((Rotorpart*)_rotorparts.get(n));
416 int Rotorgear::getEngineon()
421 float Rotor::getTorqueOfInertia()
423 return _torque_of_inertia;
426 void Rotor::setTorque(float f)
431 void Rotor::addTorque(float f)
436 void Rotor::strncpy(char *dest,const char *src,int maxlen)
439 while(src[n]&&n<(maxlen-1))
447 void Rotor::setNormal(float* normal)
450 float invsum,sqrsum=0;
451 for(i=0; i<3; i++) { sqrsum+=normal[i]*normal[i];}
453 invsum=1/Math::sqrt(sqrsum);
458 _normal_with_yaw_roll[i]=_normal[i] = normal[i]*invsum;
462 void Rotor::setForward(float* forward)
465 float invsum,sqrsum=0;
466 for(i=0; i<3; i++) { sqrsum+=forward[i]*forward[i];}
468 invsum=1/Math::sqrt(sqrsum);
471 for(i=0; i<3; i++) { _forward[i] = forward[i]*invsum; }
474 void Rotor::setForceAtPitchA(float force)
476 _force_at_pitch_a=force;
479 void Rotor::setPowerAtPitch0(float value)
481 _power_at_pitch_0=value;
484 void Rotor::setPowerAtPitchB(float value)
486 _power_at_pitch_b=value;
489 void Rotor::setPitchA(float value)
491 _pitch_a=value/180*pi;
494 void Rotor::setPitchB(float value)
496 _pitch_b=value/180*pi;
499 void Rotor::setBase(float* base)
502 for(i=0; i<3; i++) _base[i] = base[i];
505 void Rotor::setMaxCyclicail(float value)
507 _maxcyclicail=value/180*pi;
510 void Rotor::setMaxCyclicele(float value)
512 _maxcyclicele=value/180*pi;
515 void Rotor::setMinCyclicail(float value)
517 _mincyclicail=value/180*pi;
520 void Rotor::setMinCyclicele(float value)
522 _mincyclicele=value/180*pi;
525 void Rotor::setMinCollective(float value)
527 _min_pitch=value/180*pi;
530 void Rotor::setMinTiltYaw(float value)
532 _min_tilt_yaw=value/180*pi;
535 void Rotor::setMinTiltPitch(float value)
537 _min_tilt_pitch=value/180*pi;
540 void Rotor::setMinTiltRoll(float value)
542 _min_tilt_roll=value/180*pi;
545 void Rotor::setMaxTiltYaw(float value)
547 _max_tilt_yaw=value/180*pi;
550 void Rotor::setMaxTiltPitch(float value)
552 _max_tilt_pitch=value/180*pi;
555 void Rotor::setMaxTiltRoll(float value)
557 _max_tilt_roll=value/180*pi;
560 void Rotor::setTiltCenterX(float value)
562 _tilt_center[0] = value;
565 void Rotor::setTiltCenterY(float value)
567 _tilt_center[1] = value;
570 void Rotor::setTiltCenterZ(float value)
572 _tilt_center[2] = value;
575 void Rotor::setMaxCollective(float value)
577 _max_pitch=value/180*pi;
580 void Rotor::setDiameter(float value)
585 void Rotor::setWeightPerBlade(float value)
587 _weight_per_blade=value;
590 void Rotor::setNumberOfBlades(float value)
592 _number_of_blades=int(value+.5);
595 void Rotor::setRelBladeCenter(float value)
597 _rel_blade_center=value;
600 void Rotor::setDynamic(float value)
605 void Rotor::setDelta3(float value)
610 void Rotor::setDelta(float value)
615 void Rotor::setTranslift(float value)
620 void Rotor::setSharedFlapHinge(bool s)
622 _shared_flap_hinge=s;
625 void Rotor::setBalance(float b)
630 void Rotor::setC2(float value)
635 void Rotor::setStepspersecond(float steps)
637 _stepspersecond=steps;
640 void Rotor::setRPM(float value)
645 void Rotor::setPhiNull(float value)
650 void Rotor::setRelLenHinge(float value)
652 _rel_len_hinge=value;
655 void Rotor::setAlphaoutput(int i, const char *text)
657 strncpy(_alphaoutput[i],text,255);
660 void Rotor::setName(const char *text)
662 strncpy(_name,text,256);//256: some space needed for settings
665 void Rotor::setCcw(int ccw)
670 void Rotor::setNotorque(int value)
675 void Rotor::setRelLenTeeterHinge(float f)
677 _rellenteeterhinge=f;
680 void Rotor::setTeeterdamp(float f)
685 void Rotor::setMaxteeterdamp(float f)
690 void Rotor::setGlobalGround(double *global_ground, float* global_vel)
693 for(i=0; i<4; i++) _global_ground[i] = global_ground[i];
696 void Rotor::setParameter(char *parametername, float value)
698 #define p(a,b) if (strcmp(parametername,#a)==0) _##a = (value * (b)); else
700 p(translift_maxfactor,1)
701 p(ground_effect_constant,1)
702 p(vortex_state_lift_factor,1)
710 p(number_of_segments,1)
712 p(rel_len_where_incidence_is_measured,1)
715 p(airfoil_incidence_no_lift,pi/180.)
716 p(rel_len_blade_start,1)
717 p(incidence_stall_zero_speed,pi/180.)
718 p(incidence_stall_half_sonic_speed,pi/180.)
719 p(lift_factor_stall,1)
720 p(stall_change_over,pi/180.)
721 p(drag_factor_stall,1)
722 p(airfoil_lift_coefficient,1)
723 p(airfoil_drag_coefficient0,1)
724 p(airfoil_drag_coefficient1,1)
726 p(rotor_correction_factor,1)
727 SG_LOG(SG_INPUT, SG_ALERT,
728 "internal error in parameter set up for rotor: '" <<
729 parametername <<"'" << endl);
733 float Rotor::getLiftFactor()
738 void Rotorgear::setRotorBrake(float lval)
740 lval = Math::clamp(lval, 0, 1);
744 void Rotor::setTiltYaw(float lval)
746 lval = Math::clamp(lval, -1, 1);
747 _tilt_yaw = _min_tilt_yaw+(lval+1)/2*(_max_tilt_yaw-_min_tilt_yaw);
748 _directions_and_postions_dirty = true;
751 void Rotor::setTiltPitch(float lval)
753 lval = Math::clamp(lval, -1, 1);
754 _tilt_pitch = _min_tilt_pitch+(lval+1)/2*(_max_tilt_pitch-_min_tilt_pitch);
755 _directions_and_postions_dirty = true;
758 void Rotor::setTiltRoll(float lval)
760 lval = Math::clamp(lval, -1, 1);
761 _tilt_roll = _min_tilt_roll+(lval+1)/2*(_max_tilt_roll-_min_tilt_roll);
762 _directions_and_postions_dirty = true;
765 void Rotor::setCollective(float lval)
767 lval = Math::clamp(lval, -1, 1);
769 _collective=_min_pitch+(lval+1)/2*(_max_pitch-_min_pitch);
770 for(i=0; i<_rotorparts.size(); i++) {
771 ((Rotorpart*)_rotorparts.get(i))->setCollective(_collective);
775 void Rotor::setCyclicele(float lval,float rval)
777 lval = Math::clamp(lval, -1, 1);
778 _cyclicele=_mincyclicele+(lval+1)/2*(_maxcyclicele-_mincyclicele);
781 void Rotor::setCyclicail(float lval,float rval)
783 lval = Math::clamp(lval, -1, 1);
785 _cyclicail=-(_mincyclicail+(lval+1)/2*(_maxcyclicail-_mincyclicail));
788 void Rotor::getPosition(float* out)
791 for(i=0; i<3; i++) out[i] = _base[i];
794 void Rotor::calcLiftFactor(float* v, float rho, State *s)
796 //calculates _lift_factor, which is a foactor for the lift of the rotor
798 //- ground effect (_f_ge)
799 //- vortex state (_f_vs)
800 //- translational lift (_f_tl)
805 // Calculate ground effect
806 _f_ge=1+_diameter/_ground_effect_altitude*_ground_effect_constant;
808 // Now calculate translational lift
809 float v_vert = Math::dot3(v,_normal);
811 Math::cross3(v,_normal,help);
812 float v_horiz = Math::mag3(help);
813 _f_tl = ((1-Math::pow(2.7183,-v_horiz/_translift_ve))
814 *(_translift_maxfactor-1)+1)/_translift_maxfactor;
816 _lift_factor = _f_ge*_f_tl*_f_vs;
818 //store the gravity direction
819 Glue::geodUp(s->pos, _grav_direction);
822 void Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s)
824 _ground_effect_altitude=findGroundEffectAltitude(ground_cb,s,
825 _groundeffectpos[0],_groundeffectpos[1],
826 _groundeffectpos[2],_groundeffectpos[3]);
827 testForRotorGroundContact(ground_cb,s);
830 void Rotor::testForRotorGroundContact(Ground * ground_cb,State *s)
833 for (i=0;i<_num_ground_contact_pos;i++)
836 s->posLocalToGlobal(_ground_contact_pos[i], pt);
838 // Ask for the ground plane in the global coordinate system
839 double global_ground[4];
841 ground_cb->getGroundPlane(pt, global_ground, global_vel);
842 // find h, the distance to the ground
843 // The ground plane transformed to the local frame.
845 s->planeGlobalToLocal(global_ground, ground);
847 h = ground[3] - Math::dot3(_ground_contact_pos[i], ground);
848 // Now h is the distance from _ground_contact_pos[i] to ground
851 _balance1 -= (-h)/_diameter/_num_ground_contact_pos;
852 _balance1 = (_balance1<-1)?-1:_balance1;
856 float Rotor::findGroundEffectAltitude(Ground * ground_cb,State *s,
857 float *pos0,float *pos1,float *pos2,float *pos3,
858 int iteration,float a0,float a1,float a2,float a3)
872 Math::add3(p[0],p[2],p[4]);
873 Math::mul3(0.5,p[4],p[4]);//the center
875 float mina=100*_diameter;
877 for (int i=0;i<5;i++)
879 if (a[i]==-1)//in the first iteration,(iteration==0) no height is
880 //passed to this function, these missing values are
884 s->posLocalToGlobal(p[i], pt);
886 // Ask for the ground plane in the global coordinate system
887 double global_ground[4];
889 ground_cb->getGroundPlane(pt, global_ground, global_vel);
890 // find h, the distance to the ground
891 // The ground plane transformed to the local frame.
893 s->planeGlobalToLocal(global_ground, ground);
895 a[i] = ground[3] - Math::dot3(p[i], ground);
896 // Now a[i] is the distance from p[i] to ground
902 if (mina>=10*_diameter)
903 return mina; //the ground effect will be zero
905 //check if further recursion is neccessary
906 //if the height does not differ more than 20%, than
907 //we can return then mean height, if not split
908 //zhe square to four parts and calcualte the height
910 //suma * 0.2 is the mean
911 //0.15 is the maximum allowed difference from the mean
912 //to the height at the center
914 ||(Math::abs(suma*0.2-a[4])<(0.15*0.2*suma*(1<<iteration))))
917 float pc[4][3],ac[4]; //pc[i]=center of pos[i] and pos[(i+1)&3]
918 for (int i=0;i<4;i++)
920 Math::add3(p[i],p[(i+1)&3],pc[i]);
921 Math::mul3(0.5,pc[i],pc[i]);
923 s->posLocalToGlobal(pc[i], pt);
925 // Ask for the ground plane in the global coordinate system
926 double global_ground[4];
928 ground_cb->getGroundPlane(pt, global_ground, global_vel);
929 // find h, the distance to the ground
930 // The ground plane transformed to the local frame.
932 s->planeGlobalToLocal(global_ground, ground);
934 ac[i] = ground[3] - Math::dot3(p[i], ground);
935 // Now ac[i] is the distance from pc[i] to ground
938 (findGroundEffectAltitude(ground_cb,s,p[0],pc[1],p[4],pc[3],
939 iteration+1,a[0],ac[0],a[4],ac[3])
940 +findGroundEffectAltitude(ground_cb,s,p[1],pc[0],p[4],pc[1],
941 iteration+1,a[1],ac[0],a[4],ac[1])
942 +findGroundEffectAltitude(ground_cb,s,p[2],pc[1],p[4],pc[2],
943 iteration+1,a[2],ac[1],a[4],ac[2])
944 +findGroundEffectAltitude(ground_cb,s,p[3],pc[2],p[4],pc[3],
945 iteration+1,a[3],ac[2],a[4],ac[3])
949 void Rotor::getDownWash(float *pos, float *v_heli, float *downwash)
951 float pos2rotor[3],tmp[3];
952 Math::sub3(_base,pos,pos2rotor);
953 float dist=Math::dot3(pos2rotor,_normal_with_yaw_roll);
954 //calculate incidence at 0.7r;
955 float inc = _collective+_twist *0.7
956 - _twist*_rel_len_where_incidence_is_measured;
959 if (dist<0) // we are not in the downwash region
961 downwash[0]=downwash[1]=downwash[2]=0.;
965 //calculate the mean downwash speed directly beneath the rotor disk
966 float v1bar = Math::sin(inc) *_omega * 0.35 * _diameter * 0.8;
967 //0.35 * d = 0.7 *r, a good position to calcualte the mean downwashd
968 //0.8 the slip of the rotor.
970 //calculate the time the wind needed from thr rotor to here
971 if (v1bar< 1) v1bar = 1;
972 float time=dist/v1bar;
974 //calculate the pos2rotor, where the rotor was, "time" ago
975 Math::mul3(time,v_heli,tmp);
976 Math::sub3(pos2rotor,tmp,pos2rotor);
978 //and again calculate dist
979 dist=Math::dot3(pos2rotor,_normal_with_yaw_roll);
980 //missing the normal is offen not pointing to the normal of the rotor
981 //disk. Rotate the normal by yaw and tilt angle
985 if (dist<0) // we are not in the downwash region
987 downwash[0]=downwash[1]=downwash[2]=0.;
990 //of course this could be done in a runge kutta integrator, but it's such
991 //a approximation that I beleave, it would'nt be more realistic
993 //calculate the dist to the rotor-axis
994 Math::cross3(pos2rotor,_normal_with_yaw_roll,tmp);
995 float r= Math::mag3(tmp);
996 //calculate incidence at r;
997 float rel_r = r *2 /_diameter;
998 float inc_r = _collective+_twist * r /_diameter * 2
999 - _twist*_rel_len_where_incidence_is_measured;
1001 //calculate the downwash speed directly beneath the rotor disk
1004 v1 = Math::sin(inc_r) *_omega * r * 0.8;
1006 //calcualte the downwash speed in a distance "dist" to the rotor disc,
1007 //for large dist. The speed is assumed do follow a gausian distribution
1008 //with sigma increasing with dist^2:
1009 //sigma is assumed to be half of the rotor diameter directly beneath the
1010 //disc and is assumed to the rotor diameter at dist = (diameter * sqrt(2))
1012 float sigma=_diameter/2 + dist * dist / _diameter /4.;
1013 float v2 = v1bar*_diameter/ (Math::sqrt(2 * pi) * sigma)
1014 * Math::pow(2.7183,-.5*r*r/(sigma*sigma))*_diameter/2/sigma;
1016 //calculate the weight of the two downwash velocities.
1017 //Directly beneath the disc it is v1, far away it is v2
1018 float g = Math::pow(2.7183,-2*dist/_diameter);
1019 //at dist = rotor radius it is assumed to be 1/e * v1 + (1-1/e)* v2
1021 float v = g * v1 + (1-g) * v2;
1022 Math::mul3(-v,_normal_with_yaw_roll,downwash);
1023 //the downwash is calculated in the opposite direction of the normal
1026 void Rotor::euler2orient(float roll, float pitch, float hdg, float* out)
1028 // the Glue::euler2orient, inverts y<z due to different bases
1029 // therefore the negation of all "y" and "z" coeffizients
1030 Glue::euler2orient(roll,pitch,hdg,out);
1031 for (int i=3;i<9;i++) out[i]*=-1.0;
1035 void Rotor::updateDirectionsAndPositions(float *rot)
1037 if (!_directions_and_postions_dirty)
1039 rot[0]=rot[1]=rot[2]=0;
1042 rot[0]=_old_tilt_roll-_tilt_roll;
1043 rot[1]=_old_tilt_pitch-_tilt_pitch;
1044 rot[2]=_old_tilt_yaw-_tilt_yaw;
1045 _old_tilt_roll=_tilt_roll;
1046 _old_tilt_pitch=_tilt_pitch;
1047 _old_tilt_yaw=_tilt_yaw;
1049 euler2orient(_tilt_roll, _tilt_pitch, _tilt_yaw, orient);
1053 Math::sub3(_base,_tilt_center,base);
1054 Math::vmul33(orient, base, base);
1055 Math::add3(base,_tilt_center,base);
1056 Math::vmul33(orient, _forward, forward);
1057 Math::vmul33(orient, _normal, normal);
1059 #define _forward forward
1060 #define _normal normal
1061 float directions[5][3];
1062 //pointing forward, right, ... the 5th is ony for calculation
1063 directions[0][0]=_forward[0];
1064 directions[0][1]=_forward[1];
1065 directions[0][2]=_forward[2];
1070 Math::cross3(directions[i-1],_normal,directions[i]);
1072 Math::cross3(_normal,directions[i-1],directions[i]);
1074 Math::set3(directions[4],directions[0]);
1075 // now directions[0] is perpendicular to the _normal.and has a length
1076 // of 1. if _forward is already normalized and perpendicular to the
1077 // normal, directions[0] will be the same
1078 //_num_ground_contact_pos=(_number_of_parts<16)?_number_of_parts:16;
1079 for (i=0;i<_num_ground_contact_pos;i++)
1082 float s = Math::sin(pi*2*i/_num_ground_contact_pos);
1083 float c = Math::cos(pi*2*i/_num_ground_contact_pos);
1084 Math::mul3(c*_diameter*0.5,directions[0],_ground_contact_pos[i]);
1085 Math::mul3(s*_diameter*0.5,directions[1],help);
1086 Math::add3(help,_ground_contact_pos[i],_ground_contact_pos[i]);
1087 Math::add3(_base,_ground_contact_pos[i],_ground_contact_pos[i]);
1091 Math::mul3(_diameter*0.7,directions[i],_groundeffectpos[i]);
1092 Math::add3(_base,_groundeffectpos[i],_groundeffectpos[i]);
1094 for (i=0;i<_number_of_parts;i++)
1096 Rotorpart* rp = getRotorpart(i);
1097 float lpos[3],lforceattac[3],lspeed[3],dirzentforce[3];
1098 float s = Math::sin(2*pi*i/_number_of_parts);
1099 float c = Math::cos(2*pi*i/_number_of_parts);
1100 float sp = Math::sin(float(2*pi*i/_number_of_parts-pi/2.+_phi));
1101 float cp = Math::cos(float(2*pi*i/_number_of_parts-pi/2.+_phi));
1102 float direction[3],nextdirection[3],help[3],direction90deg[3];
1103 float rotorpartmass = _weight_per_blade*_number_of_blades/_number_of_parts*.453;
1104 float speed=_rotor_rpm/60*_diameter*_rel_blade_center*pi;
1105 float lentocenter=_diameter*_rel_blade_center*0.5;
1106 float lentoforceattac=_diameter*_rel_len_hinge*0.5;
1107 float zentforce=rotorpartmass*speed*speed/lentocenter;
1109 Math::mul3(c ,directions[0],help);
1110 Math::mul3(s ,directions[1],direction);
1111 Math::add3(help,direction,direction);
1113 Math::mul3(c ,directions[1],help);
1114 Math::mul3(s ,directions[2],direction90deg);
1115 Math::add3(help,direction90deg,direction90deg);
1117 Math::mul3(cp ,directions[1],help);
1118 Math::mul3(sp ,directions[2],nextdirection);
1119 Math::add3(help,nextdirection,nextdirection);
1121 Math::mul3(lentocenter,direction,lpos);
1122 Math::add3(lpos,_base,lpos);
1123 Math::mul3(lentoforceattac,nextdirection,lforceattac);
1124 //nextdirection: +90deg (gyro)!!!
1126 Math::add3(lforceattac,_base,lforceattac);
1127 Math::mul3(speed,direction90deg,lspeed);
1128 Math::mul3(1,nextdirection,dirzentforce);
1129 rp->setPosition(lpos);
1130 rp->setNormal(_normal);
1131 rp->setZentipetalForce(zentforce);
1132 rp->setPositionForceAttac(lforceattac);
1133 rp->setSpeed(lspeed);
1134 rp->setDirectionofZentipetalforce(dirzentforce);
1135 rp->setDirectionofRotorPart(direction);
1140 _directions_and_postions_dirty=false;
1143 void Rotor::compile()
1145 // Have we already been compiled?
1146 if(_rotorparts.size() != 0) return;
1148 //rotor is divided into _number_of_parts parts
1149 //each part is calcualted at _number_of_segments points
1152 //and make it a factor of 4
1153 _number_of_parts=(int(Math::clamp(_number_of_parts,4,256))>>2)<<2;
1155 _dynamic=_dynamic*(1/ //inverse of the time
1156 ( (60/_rotor_rpm)/4 //for rotating 90 deg
1157 +(60/_rotor_rpm)/(2*_number_of_blades) //+ meantime a rotorblade
1158 //will pass a given point
1160 //normalize the directions
1161 Math::unit3(_forward,_forward);
1162 Math::unit3(_normal,_normal);
1163 _num_ground_contact_pos=(_number_of_parts<16)?_number_of_parts:16;
1164 float rotorpartmass = _weight_per_blade*_number_of_blades/_number_of_parts*.453;
1165 //was pounds -> now kg
1167 _torque_of_inertia = 1/12. * ( _number_of_parts * rotorpartmass) * _diameter
1168 * _diameter * _rel_blade_center * _rel_blade_center /(0.5*0.5);
1169 float speed=_rotor_rpm/60*_diameter*_rel_blade_center*pi;
1170 float lentocenter=_diameter*_rel_blade_center*0.5;
1171 float lentoforceattac=_diameter*_rel_len_hinge*0.5;
1172 float zentforce=rotorpartmass*speed*speed/lentocenter;
1173 float pitchaforce=_force_at_pitch_a/_number_of_parts*.453*9.81;
1174 // was pounds of force, now N, devided by _number_of_parts
1175 //(so its now per rotorpart)
1177 float torque0=0,torquemax=0,torqueb=0;
1178 float omega=_rotor_rpm/60*2*pi;
1180 float omega0=omega*Math::sqrt(1/(1-_rel_len_hinge));
1181 float delta_theoretical=pitchaforce/(_pitch_a*omega*lentocenter*2*rotorpartmass);
1182 _delta*=delta_theoretical;
1184 float relamp=(omega*omega/(2*_delta*Math::sqrt(sqr(omega0*omega0-omega*omega)
1185 +4*_delta*_delta*omega*omega)))*_cyclic_factor;
1186 float relamp_theoretical=(omega*omega/(2*delta_theoretical*Math::sqrt(sqr(omega0*omega0-omega*omega)
1187 +4*delta_theoretical*delta_theoretical*omega*omega)))*_cyclic_factor;
1188 _phi=Math::acos(_rel_len_hinge);
1189 _phi-=Math::atan(_delta3);
1192 torque0=_power_at_pitch_0/_number_of_parts*1000/omega;
1193 // f*r=p/w ; p=f*s/t; r=s/t/w ; r*w*t = s
1194 torqueb=_power_at_pitch_b/_number_of_parts*1000/omega;
1195 torquemax=_power_at_pitch_b/_number_of_parts*1000/omega/_pitch_b*_max_pitch;
1205 Rotorpart* rps[256];
1207 for (i=0;i<_number_of_parts;i++)
1209 Rotorpart* rp=rps[i]=newRotorpart(zentforce,pitchaforce,_delta3,rotorpartmass,
1210 _translift,_rel_len_hinge,lentocenter);
1211 int k = i*4/_number_of_parts;
1212 rp->setAlphaoutput(_alphaoutput[k&1?k:(_ccw?k^2:k)],0);
1213 rp->setAlphaoutput(_alphaoutput[4+(k&1?k:(_ccw?k^2:k))],1+(k>1));
1214 _rotorparts.add(rp);
1215 rp->setTorque(torquemax,torque0);
1216 rp->setRelamp(relamp);
1217 rp->setTorqueOfInertia(_torque_of_inertia/_number_of_parts);
1218 rp->setDirection(2*pi*i/_number_of_parts);
1220 for (i=0;i<_number_of_parts;i++)
1222 rps[i]->setlastnextrp(rps[(i-1+_number_of_parts)%_number_of_parts],
1223 rps[(i+1)%_number_of_parts],
1224 rps[(i+_number_of_parts/2)%_number_of_parts],
1225 rps[(i-_number_of_parts/4+_number_of_parts)%_number_of_parts],
1226 rps[(i+_number_of_parts/4)%_number_of_parts]);
1229 updateDirectionsAndPositions(drot);
1230 for (i=0;i<_number_of_parts;i++)
1232 rps[i]->setCompiled();
1234 float lift[4],torque[4], v_wind[3];
1235 v_wind[0]=v_wind[1]=v_wind[2]=0;
1236 rps[0]->setOmega(_omegan);
1238 if (_airfoil_lift_coefficient==0)
1240 //calculate the lift and drag coefficients now
1244 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_a,0,0,
1245 &(torque[0]),&(lift[0])); //max_pitch a
1246 _liftcoef = pitchaforce/lift[0];
1249 rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,&(torque[0]),&(lift[0]));
1254 rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,&(torque[1]),&(lift[1]));
1259 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,&(torque[2]),&(lift[2]));
1264 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,&(torque[3]),&(lift[3]));
1269 _dragcoef1=torque0/torque[1];
1270 _dragcoef0=(torqueb-_dragcoef1*torque[3])/torque[2];
1274 _dragcoef1=(torque0/torque[0]-torqueb/torque[2])
1275 /(torque[1]/torque[0]-torque[3]/torque[2]);
1276 _dragcoef0=(torqueb-_dragcoef1*torque[3])/torque[2];
1281 _liftcoef=_airfoil_lift_coefficient/_number_of_parts*_number_of_blades;
1282 _dragcoef0=_airfoil_drag_coefficient0/_number_of_parts*_number_of_blades*_c2;
1283 _dragcoef1=_airfoil_drag_coefficient1/_number_of_parts*_number_of_blades*_c2;
1287 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_a,0,0,
1288 &(torque[0]),&(lift[0])); //pitch a
1289 rps[0]->calculateAlpha(v_wind,rho_null,_pitch_b,0,0,
1290 &(torque[1]),&(lift[1])); //pitch b
1291 rps[0]->calculateAlpha(v_wind,rho_null,0,0,0,
1292 &(torque[3]),&(lift[3])); //pitch 0
1293 SG_LOG(SG_GENERAL, SG_INFO,
1294 "Rotor: coefficients for airfoil:" << endl << setprecision(6)
1295 << " drag0: " << _dragcoef0*_number_of_parts/_number_of_blades/_c2
1296 << " drag1: " << _dragcoef1*_number_of_parts/_number_of_blades/_c2
1297 << " lift: " << _liftcoef*_number_of_parts/_number_of_blades
1299 << "at 10 deg:" << endl
1300 << "drag: " << (Math::sin(10./180*pi)*_dragcoef1+_dragcoef0)
1301 *_number_of_parts/_number_of_blades/_c2
1302 << " lift: " << Math::sin(10./180*pi)*_liftcoef*_number_of_parts/_number_of_blades
1304 << "Some results (Pitch [degree], Power [kW], Lift [N])" << endl
1305 << 0.0f << "deg " << Math::abs(torque[3]*_number_of_parts*_omegan/1000) << "kW "
1306 << lift[3]*_number_of_parts << endl
1307 << _pitch_a*180/pi << "deg " << Math::abs(torque[0]*_number_of_parts*_omegan/1000)
1308 << "kW " << lift[0]*_number_of_parts << endl
1309 << _pitch_b*180/pi << "deg " << Math::abs(torque[1]*_number_of_parts*_omegan/1000)
1310 << "kW " << lift[1]*_number_of_parts << endl << endl );
1312 //first calculation of relamp is wrong
1313 //it used pitchaforce, but this was unknown and
1314 //on the default value
1315 _delta*=lift[0]/pitchaforce;
1316 relamp=(omega*omega/(2*_delta*Math::sqrt(sqr(omega0*omega0-omega*omega)
1317 +4*_delta*_delta*omega*omega)))*_cyclic_factor;
1318 for (i=0;i<_number_of_parts;i++)
1320 rps[i]->setRelamp(relamp);
1322 rps[0]->setOmega(0);
1329 //tie the properties
1330 /* After reset these values are totally wrong. I have to find out why
1331 SGPropertyNode * node = fgGetNode("/rotors", true)->getNode(_name,true);
1332 node->tie("balance_ext",SGRawValuePointer<float>(&_balance2),false);
1333 node->tie("balance_int",SGRawValuePointer<float>(&_balance1));
1337 std::ostream & operator<<(std::ostream & out, Rotor& r)
1339 #define i(x) << #x << ":" << r.x << endl
1340 #define iv(x) << #x << ":" << r.x[0] << ";" << r.x[1] << ";" <<r.x[2] << ";" << endl
1341 out << "Writing Info on Rotor "
1344 i(_omega) i(_omegan) i(_omegarel) i(_ddt_omega) i(_omegarelneu)
1347 i( _airfoil_incidence_no_lift)
1349 i( _airfoil_lift_coefficient)
1350 i( _airfoil_drag_coefficient0)
1351 i( _airfoil_drag_coefficient1)
1353 i( _number_of_segments)
1354 i( _number_of_parts)
1356 iv( _groundeffectpos[0])iv( _groundeffectpos[1])iv( _groundeffectpos[2])iv( _groundeffectpos[3])
1357 i( _ground_effect_altitude)
1359 iv( _normal_with_yaw_roll)
1362 i( _number_of_blades)
1363 i( _weight_per_blade)
1364 i( _rel_blade_center)
1367 i( _force_at_pitch_a)
1369 i( _power_at_pitch_0)
1370 i( _power_at_pitch_b)
1387 i( _alphamin) i(_alphamax) i(_alpha0) i(_alpha0factor)
1388 i( _teeterdamp) i(_maxteeterdamp)
1389 i( _rellenteeterhinge)
1391 i( _translift_maxfactor)
1392 i( _ground_effect_constant)
1393 i( _vortex_state_lift_factor)
1394 i( _vortex_state_c1)
1395 i( _vortex_state_c2)
1396 i( _vortex_state_c3)
1397 i( _vortex_state_e1)
1398 i( _vortex_state_e2)
1399 i( _vortex_state_e3)
1400 i( _lift_factor) i(_f_ge) i(_f_vs) i(_f_tl)
1405 i( _twist) //outer incidence = inner inner incidence + _twist
1406 i( _rel_len_where_incidence_is_measured)
1407 i( _torque_of_inertia)
1408 i( _rel_len_blade_start)
1409 i( _incidence_stall_zero_speed)
1410 i( _incidence_stall_half_sonic_speed)
1411 i( _lift_factor_stall)
1412 i( _stall_change_over)
1413 i( _drag_factor_stall)
1420 i( _cyclic_factor) <<endl;
1422 for(j=0; j<r._rotorparts.size(); j++) {
1423 out << *((Rotorpart*)r._rotorparts.get(j));
1430 void Rotor:: writeInfo()
1433 std::ostringstream buffer;
1435 FILE*f=fopen("c:\\fgmsvc\\bat\\log.txt","at");
1436 if (!f) f=fopen("c:\\fgmsvc\\bat\\log.txt","wt");
1439 fprintf(f,"%s",(const char *)buffer.str().c_str());
1444 Rotorpart* Rotor::newRotorpart(float zentforce,float maxpitchforce,
1445 float delta3,float mass,float translift,float rellenhinge,float len)
1447 Rotorpart *r = new Rotorpart();
1448 r->setDelta3(delta3);
1449 r->setDynamic(_dynamic);
1450 r->setTranslift(_translift);
1453 r->setRelLenHinge(rellenhinge);
1454 r->setSharedFlapHinge(_shared_flap_hinge);
1455 r->setOmegaN(_omegan);
1456 r->setPhi(_phi_null);
1457 r->setAlpha0(_alpha0);
1458 r->setAlphamin(_alphamin);
1459 r->setAlphamax(_alphamax);
1460 r->setAlpha0factor(_alpha0factor);
1462 r->setDiameter(_diameter);
1464 #define p(a) r->setParameter(#a,_##a);
1466 p(number_of_segments)
1467 p(rel_len_where_incidence_is_measured)
1468 p(rel_len_blade_start)
1469 p(rotor_correction_factor)
1474 void Rotor::interp(float* v1, float* v2, float frac, float* out)
1476 out[0] = v1[0] + frac*(v2[0]-v1[0]);
1477 out[1] = v1[1] + frac*(v2[1]-v1[1]);
1478 out[2] = v1[2] + frac*(v2[2]-v1[2]);
1481 void Rotorgear::initRotorIteration(float *lrot,float dt)
1485 if (!_rotors.size()) return;
1486 Rotor* r0 = (Rotor*)_rotors.get(0);
1487 omegarel= r0->getOmegaRelNeu();
1488 for(i=0; i<_rotors.size(); i++) {
1489 Rotor* r = (Rotor*)_rotors.get(i);
1490 r->inititeration(dt,omegarel,0,lrot);
1494 void Rotorgear::calcForces(float* torqueOut)
1497 torqueOut[0]=torqueOut[1]=torqueOut[2]=0;
1498 // check,<if the engine can handle the torque of the rotors.
1499 // If not reduce the torque to the fueselage and change rotational
1500 // speed of the rotors instead
1503 float omegarel,omegan;
1504 Rotor* r0 = (Rotor*)_rotors.get(0);
1505 omegarel= r0->getOmegaRel();
1507 float total_torque_of_inertia=0;
1508 float total_torque=0;
1509 for(i=0; i<_rotors.size(); i++) {
1510 Rotor* r = (Rotor*)_rotors.get(i);
1511 omegan=r->getOmegan();
1512 total_torque_of_inertia+=r->getTorqueOfInertia()*omegan*omegan;
1513 //FIXME: this is constant, so this can be done in compile
1515 total_torque+=r->getTorque()*omegan;
1517 float max_torque_of_engine=0;
1518 SGPropertyNode * node=fgGetNode("/rotors/gear", true);
1521 max_torque_of_engine=_max_power_engine*_max_rel_torque;
1522 float df=_target_rel_rpm-omegarel;
1523 df/=_engine_prop_factor;
1524 df = Math::clamp(df, 0, 1);
1525 max_torque_of_engine = df * _max_power_engine*_max_rel_torque;
1529 float rel_torque_engine=1;
1530 if (total_torque<=0)
1531 rel_torque_engine=0;
1533 if (max_torque_of_engine>0)
1534 rel_torque_engine=1/max_torque_of_engine*total_torque;
1536 rel_torque_engine=0;
1538 //add the rotor brake and the gear fritcion
1540 if (r0->_rotorparts.size()) dt=((Rotorpart*)r0->_rotorparts.get(0))->getDt();
1542 float rotor_brake_torque;
1543 rotor_brake_torque=_rotorbrake*_max_power_rotor_brake+_rotorgear_friction;
1544 //clamp it to the value you need to stop the rotor
1545 //to avod accelerate the rotor to neagtive rpm:
1546 rotor_brake_torque=Math::clamp(rotor_brake_torque,0,
1547 total_torque_of_inertia/dt*omegarel);
1548 max_torque_of_engine-=rotor_brake_torque;
1550 //change the rotation of the rotors
1551 if ((max_torque_of_engine<total_torque) //decreasing rotation
1552 ||((max_torque_of_engine>total_torque)&&(omegarel<_target_rel_rpm))
1553 //increasing rotation due to engine
1554 ||(total_torque<0) ) //increasing rotation due to autorotation
1556 _ddt_omegarel=(max_torque_of_engine-total_torque)/total_torque_of_inertia;
1557 if(max_torque_of_engine>total_torque)
1559 //check if the acceleration is due to the engine. If yes,
1560 //the engine self limits the accel.
1561 float lim1=-total_torque/total_torque_of_inertia;
1562 //accel. by autorotation
1564 if (lim1<_engine_accel_limit) lim1=_engine_accel_limit;
1565 //if the accel by autorotation greater than the max. engine
1566 //accel, then this is the limit, if not: the engine is the limit
1567 if (_ddt_omegarel>lim1) _ddt_omegarel=lim1;
1569 if (_ddt_omegarel>5.5)_ddt_omegarel=5.5;
1570 //clamp it to avoid overflow. Should never be reached
1571 if (_ddt_omegarel<-5.5)_ddt_omegarel=-5.5;
1573 if (_max_power_engine<0.001) {omegarel=1;_ddt_omegarel=0;}
1574 //for debug: negative or no maxpower will result
1575 //in permanet 100% rotation
1577 omegarel+=dt*_ddt_omegarel;
1579 if (omegarel>2.5) omegarel=2.5;
1580 //clamp it to avoid overflow. Should never be reached
1581 if (omegarel<-.5) omegarel=-.5;
1583 r0->setOmegaRelNeu(omegarel);
1584 //calculate the torque, which is needed to accelerate the rotors.
1585 //Add this additional torque to the body
1586 for(j=0; j<_rotors.size(); j++) {
1587 Rotor* r = (Rotor*)_rotors.get(j);
1588 for(i=0; i<r->_rotorparts.size(); i++) {
1589 float torque_scalar=0;
1590 Rotorpart* rp = (Rotorpart*)r->_rotorparts.get(i);
1592 rp->getAccelTorque(_ddt_omegarel,torque);
1593 Math::add3(torque,torqueOut,torqueOut);
1597 _total_torque_on_engine=total_torque+_ddt_omegarel*total_torque_of_inertia;
1601 void Rotorgear::addRotor(Rotor* rotor)
1607 void Rotorgear::compile()
1610 for(int j=0; j<_rotors.size(); j++) {
1611 Rotor* r = (Rotor*)_rotors.get(j);
1616 void Rotorgear::getDownWash(float *pos, float * v_heli, float *downwash)
1619 downwash[0]=downwash[1]=downwash[2]=0;
1620 for(int i=0; i<_rotors.size(); i++) {
1621 Rotor* ro = (Rotor*)_rotors.get(i);
1622 ro->getDownWash(pos,v_heli,tmp);
1623 Math::add3(downwash,tmp,downwash); // + downwash
1627 void Rotorgear::setParameter(char *parametername, float value)
1629 #define p(a,b) if (strcmp(parametername,#a)==0) _##a = (value * (b)); else
1630 p(max_power_engine,1000)
1631 p(engine_prop_factor,1)
1632 p(yasimdragfactor,1)
1633 p(yasimliftfactor,1)
1634 p(max_power_rotor_brake,1000)
1635 p(rotorgear_friction,1000)
1636 p(engine_accel_limit,0.01)
1637 SG_LOG(SG_INPUT, SG_ALERT,
1638 "internal error in parameter set up for rotorgear: '"
1639 << parametername <<"'" << endl);
1642 int Rotorgear::getValueforFGSet(int j,char *text,float *f)
1646 sprintf(text,"/rotors/gear/total-torque");
1647 *f=_total_torque_on_engine;
1651 Rotorgear::Rotorgear()
1656 _max_power_rotor_brake=1;
1657 _rotorgear_friction=1;
1658 _max_power_engine=1000*450;
1659 _engine_prop_factor=0.05f;
1663 _engine_accel_limit=0.05f;
1664 _total_torque_on_engine=0;
1669 Rotorgear::~Rotorgear()
1671 for(int i=0; i<_rotors.size(); i++)
1672 delete (Rotor*)_rotors.get(i);
1675 }; // namespace yasim