Creating the Electric Octupole Tensor of a cubic electric octupole

[PhDeezNutz]

TL;DR Summary

Hello, I want to create an electric Octupole Tensor of the charge distribution below. I set the side length of the cube equal to 1 and the charge equal to fundamental charge (because why not). Here are the details. My approach is to use four nested for loops. 3 for the coordinates. 1 for the charge. With the charge for loop being the innermost for loop.

I created an array, where the first three entries of each column are the x,y, and z coordinates. The last entry in each column is the charge. I called this array PCQ.

l/2	l/2	-l/2	-l/2	-l/2	l/2	l/2	-l/2
-l/2	l/2	l/2	-l/2	-l/2	-l/2	l/2	l/2
l/2	l/2	l/2	l/2	-l/2	-l/2	-l/2	-l/2
q	-q	q	-q	q	-q	q	-q

The Electric Quadrupole Tensor is defined in index notation as follows, where $r_{i,\ell}, r_{j,\ell},r_{k'ell}$ could be anyone of the cartesian coordinates of the $\ell$'th charge.

$Q_{ijk} = \sum_{\ell}^{8} \sum_{i,j,k}^{3} q_{\ell} r_{i ,\ell} r_{j ,\ell} r_{k ,\ell}$

My code is

%1) Define the grid. n = 100;rmax = 1200;x = linspace(-rmax,rmax,n);y = linspace(-rmax,rmax,n);z = linspace(-rmax,rmax,n); [X,Y,Z] = meshgrid(x,y,z); %2) Need to define constants c = 3*10^(8);epsilon = 8.85*(10^(-12));mu = (4*pi)*((10)^(-7)); q = 1.602*(10^(-19)); %fundamental chargel = 1; %side length of cubel2 = 0.5*l; %half side length of cubek = 0.001*5*pi;omega = k*sqrt(mu*epsilon); %3) Create Primed Coordinates with charge pcx = [l2,l2,-l2,-l2,-l2,l2,l2,-l2];pcy = [-l2,l2,l2,-l2,-l2,-l2,l2,l2];pcz = [l2,l2,l2,l2,-l2,-l2,-l2,-l2];pcq = [q,-q,q,-q,q,-q,q,-q]; PCQ = vertcat(pcx,pcy,pcz,pcq); %4) Create the Electric Octupole Tensor Q = zeros(3,3,3); for i = 1:3 for j = 1:3 for k = 1:3 Q(i,j,k) = 0; for i1 = 1:length(PCQ) Q(i,j,k) = Q(i,j,k) + Q(i,j,k) + PCQ(4,i1)*PCQ(i,i1)*PCQ(i,i1)*PCQ(j,i1)*PCQ(k,i1); end end endend

My output is

I have no idea if this is right. And calculating these moments by hand would simply be too difficult. I am unnerved somewhat by each "column" of each "page" having the same entry. I am hopeful that someone with more experience, knowledge, and intuition could take a brief look at my code and tell me if it is correct or if it needs to be corrected (and how to correct it)

Thanks in advance, for some reason the latex is not rendering.

 

[PhDeezNutz]

Edit: I spotted an error in my code. I had

Q(i,j,k) = Q(i,j,k) + Q(i,j,k) + PC(4,i1)*PC(i,i1)*PC(j,i1)*PC(k,i1)

It should be

Q(i,j,k) = Q(i,j,k) + PC(4,i1)*PC(i,i1)*PC(j,i1)*PC(k,i1)

My output is

I still have no idea if this is right and what sanity checks I should do.

I would think that some of the entries should be positive.