1 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 --------- Copyright (C) 1999 Anthony K. Peden (apeden@earthlink.net) ---------
9 This program is free software; you can redistribute it and/or modify it under
10 the terms of the GNU Lesser General Public License as published by the Free Software
11 Foundation; either version 2 of the License, or (at your option) any later
14 This program is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
19 You should have received a copy of the GNU Lesser General Public License along with
20 this program; if not, write to the Free Software Foundation, Inc., 59 Temple
21 Place - Suite 330, Boston, MA 02111-1307, USA.
23 Further information about the GNU Lesser General Public License can also be found on
24 the world wide web at http://www.gnu.org.
27 --------------------------------------------------------------------------------
30 FUNCTIONAL DESCRIPTION
31 --------------------------------------------------------------------------------
33 This class takes the given set of IC's and finds the angle of attack, elevator,
34 and throttle setting required to fly steady level. This is currently for in-air
35 conditions only. It is implemented using an iterative, one-axis-at-a-time
38 // !!!!!!! BEWARE ALL YE WHO ENTER HERE !!!!!!!
40 /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
47 #include "models/FGAtmosphere.h"
48 #include "FGInitialCondition.h"
49 #include "models/FGAircraft.h"
50 #include "models/FGMassBalance.h"
51 #include "models/FGGroundReactions.h"
52 #include "models/FGInertial.h"
53 #include "models/FGAerodynamics.h"
54 #include "models/FGPropulsion.h"
55 #include "models/propulsion/FGEngine.h"
56 #include "math/FGColumnVector3.h"
59 #pragma warning (disable : 4786 4788)
66 static const char *IdSrc = "$Id: FGTrim.cpp,v 1.13 2010/04/23 17:23:40 dpculp Exp $";
67 static const char *IdHdr = ID_TRIM;
69 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71 FGTrim::FGTrim(FGFDMExec *FDMExec,TrimMode tt) {
75 max_sub_iterations=100;
77 A_Tolerance = Tolerance / 10;
91 if (debug_lvl & 2) cout << "Instantiated: FGTrim" << endl;
94 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96 FGTrim::~FGTrim(void) {
97 for(current_axis=0; current_axis<TrimAxes.size(); current_axis++) {
98 delete TrimAxes[current_axis];
100 delete[] sub_iterations;
103 if (debug_lvl & 2) cout << "Destroyed: FGTrim" << endl;
106 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108 void FGTrim::TrimStats() {
110 cout << endl << " Trim Statistics: " << endl;
111 cout << " Total Iterations: " << total_its << endl;
113 cout << " Sub-iterations:" << endl;
114 for (current_axis=0; current_axis<TrimAxes.size(); current_axis++) {
115 run_sum += TrimAxes[current_axis]->GetRunCount();
116 cout << " " << setw(5) << TrimAxes[current_axis]->GetStateName().c_str()
117 << ": " << setprecision(3) << sub_iterations[current_axis]
118 << " average: " << setprecision(5) << sub_iterations[current_axis]/double(total_its)
119 << " successful: " << setprecision(3) << successful[current_axis]
120 << " stability: " << setprecision(5) << TrimAxes[current_axis]->GetAvgStability()
123 cout << " Run Count: " << run_sum << endl;
127 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
129 void FGTrim::Report(void) {
130 cout << " Trim Results: " << endl;
131 for(current_axis=0; current_axis<TrimAxes.size(); current_axis++)
132 TrimAxes[current_axis]->AxisReport();
136 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138 void FGTrim::ClearStates(void) {
142 vector<FGTrimAxis*>::iterator iAxes;
143 iAxes = TrimAxes.begin();
144 while (iAxes != TrimAxes.end()) {
150 //cout << "TrimAxes.size(): " << TrimAxes.size() << endl;
153 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 bool FGTrim::AddState( State state, Control control ) {
160 vector <FGTrimAxis*>::iterator iAxes = TrimAxes.begin();
161 while (iAxes != TrimAxes.end()) {
163 if( ta->GetStateType() == state )
168 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,state,control));
169 delete[] sub_iterations;
172 sub_iterations=new double[TrimAxes.size()];
173 successful=new double[TrimAxes.size()];
174 solution=new bool[TrimAxes.size()];
179 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
181 bool FGTrim::RemoveState( State state ) {
186 vector <FGTrimAxis*>::iterator iAxes = TrimAxes.begin();
187 while (iAxes != TrimAxes.end()) {
189 if( ta->GetStateType() == state ) {
191 iAxes = TrimAxes.erase(iAxes);
198 delete[] sub_iterations;
201 sub_iterations=new double[TrimAxes.size()];
202 successful=new double[TrimAxes.size()];
203 solution=new bool[TrimAxes.size()];
208 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
210 bool FGTrim::EditState( State state, Control new_control ){
215 vector <FGTrimAxis*>::iterator iAxes = TrimAxes.begin();
216 while (iAxes != TrimAxes.end()) {
218 if( ta->GetStateType() == state ) {
219 TrimAxes.insert(iAxes,1,new FGTrimAxis(fdmex,fgic,state,new_control));
221 TrimAxes.erase(iAxes+1);
231 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
233 bool FGTrim::DoTrim(void) {
238 for(i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++){
239 fdmex->GetGroundReactions()->GetGearUnit(i)->SetReport(false);
242 fdmex->DisableOutput();
244 setEngineTrimMode(true);
246 fgic->SetPRadpsIC(0.0);
247 fgic->SetQRadpsIC(0.0);
248 fgic->SetRRadpsIC(0.0);
250 //clear the sub iterations counts & zero out the controls
251 for(current_axis=0;current_axis<TrimAxes.size();current_axis++) {
252 //cout << current_axis << " " << TrimAxes[current_axis]->GetStateName()
253 //<< " " << TrimAxes[current_axis]->GetControlName()<< endl;
254 if(TrimAxes[current_axis]->GetStateType() == tQdot) {
255 if(mode == tGround) {
256 TrimAxes[current_axis]->initTheta();
259 xlo=TrimAxes[current_axis]->GetControlMin();
260 xhi=TrimAxes[current_axis]->GetControlMax();
261 TrimAxes[current_axis]->SetControl((xlo+xhi)/2);
262 TrimAxes[current_axis]->Run();
263 //TrimAxes[current_axis]->AxisReport();
264 sub_iterations[current_axis]=0;
265 successful[current_axis]=0;
266 solution[current_axis]=false;
270 if(mode == tPullup ) {
271 cout << "Setting pitch rate and nlf... " << endl;
273 cout << "pitch rate done ... " << endl;
274 TrimAxes[0]->SetStateTarget(targetNlf);
275 cout << "nlf done" << endl;
276 } else if (mode == tTurn) {
278 //TrimAxes[0]->SetStateTarget(targetNlf);
283 for(current_axis=0;current_axis<TrimAxes.size();current_axis++) {
287 if(!solution[current_axis]) {
289 solution[current_axis]=true;
292 } else if(findInterval()) {
295 solution[current_axis]=false;
297 sub_iterations[current_axis]+=Nsub;
299 for(current_axis=0;current_axis<TrimAxes.size();current_axis++) {
300 //these checks need to be done after all the axes have run
301 if(Debug > 0) TrimAxes[current_axis]->AxisReport();
302 if(TrimAxes[current_axis]->InTolerance()) {
304 successful[current_axis]++;
309 if((axis_count == TrimAxes.size()-1) && (TrimAxes.size() > 1)) {
310 //cout << TrimAxes.size()-1 << " out of " << TrimAxes.size() << "!" << endl;
311 //At this point we can check the input limits of the failed axis
312 //and declare the trim failed if there is no sign change. If there
313 //is, keep going until success or max iteration count
315 //Oh, well: two out of three ain't bad
316 for(current_axis=0;current_axis<TrimAxes.size();current_axis++) {
317 //these checks need to be done after all the axes have run
318 if(!TrimAxes[current_axis]->InTolerance()) {
320 // special case this for now -- if other cases arise proper
321 // support can be added to FGTrimAxis
322 if( (gamma_fallback) &&
323 (TrimAxes[current_axis]->GetStateType() == tUdot) &&
324 (TrimAxes[current_axis]->GetControlType() == tThrottle)) {
325 cout << " Can't trim udot with throttle, trying flight"
326 << " path angle. (" << N << ")" << endl;
327 if(TrimAxes[current_axis]->GetState() > 0)
328 TrimAxes[current_axis]->SetControlToMin();
330 TrimAxes[current_axis]->SetControlToMax();
331 TrimAxes[current_axis]->Run();
332 delete TrimAxes[current_axis];
333 TrimAxes[current_axis]=new FGTrimAxis(fdmex,fgic,tUdot,
336 cout << " Sorry, " << TrimAxes[current_axis]->GetStateName()
337 << " doesn't appear to be trimmable" << endl;
339 trim_failed=true; //force the trim to fail
344 } //all-but-one check
346 if(N > max_iterations)
348 } while((axis_count < TrimAxes.size()) && (!trim_failed));
349 if((!trim_failed) && (axis_count >= TrimAxes.size())) {
352 cout << endl << " Trim successful" << endl;
356 cout << endl << " Trim failed" << endl;
358 for(i=0;i < fdmex->GetGroundReactions()->GetNumGearUnits();i++){
359 fdmex->GetGroundReactions()->GetGearUnit(i)->SetReport(true);
361 setEngineTrimMode(false);
362 fdmex->EnableOutput();
366 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
368 bool FGTrim::solve(void) {
370 double x1,x2,x3,f1,f2,f3,d,d0;
371 const double relax =0.9;
372 double eps=TrimAxes[current_axis]->GetSolverEps();
378 if( solutionDomain != 0) {
379 /* if(ahi > alo) { */
388 //max_sub_iterations=TrimAxes[current_axis]->GetIterationLimit();
389 while ( (TrimAxes[current_axis]->InTolerance() == false )
390 && (fabs(d) > eps) && (Nsub < max_sub_iterations)) {
393 x2=x1-d*d0*f1/(f3-f1);
394 TrimAxes[current_axis]->SetControl(x2);
395 TrimAxes[current_axis]->Run();
396 f2=TrimAxes[current_axis]->GetState();
398 cout << "FGTrim::solve Nsub,x1,x2,x3: " << Nsub << ", " << x1
399 << ", " << x2 << ", " << x3 << endl;
400 cout << " " << f1 << ", " << f2 << ", " << f3 << endl;
406 //cout << "Solution is between x1 and x2" << endl;
408 else if(f2*f3 <= 0.0) {
412 //cout << "Solution is between x2 and x3" << endl;
419 if(Nsub < max_sub_iterations) success=true;
424 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
426 produces an interval (xlo..xhi) on one side or the other of the current
427 control value in which a solution exists. This domain is, hopefully,
428 smaller than xmin..0 or 0..xmax and the solver will require fewer iterations
429 to find the solution. This is, hopefully, more efficient than having the
430 solver start from scratch every time. Maybe it isn't though...
431 This tries to take advantage of the idea that the changes from iteration to
432 iteration will be small after the first one or two top-level iterations.
434 assumes that changing the control will a produce significant change in the
435 accel i.e. checkLimits() has already been called.
437 if a solution is found above the current control, the function returns true
438 and xlo is set to the current control, xhi to the interval max it found, and
439 solutionDomain is set to 1.
440 if the solution lies below the current control, then the function returns
441 true and xlo is set to the interval min it found and xmax to the current
442 control. if no solution is found, then the function returns false.
445 in all cases, alo=accel(xlo) and ahi=accel(xhi) after the function exits.
446 no assumptions about the state of the sim after this function has run
449 bool FGTrim::findInterval(void) {
452 double current_control=TrimAxes[current_axis]->GetControl();
453 double current_accel=TrimAxes[current_axis]->GetState();;
454 double xmin=TrimAxes[current_axis]->GetControlMin();
455 double xmax=TrimAxes[current_axis]->GetControlMax();
456 double lastxlo,lastxhi,lastalo,lastahi;
458 step=0.025*fabs(xmax);
459 xlo=xhi=current_control;
460 alo=ahi=current_accel;
461 lastxlo=xlo;lastxhi=xhi;
462 lastalo=alo;lastahi=ahi;
468 if(xlo < xmin) xlo=xmin;
470 if(xhi > xmax) xhi=xmax;
471 TrimAxes[current_axis]->SetControl(xlo);
472 TrimAxes[current_axis]->Run();
473 alo=TrimAxes[current_axis]->GetState();
474 TrimAxes[current_axis]->SetControl(xhi);
475 TrimAxes[current_axis]->Run();
476 ahi=TrimAxes[current_axis]->GetState();
477 if(fabs(ahi-alo) <= TrimAxes[current_axis]->GetTolerance()) continue;
478 if(alo*ahi <=0) { //found interval with root
480 if(alo*current_accel <= 0) { //narrow interval down a bit
484 //xhi=current_control;
490 //xlo=current_control;
494 lastxlo=xlo;lastxhi=xhi;
495 lastalo=alo;lastahi=ahi;
496 if( !found && xlo==xmin && xhi==xmax ) continue;
498 cout << "FGTrim::findInterval: Nsub=" << Nsub << " Lo= " << xlo
499 << " Hi= " << xhi << " alo*ahi: " << alo*ahi << endl;
500 } while(!found && (Nsub <= max_sub_iterations) );
504 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
505 //checks to see which side of the current control value the solution is on
506 //and sets solutionDomain accordingly:
507 // 1 if solution is between the current and max
508 // -1 if solution is between the min and current
509 // 0 if there is no solution
511 //if changing the control produces no significant change in the accel then
512 //solutionDomain is set to zero and the function returns false
513 //if a solution is found, then xlo and xhi are set so that they bracket
514 //the solution, alo is set to accel(xlo), and ahi is set to accel(xhi)
515 //if there is no change or no solution then xlo=xmin, alo=accel(xmin) and
516 //xhi=xmax and ahi=accel(xmax)
517 //in all cases the sim is left such that the control=xmax and accel=ahi
519 bool FGTrim::checkLimits(void) {
521 double current_control=TrimAxes[current_axis]->GetControl();
522 double current_accel=TrimAxes[current_axis]->GetState();
523 xlo=TrimAxes[current_axis]->GetControlMin();
524 xhi=TrimAxes[current_axis]->GetControlMax();
526 TrimAxes[current_axis]->SetControl(xlo);
527 TrimAxes[current_axis]->Run();
528 alo=TrimAxes[current_axis]->GetState();
529 TrimAxes[current_axis]->SetControl(xhi);
530 TrimAxes[current_axis]->Run();
531 ahi=TrimAxes[current_axis]->GetState();
533 cout << "checkLimits() xlo,xhi,alo,ahi: " << xlo << ", " << xhi << ", "
534 << alo << ", " << ahi << endl;
536 solutionExists=false;
537 if(fabs(ahi-alo) > TrimAxes[current_axis]->GetTolerance()) {
538 if(alo*current_accel <= 0) {
543 } else if(current_accel*ahi < 0){
550 TrimAxes[current_axis]->SetControl(current_control);
551 TrimAxes[current_axis]->Run();
552 return solutionExists;
555 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
557 void FGTrim::setupPullup() {
559 g=fdmex->GetInertial()->gravity();
560 cgamma=cos(fgic->GetFlightPathAngleRadIC());
561 cout << "setPitchRateInPullup(): " << g << ", " << cgamma << ", "
562 << fgic->GetVtrueFpsIC() << endl;
563 q=g*(targetNlf-cgamma)/fgic->GetVtrueFpsIC();
564 cout << targetNlf << ", " << q << endl;
565 fgic->SetQRadpsIC(q);
566 cout << "setPitchRateInPullup() complete" << endl;
570 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
572 void FGTrim::setupTurn(void){
574 phi = fgic->GetPhiRadIC();
575 if( fabs(phi) > 0.001 && fabs(phi) < 1.56 ) {
576 targetNlf = 1 / cos(phi);
577 g = fdmex->GetInertial()->gravity();
578 psidot = g*tan(phi) / fgic->GetUBodyFpsIC();
579 cout << targetNlf << ", " << psidot << endl;
584 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
586 void FGTrim::updateRates(void){
587 if( mode == tTurn ) {
588 double phi = fgic->GetPhiRadIC();
589 double g = fdmex->GetInertial()->gravity();
591 if(fabs(phi) > 0.001 && fabs(phi) < 1.56 ) {
592 theta=fgic->GetThetaRadIC();
593 phi=fgic->GetPhiRadIC();
594 psidot = g*tan(phi) / fgic->GetUBodyFpsIC();
595 p=-psidot*sin(theta);
596 q=psidot*cos(theta)*sin(phi);
597 r=psidot*cos(theta)*cos(phi);
601 fgic->SetPRadpsIC(p);
602 fgic->SetQRadpsIC(q);
603 fgic->SetRRadpsIC(r);
604 } else if( mode == tPullup && fabs(targetNlf-1) > 0.01) {
606 g=fdmex->GetInertial()->gravity();
607 cgamma=cos(fgic->GetFlightPathAngleRadIC());
608 q=g*(targetNlf-cgamma)/fgic->GetVtrueFpsIC();
609 fgic->SetQRadpsIC(q);
613 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
615 void FGTrim::setDebug(void) {
616 if(debug_axis == tAll ||
617 TrimAxes[current_axis]->GetStateType() == debug_axis ) {
626 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
628 void FGTrim::setEngineTrimMode(bool mode) {
629 FGPropulsion* prop = fdmex->GetPropulsion();
630 for (unsigned int i=0; i<prop->GetNumEngines(); i++) {
631 prop->GetEngine(i)->SetTrimMode(mode);
635 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
637 void FGTrim::SetMode(TrimMode tt) {
643 cout << " Full Trim" << endl;
644 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tWdot,tAlpha ));
645 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tUdot,tThrottle ));
646 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tPitchTrim ));
647 //TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tHmgt,tBeta ));
648 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tVdot,tPhi ));
649 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tPdot,tAileron ));
650 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tRdot,tRudder ));
654 cout << " Longitudinal Trim" << endl;
655 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tWdot,tAlpha ));
656 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tUdot,tThrottle ));
657 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tPitchTrim ));
661 cout << " Ground Trim" << endl;
662 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tWdot,tAltAGL ));
663 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tTheta ));
664 //TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tPdot,tPhi ));
667 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tNlf,tAlpha ));
668 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tUdot,tThrottle ));
669 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tPitchTrim ));
670 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tHmgt,tBeta ));
671 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tVdot,tPhi ));
672 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tPdot,tAileron ));
673 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tRdot,tRudder ));
676 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tWdot,tAlpha ));
677 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tUdot,tThrottle ));
678 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tQdot,tPitchTrim ));
679 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tVdot,tBeta ));
680 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tPdot,tAileron ));
681 TrimAxes.push_back(new FGTrimAxis(fdmex,fgic,tRdot,tRudder ));
687 //cout << "TrimAxes.size(): " << TrimAxes.size() << endl;
688 sub_iterations=new double[TrimAxes.size()];
689 successful=new double[TrimAxes.size()];
690 solution=new bool[TrimAxes.size()];
693 //YOU WERE WARNED, BUT YOU DID IT ANYWAY.