//************************************************************************ 
// Harmonical oscillator with the Heun's method
//  x'=y
//  y'=-x, x=x(t), y=y(t), T>t>0, x(0)=0, y(0)=1 
//************************************************************************
/*
compile with	
gcc -o test_heun test_heun.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/


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

#define NX 25 
#define NY 25
//
int step = 0;
double max_steps = 40.0*3.1415926; // final time
//
double x = 0.0; //initial conditions
double y = 1.0;
//
double h = 0.05; // time step
double L=6.0;
//
double m_0x = 0.0;
double m_1x = 0.0;

double m_0y = 0.0;
double m_1y = 0.0;
//
float xv[NX][NY], yv[NX][NY], xp[NX], yp[NY];

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

double funct_y(double x, double y) {
  return -x;
}

main ()
{ 

  int i, j, iret , iclr;

  double xstep, ystep;
// define grid for a vector plot
  xstep = 2*L / (NX - 1);; 
  ystep = 2*L / (NY - 1);; 
  for (i = 0;  i < NX; i++)
  { xp[i] = (float) (-L +  i * xstep);  
    for (j = 0; j < NY; j++) 
    { yp[j] = (float) (-L +  j *ystep); 
      xv[i][j] = yp[j];
      yv[i][j] = -xp[i];
    }
  }
  

  winsiz (600, 600);
  page (2600, 2600);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini ();
  pagera ();
  hwfont ();
  name ("x", "x");
  name ("y", "y");

  axspos (350, 2300);
  axslen (2000, 2000);
  titlin ("Phase diagram", 4);
  graf (-L, L, -L, 1., -L, L, -L, 1.); 

  height (50);
  title ();

 setrgb(0,0,1); // set blue color

 vecmat ((float *) xv, (float *) yv, NX, NY, xp, yp, -1); // plot vector field
 
//-----------------------------------------------------------------
// Heun's method
for ( step = 0; step < max_steps; step++ ) {

    m_0x = h * funct_x(x, y);
    m_0y = h * funct_y(x, y);
    
    m_1x = h * funct_x(x + m_0x, y + m_0y);
    m_1y = h * funct_y(x + m_0x, y + m_0y);

  
    x += (m_0x + m_1x) / 2.0;
    y += (m_0y + m_1y) / 2.0;

    setrgb(0,0,0);
    hsymbl(10);
    rlsymb(21,x,y);
   }
//-----------------------------------------------------------------

 /*stmmod ("rk4", "integration");
  stmmod ("on", "close");
  color ("blue");
  iclr = intrgb (0.0, 0.0, 0.0);
  vecclr (iclr);
  stream ((float *) xv, (float *) yv, NX, NY, xp, yp, NULL, NULL, 0); //plot the stream
*/
 
  disfin ();
}