# Step P1: Import needed NRPy+ core modules:
from outputC import lhrh         # NRPy+: Core C code output module
import NRPy_param_funcs as par   # NRPy+: Parameter interface
import finite_difference as fin  # NRPy+: Finite difference C code generation module
import grid as gri               # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support

# Set the spatial dimension to 1
par.set_paramsvals_value("grid::DIM = 1")

# Register the input gridfunction "phi" and the gridfunction to which data are output, "output":
phi, output = gri.register_gridfunctions("AUX",["phi","output"])

# Declare phi_dDD as a rank-2 indexed expression: phi_dDD[i][j] = \partial_i \partial_j phi
phi_dDD = ixp.declarerank2("phi_dDD","nosym")

# Set output to \partial_0^2 phi
output = phi_dDD[0][0]

# Output to the screen the core C code for evaluating the finite difference derivative
fin.FD_outputC("stdout",lhrh(lhs=gri.gfaccess("out_gf","output"),rhs=output))

# Now, output the above C code to a file named "finite_diff_tutorial-second_deriv.h".
fin.FD_outputC("finite_diff_tutorial-second_deriv.h",lhrh(lhs=gri.gfaccess("aux_gfs","output"),rhs=output))

%%writefile finite_difference_playground.c

// Part P1: Import needed header files
#include "stdio.h"  // Provides printf()
#include "stdlib.h" // Provides malloc() and free()
#include "math.h"   // Provides sin()

// Part P2: Declare the IDX2S(gf,i) macro, which enables us to store 2-dimensions of
//          data in a 1D array. In this case, consecutive values of "i"
//          ("gf" held to a fixed value) are consecutive in memory, where
//          consecutive values of "gf" (fixing "i") are separated by N elements in
//          memory.
#define IDX2S(gf, i) ( (i) + Npts_in_stencil * (gf) )

// Part P3: Set PHIGF and OUTPUTGF macros
#define PHIGF    0
#define OUTPUTGF 1

// Part P4: Import code generated by NRPy+ to compute f''(x)
//          as a finite difference derivative.
void f_dDD_FD(double *in_gfs,double *aux_gfs,const int i0,const int Npts_in_stencil,const double invdx0)  {
#include "finite_diff_tutorial-second_deriv.h"
}

// Part P5: Define the function we wish to differentiate, as well as its exact second derivative:
double f(const double x)           { return  sin(x); } // f(x)
double f_dDD_exact(const double x) { return -sin(x); } // f''(x)

// Part P6: Define x_i = (x_0 + i*Delta_x)
double x_i(const double x_0,const int i,const double Delta_x) {
    return (x_0 + (double)i*Delta_x);
}


// main() function
int main(int argc,char *argv[]) {
    // Step 0: Read command-line arguments (TODO)

    // Step 1: Set some needed constants
    const int Npts_in_stencil = 5; // Equal to the finite difference order, plus one.  '+str(par.parval_from_str("finite_difference::FD_CENTDERIVS_ORDER"))+'+1;
    const double PI = 3.14159265358979323846264338327950288; // The scale over which the sine function varies.
    const double x_eval = PI/4.0; // x_0 = desired x at which we wish to compute f(x)

    // Step 2: Evaluate f''(x_eval) using the exact expression:
    double EX = f_dDD_exact(x_eval);

    // Step 3: Allocate space for two gridfunctions
    double *in_gfs = (double *)malloc(sizeof(double)*Npts_in_stencil*2);

    // Step 4: Loop over grid spacings
    for(double Delta_x = 1e-3*(2*PI);Delta_x<=1.5e-1*(2*PI);Delta_x*=1.1) {

        // Step 4a: x_eval is the center point of the finite differencing stencil,
        //          thus x_0 = x_eval - 2*dx for fourth-order-accurate first & second finite difference derivs,
        //          and  x_0 = x_eval - 3*dx for sixth-order-accurate first & second finite difference derivs, etc.
        //          In general, for the integer Npts_in_stencil, we have
        //          x_0 = x_eval - (double)(Npts_in_stencil/2)*Delta_x,
        //          where we rely upon integer arithmetic (which always rounds down) to ensure
        //          Npts_in_stencil/2 = 5/2 = 2 for fourth-order-accurate first & second finite difference derivs:
        const double x_0 = x_eval - (double)(Npts_in_stencil/2)*Delta_x;

        // Step 4b: Set \phi=PHIGF to be f(x) as defined in the
        //          f(const double x) function above, where x_i = stencil_start_x + i*Delta_x:
        for(int ii=0;ii<Npts_in_stencil;ii++) {
            in_gfs[IDX2S(PHIGF, ii)] = f(x_i(x_0,ii,Delta_x));
        }

        // Step 4c: Set invdx0, which is needed by the NRPy+-generated "finite_diff_tutorial-second_deriv.h"
        const double invdx0 = 1.0/Delta_x;

        // Step 4d: Evaluate the finite-difference second derivative of f(x):
        const int i0 = Npts_in_stencil/2; // The derivative is evaluated at the center of the stencil.
        f_dDD_FD(in_gfs,in_gfs,i0,Npts_in_stencil,invdx0);
        double FD = in_gfs[IDX2S(OUTPUTGF,i0)];

        // Step 4e: Print log_10(\Delta x) and log_10([relative error])
        printf("%e\t%.15e\n",log10(Delta_x),log10(fabs((EX-FD)/(EX))));
    }

    // Step 5: Free the allocated memory for the gridfunctions.
    free(in_gfs);

    return 0;
}

import cmdline_helper as cmd
cmd.C_compile("finite_difference_playground.c", "fdp")
cmd.delete_existing_files("data.txt")
cmd.Execute("fdp", "", "data.txt")

%matplotlib inline
import matplotlib.pyplot as plt

# from https://stackoverflow.com/questions/52386747/matplotlib-check-for-empty-plot
import numpy
x, y = numpy.loadtxt(fname='data.txt', delimiter='\t', unpack=True)
fig, ax = plt.subplots()
ax.set_xlabel('log10(Delta_x)')
ax.set_ylabel('log10([Relative Error])')
ax.plot(x, y, color='b')

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish-Finite_Difference_Playground")