/* CAPD C++ file generated Automatically from HyLink */
 #include <iostream> 
 #include "capd/capdlib.h" 
 using namespace std; 
 using namespace capd; 
 
int getNextMode(int curMode, interval curModeTime);
main(){ 
 
  cout.precision(10);
  try{ 
 /* Differential equation for mode M1 Testing */ 
    IMap mode1("par:c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25,w,theta_bar;var:p,l,pe,i;fun:c1*(2*theta_bar*(c20*p*p+c21*p+c22)-c12*(c2+c3*w*p+c4*w*p*p+c5*w*w*p)),4*(c15+c16*c25*((1/c11)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c17*c25*c25*((1/c11)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))*((1/c11)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c18*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w)) +c19*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w))*c25*((1/c11)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))-l), c1*(2*c23*(theta_bar)*(c20*p*p+c21*p+c22)-(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w)),0;");
 /* Differential equation for mode M2 Testing */ 
    IMap mode2("par:c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25,w,theta_bar;var:p,l,pe,i;fun:c1*(2*theta_bar*(c20*p*p+c21*p+c22)-c12*(c2+c3*w*p+c4*w*p*p+c5*w*w*p)),4*(c15+c16*c25*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c17*c25*c25*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c18*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w)) +c19*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w)) *c25*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))-l), c1*(2*c23*(theta_bar)*(c20*p*p+c21*p+c22)-(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w)),c14*(c24*l-c11);");
/*Differential equations for */
    IMap mode3("par:c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25,w,theta_bar;var:p,l,pe,i;fun:c1*(2*theta_bar*(c20*p*p+c21*p+c22)-c12*(c2+c3*w*p+c4*w*p*p+c5*w*w*p)),4*(c15+c16*c25*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c17*c25*c25*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c18*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w)) +c19*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w)) *c25*((1/c11)*(1+i+c13*(c24*l-c11))*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))-l), c1*(2*c23*(theta_bar)*(c20*p*p+c21*p+c22)-(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w)),c14*(c24*l-c11);");
    IMap mode4("par:c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25,w,theta_bar;var:p,l,pe,i;fun:c1*(2*theta_bar*(c20*p*p+c21*p+c22)-c12*(c2+c3*w*p+c4*w*p*p+c5*w*w*p)),4*(c15+c16*c25*((1/12.5)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c17*c25*c25*((1/12.5)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))*((1/12.5)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))+c18*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w)) +c19*(c12*(c2+c3*w*p+c4*w*p*p+c5*p*w*w)) *c25*((1/12.5)*(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w))-l), c1*(2*c23*(theta_bar)*(c20*p*p+c21*p+c22)-(c2+c3*w*pe+c4*w*pe*pe+c5*pe*w*w)),0;");

    ITaylor* solvers[] = {
     new ITaylor(mode1,10),
     new ITaylor(mode2,10),
     new ITaylor(mode3,10),
     new ITaylor(mode4,10),
    };

    mode1.setParameter("c1",0.41328);
    mode1.setParameter("c2",-0.366);
    mode1.setParameter("c3",0.08979);
    mode1.setParameter("c4",-0.0337);
    mode1.setParameter("c5",0.0001);
    mode1.setParameter("c6",2.821);
    mode1.setParameter("c7",-0.05231);
    mode1.setParameter("c8",0.10299);
    mode1.setParameter("c9",-0.00063);
    mode1.setParameter("c10",1.0);
    mode1.setParameter("c11",14.7);
    mode1.setParameter("c12",0.9);
    mode1.setParameter("c13",0.04);
    mode1.setParameter("c14",0.14);
    mode1.setParameter("c15",13.893);
    mode1.setParameter("c16",-35.2518);
    mode1.setParameter("c17",20.7364);
    mode1.setParameter("c18",2.6287);
    mode1.setParameter("c19",-1.592);
    mode1.setParameter("c20",-2.3421);
    mode1.setParameter("c21",2.7799);
    mode1.setParameter("c22",-0.3273);
    mode1.setParameter("c23",1.0);
    mode1.setParameter("c24",1.0);
    mode1.setParameter("c25",1.0);
    mode1.setParameter("w",(1000*3.1415926535)/30);
    mode1.setParameter("theta_bar",247.0);

    mode2.setParameter("c1",0.41328);
    mode2.setParameter("c2",-0.366);
    mode2.setParameter("c3",0.08979);
    mode2.setParameter("c4",-0.0337);
    mode2.setParameter("c5",0.0001);
    mode2.setParameter("c6",2.821);
    mode2.setParameter("c7",-0.05231);
    mode2.setParameter("c8",0.10299);
    mode2.setParameter("c9",-0.00063);
    mode2.setParameter("c10",1.0);
    mode2.setParameter("c11",14.7);
    mode2.setParameter("c12",0.9);
    mode2.setParameter("c13",0.04);
    mode2.setParameter("c14",0.14);
    mode2.setParameter("c15",13.893);
    mode2.setParameter("c16",-35.2518);
    mode2.setParameter("c17",20.7364);
    mode2.setParameter("c18",2.6287);
    mode2.setParameter("c19",-1.592);
    mode2.setParameter("c20",-2.3421);
    mode2.setParameter("c21",2.7799);
    mode2.setParameter("c22",-0.3273);
    mode2.setParameter("c23",1.0);
    mode2.setParameter("c24",1.0);
    mode2.setParameter("c25",1.0);
    mode2.setParameter("w",(1000*3.1415926535)/30);
    mode2.setParameter("theta_bar",247.0);

    mode3.setParameter("c1",0.41328);
    mode3.setParameter("c2",-0.366);
    mode3.setParameter("c3",0.08979);
    mode3.setParameter("c4",-0.0337);
    mode3.setParameter("c5",0.0001);
    mode3.setParameter("c6",2.821);
    mode3.setParameter("c7",-0.05231);
    mode3.setParameter("c8",0.10299);
    mode3.setParameter("c9",-0.00063);
    mode3.setParameter("c10",1.0);
    mode3.setParameter("c11",14.7);
    mode3.setParameter("c12",0.9);
    mode3.setParameter("c13",0.04);
    mode3.setParameter("c14",0.14);
    mode3.setParameter("c15",13.893);
    mode3.setParameter("c16",-35.2518);
    mode3.setParameter("c17",20.7364);
    mode3.setParameter("c18",2.6287);
    mode3.setParameter("c19",-1.592);
    mode3.setParameter("c20",-2.3421);
    mode3.setParameter("c21",2.7799);
    mode3.setParameter("c22",-0.3273);
    mode3.setParameter("c23",1.0);
    mode3.setParameter("c24",1.0);
    mode3.setParameter("c25",1.0);
    mode3.setParameter("w",(1000*3.1415926535)/30);
    mode3.setParameter("theta_bar",10.0);

    mode4.setParameter("c1",0.41328);
    mode4.setParameter("c2",-0.366);
    mode4.setParameter("c3",0.08979);
    mode4.setParameter("c4",-0.0337);
    mode4.setParameter("c5",0.0001);
    mode4.setParameter("c6",2.821);
    mode4.setParameter("c7",-0.05231);
    mode4.setParameter("c8",0.10299);
    mode4.setParameter("c9",-0.00063);
    mode4.setParameter("c10",1.0);
    mode4.setParameter("c11",14.7);
    mode4.setParameter("c12",0.9);
    mode4.setParameter("c13",0.04);
    mode4.setParameter("c14",0.14);
    mode4.setParameter("c15",13.893);
    mode4.setParameter("c16",-35.2518);
    mode4.setParameter("c17",20.7364);
    mode4.setParameter("c18",2.6287);
    mode4.setParameter("c19",-1.592);
    mode4.setParameter("c20",-2.3421);
    mode4.setParameter("c21",2.7799);
    mode4.setParameter("c22",-0.3273);
    mode4.setParameter("c23",1.0);
    mode4.setParameter("c24",1.0);
    mode4.setParameter("c25",1.0);
    mode4.setParameter("w",(1000*3.1415926535)/30);
    mode4.setParameter("theta_bar",369.1);

    double initialState[4];
    IVector IState(4);
    int i;
    double absErr, relErr;
    double tstep,Gt;
    double curTime;
    int curMode, nextMode;
    std::cin >> curTime;
    for(i=0;i<4;i++){
      std::cin >> initialState[i]; IState[i] = initialState[i];
    }
    std::cin >> absErr >> relErr >> tstep >> Gt;
    std::cin >> curMode;
    nextMode = curMode;
    IVector SimState(4);
    IVector PrevSimState(4);
    PrevSimState = IState;
    interval tstepi(tstep);
    interval Gti(Gt);
    for(i=0;i<2;i++)
      solvers[i]->setStep(tstepi);
    interval currTime(curTime), currModeTime(0.0);
    C0HORect2Set SimSet(IState);
    while(currTime < Gti){
      std::cout << " " << currTime.leftBound();
      for(i=0;i<4;i++){
        std::cout << " " << PrevSimState[i].leftBound();
      }
      std::cout << endl;
      SimSet.move(*solvers[curMode-1]);
      SimState = (IVector)SimSet;
      currTime+=tstepi;currModeTime+=tstepi;
      std::cout << " " << currTime.leftBound();
      for(i=0;i<4;i++){
        std::cout << " " << SimState[i].leftBound();
      }
      std::cout << endl;
      PrevSimState = SimState;
      nextMode = getNextMode(curMode,currModeTime);
      if(nextMode != curMode){
        curMode = nextMode; currModeTime = 0.0;
      }
    }
  }catch(exception& e){
    cout << "Exception caught!" << e.what() << endl << endl;
  }
}

int getNextMode(int curMode, interval curModeTime){
  return curMode;
}
