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

double Z_x = 0.5;
double Z_y = 2.0;
double A_x = 3.0;
double A_y = .5;

double scale_factor = 0.05;

void func(double x, double y, double *field_x, double *field_y) {
  
  *field_x = Z_x * x * y - A_x * x;
  *field_y = Z_y * y - A_y * x * y ;

}

void euler_step(double delta_t, double *x, double *y, void (*field)(double, double, double *, double *)) {
  
  double delta_x = 0.0;
  double delta_y = 0.0;

  field(*x, *y, &delta_x, &delta_y);

  *x += delta_t * delta_x;
  *y += delta_t * delta_y;

}

void plot_field(double x_0, double x_1, double y_0, double y_1, int gridpoints, void (*field)(double, double, double *, double *)) {

  int i = 0;
  int j = 0;
  
  double x_step = (x_1 - x_0) / gridpoints;
  double y_step = (y_1 - y_0) / gridpoints;

  double f_x = 0.0;
  double f_y = 0.0;

  for ( i = 0; i < gridpoints; i++) {
    for ( j = 0; j < gridpoints; j++) {
      
      field(x_0 + i * x_step, y_0 + j * y_step, &f_x, &f_y);

      rlvec(x_0 + i * x_step, y_0 + j * y_step, x_0 + i * x_step +  scale_factor * f_x, y_0 + j * y_step + scale_factor * f_y, 1001);
    }
  }
  
}

int main () {
  
  double x_0 = -20.0;
  double x_1 = 20.0;
  double y_0 = -20.0;
  double y_1 = 20.0;

  int step = 0;
  int max_step = 10000;

  double delta_t = 0.001;

  double x = 12.0;
  double y = 10.0;

  metafl("CONS");
  x11mod("STORE");

  disini();
  
  graf(x_0, x_1, x_0, (x_1 - x_0) / 10.0, y_0, y_1, y_0, (y_1 - y_0) / 10.0 );
  plot_field(x_0, x_1, y_0, y_1, 20, &func);

  for ( step; step < max_step; step++ ) {

    euler_step(delta_t, &x, &y, &func);
    rlsymb(21,x, y);
    sendbf();
    printf("Step: %d, (%f,%f)\n",step, x,y);
  }
  
  disfin();
  
  return 0;
}