trottersuzuki.h

#ifndef __TROTTERSUZUKI_H
#define __TROTTERSUZUKI_H

#include <string>
#define _USE_MATH_DEFINES
#include <cfloat>
#include <complex>
#if HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_MPI
#include <mpi.h>
#endif

using namespace std;

class Lattice {
public:
    int mpi_rank;    
    int mpi_procs;    
    double length_x, length_y;    
    double delta_x, delta_y;    
    int dim_x, dim_y;    
    int global_no_halo_dim_x, global_no_halo_dim_y;    
    int global_dim_x, global_dim_y;    
    int periods[2];    
    string coordinate_system;   

    // Computational topology
    int halo_x, halo_y;    
    int start_x, start_y;    
    int end_x, end_y;    
    int inner_start_x, inner_start_y;    
    int inner_end_x, inner_end_y;    
    int mpi_coords[2], mpi_dims[2];    
#ifdef HAVE_MPI
    MPI_Comm cartcomm;    
#endif
};

class Lattice1D: public Lattice {
public:
    Lattice1D(int dim, double length, bool periodic_x_axis = false, string coordinate_system = "cartesian");
};

class Lattice2D: public Lattice {
public:
    Lattice2D(int dim, double length,
              bool periodic_x_axis = false, bool periodic_y_axis = false,
              double angular_velocity = 0., string coordinate_system = "cartesian");
    Lattice2D(int dim_x, double length_x, int dim_y, double length_y,
              bool periodic_x_axis = false, bool periodic_y_axis = false,
              double angular_velocity = 0., string coordinate_system = "cartesian");
private:
    void init(int dim_x, double length_x, int dim_y, double length_y,
              bool periodic_x_axis = false, bool periodic_y_axis = false,
              double angular_velocity = 0., string coordinate_system = "cartesian");
};

class State {
public:
    double *p_real;    
    double *p_imag;    
    Lattice *grid;    
    int angular_momentum;   

    State(Lattice *grid, int angular_momentum = 0, double *p_real = 0, double *p_imag = 0);
    State(const State &obj );    
    ~State();    
    void init_state(complex<double> (*ini_state)(double x) ); 
    void init_state(complex<double> (*ini_state)(double x, double y) );    
    void loadtxt(char *file_name);    

    void imprint(complex<double> (*function)(double x) );    
    void imprint(complex<double> (*function)(double x, double y) );    
    double *get_particle_density(double *density = 0 );  
    double *get_phase(double *phase = 0 );  
    double get_expected_value(string _operator);    
    double get_squared_norm(void);    
    double get_mean_x(void);    
    double get_mean_xx(void);    
    double get_mean_y(void);    
    double get_mean_yy(void);    
    double get_mean_px(void);    
    double get_mean_pxpx(void);    
    double get_mean_py(void);    
    double get_mean_pypy(void);    
    double get_mean_angular_momentum(void);    

    void write_to_file(string fileprefix );    
    void write_particle_density(string fileprefix );    
    void write_phase(string fileprefix );    
    bool expected_values_updated;    

protected:
    bool self_init;    
    void calculate_expected_values(void);    
    double mean_X, mean_XX;    
    double mean_Y, mean_YY;    
    double mean_Px, mean_PxPx;    
    double mean_Py, mean_PyPy;    
    double mean_angular_momentum;    
    double norm2;    
};

class BesselState: public State {
public:
    BesselState(Lattice1D *grid, int angular_momentum = 0, int zeros = 1, double norm = 1, double phase = 0, double *p_real = 0, double *p_imag = 0);

    BesselState(Lattice2D *grid, int angular_momentum = 0, int zeros = 1, int n_y = 0, double norm = 1, double phase = 0, double *p_real = 0, double *p_imag = 0);

private:
    int angular_momentum;   
    int zeros, n_y;    
    double zero;    
    double normalization;    
    double norm, phase;    
    complex<double> bessel_state1D(double x);    
    complex<double> bessel_state2D(double x, double y);    
};

class ExponentialState: public State {
public:
    ExponentialState(Lattice1D *grid, int n_x = 1, double norm = 1, double phase = 0, double *p_real = 0, double *p_imag = 0);

    ExponentialState(Lattice2D *grid, int n_x = 1, int n_y = 1, double norm = 1, double phase = 0, double *p_real = 0, double *p_imag = 0);

private:
    int n_x, n_y;    
    double norm, phase;    
    complex<double> exp_state(double x, double y);    
};

class GaussianState: public State {
public:
    GaussianState(Lattice1D *grid, double omega_x, double mean_x = 0, double norm = 1, double phase = 0,
                  double *p_real = 0, double *p_imag = 0);

    GaussianState(Lattice2D *grid, double omega_x, double omega_y = -1., double mean_x = 0, double mean_y = 0, double norm = 1, double phase = 0,
                  double *p_real = 0, double *p_imag = 0);

private:
    double mean_x;    
    double mean_y;    
    double omega_x;    
    double omega_y;    
    double norm;    
    double phase;    
    complex<double> gauss_state(double x, double y);    
};

class SinusoidState: public State {
public:
    SinusoidState(Lattice1D *grid, int n_x = 1, double norm = 1, double phase = 0, double *p_real = 0, double *p_imag = 0);

    SinusoidState(Lattice2D *grid, int n_x = 1, int n_y = 1, double norm = 1, double phase = 0, double *p_real = 0, double *p_imag = 0);

private:
    int n_x, n_y;    
    double norm, phase;    
    complex<double> sinusoid_state(double x, double y);    
};

class Potential {
public:
    Lattice *grid;    
    double *matrix;    

    Potential(Lattice *grid, char *filename);
    Potential(Lattice *grid, double *external_pot = 0);
    Potential(Lattice *grid, double (*potential_function)(double x, double y));
    Potential(Lattice *grid, double (*potential_function)(double x, double y, double t), int t = 0);
    virtual ~Potential();
    virtual double get_value(int x); 
    virtual double get_value(int x, int y);    
    bool update(double t);    
    bool updated_potential_matrix;
protected:
    double current_evolution_time;    
    double (*static_potential)(double x, double y);    
    double (*evolving_potential)(double x, double y, double t);    
    bool self_init;    
    bool is_static;    
};

class HarmonicPotential: public Potential {
public:
    HarmonicPotential(Lattice2D *grid, double omegax, double omegay, double mass = 1., double mean_x = 0., double mean_y = 0.);
    ~HarmonicPotential();
    double get_value(int x, int y);    

private:
    double omegax, omegay;    
    double mass;    
    double mean_x, mean_y;    
};

class Hamiltonian {
public:
    Potential *potential;   
    double mass;    
    double coupling_a;    
    double LeeHuangYang_coupling_a;    
    double angular_velocity;    
    double rot_coord_x;    
    double rot_coord_y;    
    double azimuthal_potential(double x, int angular_momentum);

    Hamiltonian(Lattice *grid, Potential *potential = 0, double mass = 1., double coupling_a = 0., double LeeHuangYang_coupling_a = 0.,
                double angular_velocity = 0.,
                double rot_coord_x = 0, double rot_coord_y = 0);
    ~Hamiltonian();

protected:
    bool self_init;    
    Lattice *grid;    
};

class Hamiltonian2Component: public Hamiltonian {
public:
    double mass_b;    
    double coupling_ab;    
    double coupling_b;    
    //double LeeHuangYang_coupling_b;    ///< Coupling constant of the Lee-Huang-Yang term of the second component.
    double omega_r;    
    double omega_i;    
    Potential *potential_b;    
    double azimuthal_potential_b(double x, int angular_momentum);

    Hamiltonian2Component(Lattice *grid, Potential *potential = 0,
                          Potential *potential_b = 0,
                          double mass = 1., double mass_b = 1.,
                          double coupling_a = 0.,/* double LeeHuangYang_coupling_a = 0.,*/ double coupling_ab = 0.,
                          double coupling_b = 0.,/* double LeeHuangYang_coupling_b = 0.,*/
                          double omega_r = 0, double omega_i = 0,
                          double angular_velocity = 0.,
                          double rot_coord_x = 0,
                          double rot_coord_y = 0);
    ~Hamiltonian2Component();
};

class ITrotterKernel {
public:
    virtual ~ITrotterKernel() {};
    virtual void run_kernel() = 0;    
    virtual void run_kernel_on_halo() = 0;    
    virtual void wait_for_completion() = 0;    
    virtual void get_sample(size_t dest_stride, size_t x, size_t y, size_t width, size_t height, double * dest_real, double * dest_imag, double * dest_real2 = 0, double * dest_imag2 = 0) const = 0; 
    virtual void normalization() = 0;    
    virtual void rabi_coupling(double var, double delta_t) = 0;    
    virtual double calculate_squared_norm(bool global = true) const = 0;  
    virtual bool runs_in_place() const = 0;
    virtual string get_name() const = 0;                
    virtual void update_potential(double *_external_pot_real, double *_external_pot_imag, int which) = 0;    
    virtual void cpy_first_positive_to_first_negative() = 0;    

    virtual void start_halo_exchange() = 0;                 
    virtual void finish_halo_exchange() = 0;                

};

class Solver {
public:
    Lattice *grid;    
    State *state;    
    State *state_b;    
    Hamiltonian *hamiltonian;    
    double current_evolution_time;    

    Solver(Lattice *grid, State *state, Hamiltonian *hamiltonian, double delta_t,
           string kernel_type = "cpu");
    Solver(Lattice *grid, State *state1, State *state2,
           Hamiltonian2Component *hamiltonian,
           double delta_t, string kernel_type = "cpu");
    ~Solver();
    void evolve(int iterations, bool imag_time = false);  
    void update_parameters();  
    double get_total_energy(void);    
    double get_squared_norm(size_t which = 3 );  
    double get_kinetic_energy(size_t which = 3 );  
    double get_potential_energy(size_t which = 3 );  
    double get_rotational_energy(size_t which = 3 );  
    double get_intra_species_energy(size_t which = 3 );  
    double get_LeeHuangYang_energy(void);    
    double get_inter_species_energy(void);    
    double get_rabi_energy(void);    
    void set_exp_potential(double *real, int real_length, double *imag,
                           int imag_length, int which); 
private:
    bool imag_time;    
    double **external_pot_real;    
    double **external_pot_imag;    
    double delta_t;    
    double norm2[2];    
    bool single_component;    
    string kernel_type;    
    ITrotterKernel * kernel;    
    void initialize_exp_potential(double time_single_it, int which);    
    void init_kernel();    
    double total_energy;    
    double kinetic_energy[2];    
    double tot_kinetic_energy;    
    double potential_energy[2];    
    double tot_potential_energy;    
    double rotational_energy[2];    
    double tot_rotational_energy;    
    double intra_species_energy[2];    
    double tot_intra_species_energy;    
    double LeeHuangYang_energy;    
    double inter_species_energy;    
    double rabi_energy;    
    bool has_parameters_changed;   
    bool energy_expected_values_updated;    
    void calculate_energy_expected_values(void);    
    bool is_python;
};

double const_potential(double x);    
double const_potential(double x, double y);    
void map_lattice_to_coordinate_space(Lattice *grid, int x_in, double *x_out);  
void map_lattice_to_coordinate_space(Lattice *grid, int x_in, int y_in, double *x_out, double *y_out); 
#endif // __TROTTERSUZUKI_H