//************************************************************************ 
// Forced oscillations: Pohl's pendulum with the classical RK4 method
//  J*phi''+K*phi'+D*phi-N*sin(phi)=F*sin(omega t)  (1)
//  We write down Eq.(1) as a system of three ODE's of the first order:
//   x'=y,
//   y'=-a*y-b*x+c*sin(x)+d*sin(omega*z),
//   z'=1,
// where a=K/J, b=D/J, c=N/J, d=F/J.
//************************************************************************
/* 
compile with	
gcc -o pohl pohl.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/

#include "/usr/local/dislin/dislin.h"
#include <stdio.h>
#include <math.h>

int step = 0.;
int max_steps = 6000.;

double h = 0.025;//time step

//initial conditions
double x = 0.0;
double y = 0.0;
double z = 0.0;

//constants
double a = 0.799;
double b=9.44;
double c=14.68;
double d=2.1;
//*************************************************
double omega=2.25; // control parameter, omega=2.5 (period-one), 2.32 (period-two), 2.3 (period four), 2.25 (chaos).
//*************************************************
double m_0x = 0.0;
double m_1x = 0.0;
double m_2x = 0.0;
double m_3x = 0.0;
double m_0y = 0.0;
double m_1y = 0.0;
double m_2y = 0.0;
double m_3y = 0.0;
double m_0z = 0.0;
double m_1z = 0.0;
double m_2z = 0.0;
double m_3z = 0.0;

//define r.h.s
double funct_x(double x, double y, double z) {
  return y;
}

double funct_y(double x, double y, double z) {
  return -a*y-b*x+c*sin(x)+d*sin(omega*z);
}

double funct_z(double x, double y, double z) {
  return 1.0;
}

main() {

  double y_alt = y;
  double x_min_alt = 0.0;
  double x_max_alt = 0.0;

  winsiz (800, 800);
  page (2800, 2800);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini ();
  pagera ();
  hwfont ();
  texmod ("on");


  // Phase diagramm
  axspos(250,1300);
  axslen(1000,1000);
  titlin ("Phase diagram", 4);
  graf(0.0,2.5,0.0,0.5,-3.0,3.0,-3.0,1.0);
  title();
  endgrf();

  // phi(t)
  axspos(1600,1300);
  axslen(1000,1000);
  titlin ("$\\phi(t)$", 4);
  graf(0.0,max_steps,0.0,max_steps/2.0,0.0,3.0,0.0,1.0);
  title();
  endgrf();


// Return map
  axspos(900,2600);
  axslen(1000,1000);
  titlin ("Return Map", 4);
  graf(2.0,2.3,2.0,0.1,2.0,2.3,2.0,0.1);
  title();
  endgrf();

  for ( step = 0; step < max_steps; step++ ) {
// classical RK4 method
    m_0x = h * funct_x(x, y, z);
    m_0y = h * funct_y(x, y, z);
    m_0z = h * funct_z(x, y, z);
    
    m_1x = h * funct_x(x + 0.5 * m_0x, y + 0.5 * m_0y, z + 0.5 * m_0z);
    m_1y = h * funct_y(x + 0.5 * m_0x, y + 0.5 * m_0y, z + 0.5 * m_0z);
    m_1z = h * funct_z(x + 0.5 * m_0x, y + 0.5 * m_0y, z + 0.5 * m_0z);

    m_2x = h * funct_x(x + 0.5 * m_1x, y + 0.5 * m_1y, z + 0.5 * m_1z);
    m_2y = h * funct_y(x + 0.5 * m_1x, y + 0.5 * m_1y, z + 0.5 * m_1z);
    m_2z = h * funct_z(x + 0.5 * m_1x, y + 0.5 * m_1y, z + 0.5 * m_1z);

    m_3x = h * funct_x(x + m_2x, y + m_2y, z + m_2z);
    m_3y = h * funct_y(x + m_2x, y + m_2y, z + m_2z);
    m_3z = h * funct_z(x + m_2x, y + m_2y, z + m_2z);

    y_alt = y;
    x += (m_0x + 2.0 * m_1x + 2.0 * m_2x + m_3x) / 6.0;
    y += (m_0y + 2.0 * m_1y + 2.0 * m_2y + m_3y) / 6.0;
    z += (m_0z + 2.0 * m_1z + 2.0 * m_2z + m_3z) / 6.0;
// plot point on the phase diagram (x,y)
    setrgb(0,0,0);
    axspos(250,1300);
    graf(0.0,2.5,0.0,0.5,-3.0,3.0,-3.0,1.0);
    setrgb(0,0,1);
    hsymbl(10);
    rlsymb(21,x,y);
   endgrf();
//----------- plot the oscillation curve x(t) ---------------------
   setrgb(0,0,0);
   axspos(1600,1300);
   graf(0.0,max_steps,0.0,max_steps/2.0,0.0,3.0,0.0,1.0);
   setrgb(0,0,1);
   hsymbl(10);
   rlsymb(21,step,x);
   endgrf();
//-------------find x_min or x_max and plot the return map x_n(x_n-1)-------
/*if ( (y_alt < 0.0) && (y > 0.0) ) {
	if (x_max_alt) {
          axspos(900,2600);
          graf(1.9,2.3,1.9,0.1,2.0,2.3,2.0,0.1);
	  rlsymb(21,x_max_alt,x);
	}
	x_max_alt = x;
      }
           endgrf();
*/
if ( (y_alt > 0.0) && (y < 0.0) ) {
	if (x_min_alt) {
          setrgb(0,0,0);
	  axspos(900,2600);
          graf(2.0,2.3,2.0,0.1,2.0,2.3,2.0,0.1);
          setrgb(0,0,1);
          hsymbl(12);
	  rlsymb(21,x_min_alt,x);
	}
	x_min_alt = x;
      }
  endgrf();
}
  disfin();
  
}