Method of Images - Point Charge with Semi-infinite Dielectric To calculate the electric field for the system where a point charge is placed at a distance d from a semi-infinite dielectric medium (eg. http://farside.ph.utexas.edu/teaching/jk1/lectures/node42.html), we use the method of images: one image charge q1 in the dielectric medium with distance a to the interface for calculating the field in vacuum; one image charge q2 in vacuum with distance b to the interface for calculating the field in the medium.
However, how do we determine the locations (a and b) of the image charges?
The sources I've seen (including the one above) just simply state that the a=b=d without proof. (This one gives a proof but I think its calculation details are incorrect: https://physics.princeton.edu/~mcdonald/examples/image.pdf)  
I've tried to use the boundary conditions to determine the locations but we have 4 unknowns (q1, q2, a, b) and only 3 equations from the boundary conditions:
(1) the potential is continuous at the interface
(2) the tangential components of the electric fields are the same at the interface
(3) the normal components of the displacement fields (D=epsilon*E) are the same at the interface.
Therefore, I only found the relationship between a, b, and d instead of exact values, and can express q1 and q2 in terms of a, b, q and d.
So I guess there must be some reasons that I'm not aware of for choosing a=b=d?
 A: 
In the notion of Jackson 4.4, assume the original charge located at $z = d > 0$. For the right-half plane, the image charge of the left-hand side ($\epsilon_2$) is $q_2$ while for the left-half plane, the image charge of the right-hand side ($\epsilon_1$) is $q_1$.
Let's first pretend that we don't know the exact position of $q_1$ and $q_2$ and therefore denote them as $z=d_1$ and $z=d_2$.
Matching all the boundary condition ($D_z$, $E_\rho$ continuous) yields
$$
\frac{qd}{(\sqrt{\rho^2 + d^2})^3} - \frac{q_2d_2}{(\sqrt{\rho^2 + d_2^2})^3} =  - \frac{q_1d_1}{\sqrt{(\rho^2 + d_1^2})^3}
\\
\frac{1}{\epsilon_1} \left( - \frac{q}{(\sqrt{\rho^2 + d^2})^3} - \frac{q_2}{(\sqrt{\rho^2 + d_2^2})^3} \right) = \frac{1}{\epsilon_2} \left( - \frac{q_1}{\sqrt{(\rho^2 + d_1^2})^3} \right)
$$
Looks like we are running out of criteria to solve these two equations. However, if we want the method of image charge to work, there's a hidden requirement: The set of $q_1, q_2, d_1, d_2$ should work across the entire region of interest, i.e. the final equation we are looking for should not contain $\rho$ and $z$ dependence, otherwise, it means the method of image charge fails because the image charge setup we presumed is simply wrong (or method of image charge is simply not applicable to this problem).
Back to this problem, we need to find a combination of $q_1, q_2, d_1, d_2$ that gets rid of $\rho$ and $z$ dependence. This is the third criterion besides the two boundary conditions. And apparently, setting $d_1=d_2=d$ works so all denominators cancel out.
A: I finally figured out what constraints are missing apart from the boundary conditions: All values of the imaginary charges should be constants, i.e. no dependence on the coordinates e.g. the radial component "r".
For example, solving the boundary conditions only gives the values of the imaginary charges as functions of {the source charge, source charge location, locations of imaginary charges, r }. If one imposes the constraint mentioned before, the solution that all distances are equal naturally drops out.
