#include <stdio.h>
#include <math.h>
#include   <string.h>
#include   <stdlib.h>
#include   <alloc.h>
#define TRUE 1
#define FALSE 0
#define PI acos(-1)

int k_index;
float th=0.3;
float u_plant,y_plant,y_plant_old,ydot,ydot_old;
void rmodel();
float rmodel1(float,float);
void main (void)
{
  int i;

  y_plant=ydot=0.0;
  printf("t=0:0.3:60;\n");
  for(k_index=1; k_index<=201; k_index++)
   {
     u_plant=1.0;
     rmodel();
     printf("y(%d)=%f;\n",k_index,y_plant);
   }
  printf("\n plot(t,y);\n");
}

void rmodel(void)
{
     float a1,a2;
     float b1,b2,b3,b4,c1,c2,c3,c4;

     a1 = y_plant;
     a2 = ydot;
     c1 = a2;
     b1 = rmodel1(a1,c1);
     c2 = a2+th/2*b1;
     b2 = rmodel1(a1+th/2*c1,c2);
     c3 = a2+th/2*b2;
     b3 = rmodel1(a1+th/2*c2,c3);
     c4 = a2+th/2*b3;
     b4 = rmodel1(a1+th*c3,c4);
     y_plant += th/6*(c1+2*c2+2*c3+c4);
     ydot += th/6*(b1+2*b2+2*b3+b4);
}

float rmodel1(float ae,float af)
{
      float den2;

      den2 = 19.54*u_plant-0.54*ae-0.4*af;
      return(den2);
}


/*  void forward()
{

 x1=0.9765*x1_old+0.2846*x2_old+0.047*u_plant;
 x2= -0.1537*x1_old+0.8911*x2_old+0.3073*u_plant;
 y_plant=x1;

 x1_old=x1;
 x2_old=x2;
}   */