# Help understanding this numerical surface integration technique?

I'm attempting to write a FORTRAN program that calculates the magnetic field, B, at any point outside of a bar magnet.

I'm going to use a first order euler scheme, where each side of the bar magnet is split into small cells, each with centres at (xi,yi,zi). I know I can ignore all of the sides that have any z values, and just focus on the top and bottom sides that are orientated in the x-y plane. So the method says this:

$\int f(x,y,z)dS = \Delta S \cdot \sum f(x_{i},y_{i}, z_{i})$

where the integral is over the surface S, and the summation is over i.

Delta S, each area, is given by $\hat{n}\cdot d\vec{S}$ , so if the cell is oriented in the x-y plane it's just $\Delta x\cdot \Delta y$ .

Here is a screenshot of the specific method instructions with a figure that demonstrates it

The function for the magnetic field is this: $\vec{B}(\vec{r}) = \frac{\mu _{0}}{4\pi }\cdot \int\frac{(r-{r}')\cdot M(r)\cdot \hat{n}}{|(r-{r}')|^{3}}$

Where the integral is over the surface S

I'm struggling to understand this method. I've tried to construct a flow chart, but can't get very far so I figured the problem is with the mathematics. Any help to understand it would be appreciated, and also any help with the flow chart would be fantastic. Here is my flowchart

Thanks