Matrix of Objects: Boundary conditions

[rolotomassi]

I have created a matrix with a class called Lattice. The lattice is filled with objects of type 'Dipole' which is created with another class. The problem I am having is with boundary conditions when I look for a neighbour. e.g If i pick a dipole on the top row, I want its 'above' neighbour to be the site on the bottom row etc...

Below is the .cpp and .h file which contains the functions im having problems with. Note that the value passed to both functions can be negative. I thought the %xSize should take care of this but still running into problems. I believe the dipole functions are working correctly.

Im running a Monte Carlo simulation by the way.

When it crashes it is always at an edge of the matrix. Usually the top row.
If my dipole is top row, it crashes when I find the 'above' neighbour.
When returning properties of this neighbour, "type" and "direction," defined in the dipole class to be 0,1 or 0,1,2,3 it returns large +ve or -ve values for the neighbour.

When I explicitly call GetDipole(a,b) and enter a and b to be negative values outside the size of the matrix, it always returns zero. If a and b are positive values outside the size of the matrix, it locates and returns the correct values.

PLEASE help as I cant see any problems :S

#ifndef lattice
#define lattice
#include "dipole.h"

class Lattice{

    private:

    Dipole* dipoleArray;

    double E;    // pointing in +ve y axis V/m

    public:
        const int xSize;
        const int ySize;
        Lattice(const int, const int);
        ~Lattice();
        Dipole* GetDipole(int x, int y);  //  return a dipole
        double calcLocalEnergy( int x, int y);


};



#endif


///////////////////////////////////////////////////////////


#include "dipole.h"
#include "lattice.h"
#include <cmath>
#include <cstdlib>
#include<iostream>
Lattice::Lattice(const int xSize, const int ySize): xSize(xSize), ySize(ySize){
    dipoleArray =  new Dipole[xSize*ySize];    //store 1D array list

    for(int i=0; i<xSize; i++){
        for(int j=0; j<ySize; j++){
            //loop through the entire array and create dipoles
            dipoleArray[ i+j*xSize ] = Dipole();  // i + j*xSize to mimic 2D lattice form
        
        }
    }
}

Dipole* Lattice::GetDipole(int x, int y){

    return &dipoleArray[ (x%xSize) + (y%ySize)*xSize  ];    // modulo is for loop behaviour
}
double Lattice::calcLocalEnergy( int x, int y){
    
    double k = 1/(4*3.14158*8.85E-12);     double p1=0;     double p2 =0;        double r  = 1E-8;     double energy = 0;
    double E =1e4;    double pZero = 1e-27;    double pOne = 1e-29;

//    std::cout<<"0"<<std::endl;
    Dipole* central = GetDipole(x,y);  // define central dipole in the lattice at position [x,y] given by the argument
                                                    // take modulo because x,y can potentially  = -1 OR xSize for neighbour dipoles
                                                    // modulo ensures neighbours are found correctly
                                                    // for 0 <=  x,y < xSize modulo has no effect
                                                    
    Dipole* neighbour[4];
    neighbour[0] = GetDipole(x,y-1);         //up = 0, right = 1, down = 2, left = 3
    neighbour[1] = GetDipole(x+1,y);
    neighbour[2] = GetDipole(x,y+1);
    neighbour[3] = GetDipole(x-1,y);
    int myDirection =  central->GetDirection() ;
    int type = central->GetType();
    if(type==0){p1 = pZero;    }    else{p1 = pOne;    }

    for(int i=0; i<4; i++){
    
        int neighbourDirection =  neighbour[i]->GetDirection();
        int Ntype = neighbour[i]->GetType();
        if(Ntype==1){    p2 = pOne;    }    else{    p2 = pZero;    }
    //    std::cout<<"1"<<std::endl;
        double    Energy[4][4]=
               {
                {-2*k*p1*p2/r/r/r - (p1+p2)*E, -p1*E, 2*k*p1*p2/r/r/r - (p1-p2)*E, -p1*E },    // ABOVE OR BELOW DIPOLES
                {-p2*E, 2*k*p1*p2/3/r/r/r, p2*E, -2*k*p1*p2/3/r/r/r    },
                {2*k*p1*p2/r/r/r + (p1-p2)*E, p1*E, -2*k*p1*p2/r/r/r + (p1+p2)*E, p1*E},
                {-p2*E, -2*k*p1*p2/3/r/r/r, p2*E, 2*k*p1*p2/3/r/r/r}
               };
        //       std::cout<<"2"<<std::endl;
       double    Energy2[4][4] =
               {
                { 2*k*p1*p2/3/r/r/r -(p1+p2)*E, -p1*E, -2*k*p1*p2/3/r/r/r -(p1-p2)*E, -p1*E },  // RIGHT OR LEFT ENERGIES
                { -p2*E, -2*k*p1*p2/r/r/r, p2*E, 2*k*p1*p2/r/r/r },
                { -2*k*p1*p2/3/r/r/r +(p1-p2)*E, p1*E, 2*k*p1*p2/3/r/r/r +(p1+p2)*E, p1*E },
                {-p2*E, 2*k*p1*p2/r/r/r, p2*E, -2*k*p1*p2/r/r/r }
               };
    //std::cout<<"3"<<std::endl;
    std::cout<<"i = "<< i << "\t" << type << "\t "<<Ntype <<"\t"<< neighbourDirection<<"\t"<<myDirection<<std::endl;
    std::cout<<"x = " << x << "y = "<< y <<std::endl;
            if(i==0 || i==2){
                energy += Energy[myDirection][neighbourDirection];
            }    
        
            else if(i==1 || i==3){         
                energy += Energy2[myDirection][neighbourDirection];        
            }    
    }//std::cout<<"4"<<std::endl;
    std::cout<<"4end clac local E"<<std::endl;
     return energy;
  
}

Lattice::~Lattice(){

    delete[] dipoleArray;
}