Huygens wave theory not applicable to lasers or parallel beams of light? According to Huygens wave theory, every point on a wavefront acts as a secondary source of waves.

Using this principle we can never have pretty narrow parallel beams of light right? Like lasers? There will always be "leakage"(wrong terminology?) of waves from the edge of the beam right?

 A: This is wrong. A wide region (width $\gg\lambda$) of many, evenly distributed in-phase spherical wave emitters, like the Huygens equivalent of a laser beam's cross section, behaves like a phased antenna array, with the result that destructive interference between the emitters cancels radiation deviating significantly from the beam and reinforces radiation along the beam. So you do indeed get arbitrarily narrow beams from Huygens's theory.
What you might be thinking is that this theory does not cope well with a unidirectionally propagating beam. If we replace the laser beam's cross section in-phase with spherical emitters, a narrow beam is radiated both forwards and backwards (a phased antenna array does this, for example). This is where one must call on an obliquity factor to restore unidirectionality. Fresnel-Huygens theory multiplies the spherical radiation pattern by $\frac{1}{2}(1+\cos\theta)$, where $\theta$ is the angle between the direction of propagation and the line joining the centre of the Huygens wavefront and the point in question. 
If you look at the Huygens-Fresnel diffraction integral in the Wikipedia link above, then consider the beam cross-section $B$ in the $x-y$ plane and then also a point $P$ with co-ordinates $(x, y, z)$ in the farfield (i.e. so that $R = \sqrt{z^2+y^2+z^2} \gg w$, where $w$ is the maximum width of the beam cross section, then a Huygens superposition of spherical emitters spread over $B$ will beget a disturbance at $P$ of:
$$\int_B \frac{\exp\left(i\,k\,\sqrt{\left(x-x^\prime\right)^2+\left(y-y^\prime\right)^2+z^2}\right)}{\sqrt{\left(x-x^\prime\right)^2+\left(y-y^\prime\right)^2+z^2}}\,h\left(x^\prime,y^\prime\right)\,\mathrm{d}x^\prime\mathrm{d}y^\prime\approx $$
$$\frac{e^{i\,k\,R}}{R} \int_B \exp\left(-i \frac{k}{R}\left(x\,x^\prime+y\,y^\prime\right)\right)\,h\left(x^\prime,y^\prime\right)\,\mathrm{d}x^\prime\mathrm{d}y^\prime$$
where $h\left(x^\prime,y^\prime\right)$ represents any beam apodisation and aberration (intensity variation and phase, respectively) and I have simply pulled the denominator out of the integral: this is justified because, as a proportion of itself, $\sqrt{\left(x-x^\prime\right)^2+\left(y-y^\prime\right)^2+z^2}$ does not vary much for $\left(x^\prime,y^\prime\right)\in B$ as $R\rightarrow\infty$, but, as a number of wavelengths it varies a great deal, so we must keep the numerator. Thus we see two results:


*

*The expression is symmetric in $z$, i.e. the same disturbance propagates backwards as it does forwards, so that the beam is not unidirectionally propagating. This is why Fresnel introduced his obliquity factor $K(\chi)$ to remove the backwards propagating wave;

*The farfield beam cross section is given by the Fourier transform of the beam cross section at $B$, but with a scale factor $\frac{k}{R}$ in the Fourier transform's argument, i.e. once we get a reasonable distance from the cross section at $z=0$, the "shape" of the projected beam does not change - it is always equal to a Fourier transform of the beam at $z=0$, but the same shape is dilated by a factor proportional to $R$.


This diffraction scheme is called Fraunhoffer diffraction. So the diffracted beam can be made tightly focussed in the farfield if the cross section at $z=0$ is wide. By dint of the Fourier transform, there is an inverse proportionality between the width of the beam at $z=0$ and that in the farfield. This is exactly the behavior of a high quality laser beam - you will always get a dilation proportional to $R$ - the only way to avoid this is to have a plane wave of infinite extent.
To get some intuition for the result, witness that the above theory gives, for the diffraction of the Gaussian beam defined by:
$$h\left(x^\prime,y^\prime\right) = \exp\left(-\frac{{x^\prime}^2+{y^\prime}^2}{2\sigma^2}\right)$$
the diffracted result at the point $P=(x,y,z)$:
$$\frac{e^{i\,k\,R}}{R} \int_B \exp\left(-i \frac{k}{R}\left(x\,x^\prime+y\,y^\prime\right)\right)\,h\left(x^\prime,y^\prime\right)\,\mathrm{d}x^\prime\mathrm{d}y^\prime  = \frac{\sigma^2}{R} \exp\left(-\frac{k^2 \sigma ^2 \left(x^2+y^2\right)}{2 R^2}\right)$$
