# include <cmath>
# include <cstdlib>
# include <ctime>
# include <fstream>
# include <iomanip>
# include <iostream>

using namespace std;

int main ( );
int i4_modp ( int i, int j );
int i4_wrap ( int ival, int ilo, int ihi );
double *lorenz96_rhs ( double t, int m, double x[] );
double r8_normal_01 ( );
double *r8vec_linspace_new ( int n, double a, double b );
double *rk4vec ( double t0, int n, double u0[], double dt, 
  double *f ( double t, int n, double u[] ) );
void timestamp ( );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    lorenz96_ode() solves the Lorenz96 ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    30 September 2025
//
//  Author:
//
//    John Burkardt
//
{
  string command_filename = "lorenz96_ode_commands.txt";
  ofstream command_unit;
  string data_filename = "lorenz96_ode_data.txt";
  ofstream data_unit;
  double dt;
  double force = 8.0;
  int i;
  int j;
  int n = 4;
  int nt = 2000;
  double perturb = 0.001;
  double *t;
  double t0 = 0.0;
  double tstop = 30.0;
  double *y;
  double *ynew;

  timestamp ( );
  cout << "\n";
  cout << "lorenz96_ode():\n";
  cout << "  C++ version\n";
  cout << "  Compute solutions of the Lorenz96 system.\n";
  cout << "  Write data to a file for use by gnuplot.\n";
//
//  Initialize the random number generator.
//
  srand48 ( time ( NULL ) );
//
//  Data
//
  dt = ( tstop - t0 ) / ( double ) ( nt );
//
//  Store the initial conditions in entry 0.
//
  t = r8vec_linspace_new ( nt + 1, t0, tstop );

  y = new double[n*(nt+1)];
  for ( i = 0; i < n; i++ )
  {
    y[i+0*nt] = force + perturb * r8_normal_01 ( );
  }
//
//  Compute the approximate solution at equally spaced times.
//
  for ( j = 0; j < nt; j++ )
  {
    ynew = rk4vec ( t[j], n, y+j*n, dt, lorenz96_rhs );
    for ( i = 0; i < n; i++ )
    {
      y[i+(j+1)*n] = ynew[i];
    }
    delete [] ynew; 
  }
//
//  Create the plot data file.
//
  data_unit.open ( data_filename.c_str ( ) );
  for ( j = 0; j <= nt; j++ )
  {
    data_unit << "  " << setw(14) << t[j]
              << "  " << setw(14) << y[0+j*n]
              << "  " << setw(14) << y[1+j*n]
              << "  " << setw(14) << y[2+j*n]
              << "  " << setw(14) << y[3+j*n] << "\n";
  }
  data_unit.close ( );
  cout << "  Created data file \"" << data_filename << "\".\n";
/*
  Create the plot command file.
*/
  command_unit.open ( command_filename.c_str ( ) );
  command_unit << "# " << command_filename << "\n";
  command_unit << "#\n";
  command_unit << "# Usage:\n";
  command_unit << "#  gnuplot < " << command_filename << "\n";
  command_unit << "#\n";
  command_unit << "set term png\n";
  command_unit << "set output 'lorenz96_time.png'\n";
  command_unit << "set xlabel '<--- T --->'\n";
  command_unit << "set ylabel '<--- Y0(T), Y1(T), Y2(T) --->'\n";
  command_unit << "set title 'Y0, Y1, Y2 versus Time'\n";
  command_unit << "set grid\n";
  command_unit << "set style data lines\n";
  command_unit << "plot '" << data_filename 
               << "' using 1:2 lw 3 linecolor rgb 'blue',";
  command_unit << "'' using 1:3 lw 3 linecolor rgb 'red',";
  command_unit << "'' using 1:4 lw 3 linecolor rgb 'green'\n";
  command_unit << "set output 'lorenz96_3d.png'\n";
  command_unit << "set xlabel '<--- Y0(T) --->'\n";
  command_unit << "set ylabel '<--- Y1(T) --->'\n";
  command_unit << "set zlabel '<--- Y2(T) --->'\n";
  command_unit << "set title '(Y0(T),Y1(T),Y2(T)) trajectory'\n";
  command_unit << "set grid\n";
  command_unit << "set style data lines\n";
  command_unit << "splot '" << data_filename 
               << "' using 2:3:4 lw 1 linecolor rgb 'blue'\n";
  command_unit << "quit\n";
  command_unit.close ( );
  cout << "  Created command file '" << command_filename << "'\n";
//
//  Terminate.
//
  cout << "\n";
  cout << "lorenz96_ode():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

int i4_modp ( int i, int j )

//****************************************************************************80
//
//  Purpose:
//
//    i4_modp() returns the nonnegative remainder of I4 division.
//
//  Discussion:
//
//    If
//      NREM = I4_MODP ( I, J )
//      NMULT = ( I - NREM ) / J
//    then
//      I = J * NMULT + NREM
//    where NREM is always nonnegative.
//
//    The MOD function computes a result with the same sign as the
//    quantity being divided.  Thus, suppose you had an angle A,
//    and you wanted to ensure that it was between 0 and 360.
//    Then mod(A,360) would do, if A was positive, but if A
//    was negative, your result would be between -360 and 0.
//
//    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
//
//        I         J     MOD  I4_MODP   I4_MODP Factorization
//
//      107        50       7       7    107 =  2 *  50 + 7
//      107       -50       7       7    107 = -2 * -50 + 7
//     -107        50      -7      43   -107 = -3 *  50 + 43
//     -107       -50      -7      43   -107 =  3 * -50 + 43
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 May 1999
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int I, the number to be divided.
//
//    int J, the number that divides I.
//
//  Output:
//
//    int I4_MODP, the nonnegative remainder when I is divided by J.
//
{
  int value;

  if ( j == 0 )
  {
    cerr << "\n";
    cerr << "I4_MODP(): Fatal error!\n";
    cerr << "  I4_MODP ( I, J ) called with J = " << j << "\n";
    exit ( 1 );
  }

  value = i % j;

  if ( value < 0 )
  {
    value = value + abs ( j );
  }

  return value;
}
//****************************************************************************80

int i4_wrap ( int ival, int ilo, int ihi )

//****************************************************************************80
//
//  Purpose:
//
//    i4_wrap() forces an I4 to lie between given limits by wrapping.
//
//  Example:
//
//    ILO = 4, IHI = 8
//
//    I   Value
//
//    -2     8
//    -1     4
//     0     5
//     1     6
//     2     7
//     3     8
//     4     4
//     5     5
//     6     6
//     7     7
//     8     8
//     9     4
//    10     5
//    11     6
//    12     7
//    13     8
//    14     4
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 December 2012
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int IVAL, an integer value.
//
//    int ILO, IHI, the desired bounds for the integer value.
//
//  Output:
//
//    int I4_WRAP, a "wrapped" version of IVAL.
//
{
  int jhi;
  int jlo;
  int value;
  int wide;

  if ( ilo <= ihi )
  {
    jlo = ilo;
    jhi = ihi;
  }
  else
  {
    jlo = ihi;
    jhi = ilo;
  }

  wide = jhi + 1 - jlo;

  if ( wide == 1 )
  {
    value = jlo;
  }
  else
  {
    value = jlo + i4_modp ( ival - jlo, wide );
  }

  return value;
}
//****************************************************************************80

double *lorenz96_rhs ( double t, int n, double y[] )

//****************************************************************************80
//
//  Purpose:
//
//    lorenz96_rhs() evaluates the right hand side of the Lorenz ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    30 September 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double t, the value of the independent variable.
//
//    int n, the spatial dimension.
//
//    double y[n], the values of the dependent variables
//    at time T.
//
//  Output:
//
//    double dydt[n], the values of the derivatives
//    of the dependent variables at time T.
//
{
  double *dydt;
  double force = 8.0;
  int i;
  int im1;
  int im2;
  int ip1;

  dydt = new double[n];

  for ( i = 0; i < n; i++ )
  {
    ip1 = i4_wrap ( i + 1, 0, n - 1 );
    im1 = i4_wrap ( i - 1, 0, n - 1 );
    im2 = i4_wrap ( i - 2, 0, n - 1 );
    printf ( "im2, im1, i, ip1: %d %d %d %d\n", im2, im1, i, ip1 );
    dydt[i] = ( y[ip1] - y[im2] ) * y[im1] - y[i] + force;
  }

  return dydt;
}
//****************************************************************************80

double r8_normal_01 ( )

//****************************************************************************80
//
//  Purpose:
//
//    r8_normal_01() returns a unit pseudonormal R8.
//
//  Discussion:
//
//    The standard normal probability distribution function (PDF) has
//    mean 0 and standard deviation 1.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 September 2022
//
//  Author:
//
//    John Burkardt
//
//  Output:
//
//    double R8_NORMAL_01, a normally distributed random value.
//
{
  double r1;
  double r2;
  const double r8_pi = 3.141592653589793;
  double x;

  r1 = drand48 ( );
  r2 = drand48 ( );
  x = sqrt ( - 2.0 * log ( r1 ) ) * cos ( 2.0 * r8_pi * r2 );

  return x;
}
//****************************************************************************80

double *r8vec_linspace_new ( int n, double a_first, double a_last )

//****************************************************************************80
//
//  Purpose:
//
//    r8vec_linspace_new() creates a vector of linearly spaced values.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//    4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
//
//    In other words, the interval is divided into N-1 even subintervals,
//    and the endpoints of intervals are used as the points.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2011
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of entries in the vector.
//
//    double A_FIRST, A_LAST, the first and last entries.
//
//  Output:
//
//    double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data.
//
{
  double *a;
  int i;

  a = new double[n];

  if ( n == 1 )
  {
    a[0] = ( a_first + a_last ) / 2.0;
  }
  else
  {
    for ( i = 0; i < n; i++ )
    {
      a[i] = ( ( double ) ( n - 1 - i ) * a_first 
             + ( double ) (         i ) * a_last ) 
             / ( double ) ( n - 1     );
    }
  }
  return a;
}
//****************************************************************************80

double *rk4vec ( double t0, int m, double u0[], double dt, 
  double *f ( double t, int m, double u[] ) )

//****************************************************************************80
//
//  Purpose:
//
//    rk4vec() takes one Runge-Kutta step for a vector ODE.
//
//  Discussion:
//
//    It is assumed that an initial value problem, of the form
//
//      du/dt = f ( t, u )
//      u(t0) = u0
//
//    is being solved.
//
//    If the user can supply current values of t, u, a stepsize dt, and a
//    function to evaluate the derivative, this function can compute the
//    fourth-order Runge Kutta estimate to the solution at time t+dt.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    09 October 2013
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double T0, the current time.
//
//    int M, the spatial dimension.
//
//    double U0[M], the solution estimate at the current time.
//
//    double DT, the time step.
//
//    double *F ( double T, int M, double U[] ), a function which evaluates
//    the derivative, or right hand side of the problem.
//
//  Output:
//
//    double RK4VEC[M], the fourth-order Runge-Kutta solution estimate
//    at time T0+DT.
//
{
  double *f0;
  double *f1;
  double *f2;
  double *f3;
  int i;
  double t1;
  double t2;
  double t3;
  double *u;
  double *u1;
  double *u2;
  double *u3;
//
//  Get four sample values of the derivative.
//
  f0 = f ( t0, m, u0 );

  t1 = t0 + dt / 2.0;
  u1 = new double[m];
  for ( i = 0; i < m; i++ )
  {
    u1[i] = u0[i] + dt * f0[i] / 2.0;
  }
  f1 = f ( t1, m, u1 );

  t2 = t0 + dt / 2.0;
  u2 = new double[m];
  for ( i = 0; i < m; i++ )
  {
    u2[i] = u0[i] + dt * f1[i] / 2.0;
  }
  f2 = f ( t2, m, u2 );

  t3 = t0 + dt;
  u3 = new double[m];
  for ( i = 0; i < m; i++ )
  {
     u3[i] = u0[i] + dt * f2[i];
  }
  f3 = f ( t3, m, u3 );
//
//  Combine them to estimate the solution.
//
  u = new double[m];
  for ( i = 0; i < m; i++ )
  {
     u[i] = u0[i] + dt * ( f0[i] + 2.0 * f1[i] + 2.0 * f2[i] + f3[i] ) / 6.0;
  }
//
//  Free memory.
//
  delete [] f0;
  delete [] f1;
  delete [] f2;
  delete [] f3;
  delete [] u1;
  delete [] u2;
  delete [] u3;

  return u;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}