Intensity of the diffraction pattern of the double slit

Question

I am trying another approach for my last unanswered question. (Bounty still on for 3 days. Anyone? Please?) Note that this is not the same question but a greatly simplified version concerning a much more general situation. As this may also be viewed as a problem of probability distributions, I am also going for the quantum tags Many textbooks talk about the "brightness" distribution of the double slit. The one here (Lang/Pucker, Mathematische Methoden der Physik) gives the following example:

$$ A((x,y),(x_0,y_0))=\frac{i k \sqrt{(x-x_0)^2+(y-y_0)^2}}{\sqrt{(x-x_0)^2+(y-y_0)^2}} $$ is a spherical (circular) wave.

The double slit is modelt by two of these wave at $(0,-0.5)$, $(0,0.5)$. Let $k=20$, Distance to the screen $a=10$. Then $|A((a,y),(0,-0.5))+A((a,y),(0,0.5))|^2$ is claimed to be the "intensity" $I(y)$ on the screen. The distance of the maxima with these parameters is around 3, with notable values at y>10.

I have multiple issues with that:

1. "Intensity" is a directed quantity. "Intensity on the screen" implies that I could integrate over $I(y)$ to get a power. This doesn't seem to be the case, as the flux is not perpendicular to the screen. This can be easily shown by integrating over a closed circle around a single wave and varying the center. The "power" will vary notably, too. This can be fixed using a directed surface/line element in the integral, but the theory does not give directions for a double slit (for the single slit the direction is trivial).

2. This is modeled 2D, presumably because the slits are assumed to be long. Then I would also assume an intensity dependency of 1/r, not 1/r^2 in this example. Thus, I propose a power of $1/2$ to the denominator of A.

3. Even if I choose the distance to the screen much larger than the distance between the slits (I assume then I may take the flux to be perpendicular to the circumference of a circle around the center between the slits), the result does not equal the power transmitted by the two waves separately (which would be $4\pi$). With the parameters above, I get a full power of $4.6\pi$ (after the correction proposed in 2., <<1 without it), but I can change it arbitrarily when I vary the parameters (e.g $0.76\pi$ for a slit distance of .5 instead of 1). The code I used is shown below.

In the end, I would like to place another slit on the screen and calculate the pattern behind. Consequently, understanding what is happening here is crucial for me.

Here is my code in Mathematica. Note that it uses a circular screen.

r = 10;
L = .31415;
SlotRange = {-0.5, 0.5};
SlotStepW = (SlotRange[[2]] - SlotRange[[1]]);
ImageRange = {-Pi, Pi};
ImageStepW = (ImageRange[[2]] - ImageRange[[1]])/1000;
beam = Table[
   Sum[
     N[
       Exp[
        I (2 Pi/L Sqrt[(r Sin[Phi] - 
                x0)^2 + (r Cos[Phi])^2])]/(Sqrt[(r Sin[Phi] - 
              x0)^2 + (r Cos[Phi])^2])^(1/2)
      ]*SlotStepW,
    {x0, SlotRange[[1]], SlotRange[[2]], SlotStepW}
    ],
   {Phi, ImageRange[[1]], ImageRange[[2]], ImageStepW}
   ];
ListPlot[Abs[#]^2 & /@ beam, DataRange -> ImageRange, PlotRange -> All]
energy = (Plus @@ (Abs[#]^2 & /@ beam))*ImageStepW*r/Pi