- #1

- 53

- 0

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.

Code:

```
#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;
}
```

EDIT:: 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

PLEASE help as I cant see any problems :S