main() {
int n,N;
double d,V,U,shift;
// Scan grid size (n) and grid spacing (d)
printf("n: "); scanf("%d",&n);
printf("dx: "); scanf("%lf",&d);
// Compute total number of points in and volume of cell
N=n*n*n; V=pow(n*d,3);
printf("%d^3 grid of size %f^3=%lf (bohr^3)...\n",n,n*d,V);
// Initialize physics vectors on n^3 grid with spacing d in each direction
pvector ones(i3vec(n,n,n),d3vec(d,d,d));
pvector v1(i3vec(n,n,n),d3vec(d,d,d));
pvector v2(i3vec(n,n,n),d3vec(d,d,d));
pvector v3(i3vec(n,n,n),d3vec(d,d,d));
// Set pvector ones to all 1's.
ones.evalr(constant,1.);
// Does b_{G=0}(x_i)=ones?
v1.evalr(constant,0.); // Set v1=all zeros
v1.dat[0][0][0]=Complex(1.,0.); // Set G=[0][0][0] component to 1
v2=I(v1)-ones; // Compute difference
printf("b_{G=(0,0,0)}(x_i)=ones? %e\n",v2^v2); // Print magnitude of error
}