Proper atmospheric laser propagation including reflections

I'm attempting to model the propagation of a laser beam (specifically, its beam radius) through the atmosphere, bouncing off a (curved) surface. This is fairly simple up to the surface itself, but I'm stuck at what to do right after the reflection. I have some ideas, but I'm uncertain about whether they're correct. Please scroll down to the bottom for hard, concrete questions!

Starting simple, for a collimated Gaussian beam the radius is given by

 $$R^2 (z)=R_0^2 +\theta(z)^2z^2$$

With $$R_0$$ the radius at the laser aperture, $$z$$ the propagation range and $$\theta$$ the beam divergence. $$\theta$$ is composed of multiple subcomponents through $$\theta^2 = \theta_{diffraction}^2+\theta_{quality}^2+\theta_{turbulence}^2+...$$. Note that some terms, such as $$\theta_{turbulence}$$ are a function of $$z$$, too.

Introducing a focus spot in this collimated beam, a small addition of the term $$(1-\frac{z}{f})^2$$ is required:

 $$R(z)^2 = R_0^2(1-\frac{z}{f})^2+\theta^2z^2$$

Now, I want to specularly reflect this beam off an obstacle that may or may not be curved. This causes two effects:

• Additional divergence term $$\theta_{refl}$$ caused by surface roughness
• (De)focusing caused by surface nonflatness (making it act as a curved mirror)

I'd implement the additional divergence term as follows:

 $$R(z)^2 = R_0^2(1-\frac{z}{f})^2+\theta^2z^2 +\theta_{refl}^2z_{refl}^2$$

Here, $$\theta_{refl}=0$$ for ranges shorter than the source-to-target range. $$z_{ref}$$ is the propagation range AFTER reflection has taken place.

However, the additional focus term caused by the curved surface seems to be harder to implement as you can't just slap another $$(1-\frac{z^*}{f^*})^2$$ on the equation. I have considered:

Switching to the ABCD-matrix approach

When using Gaussian beams in these matrix calculations, you need the complex beam parameter q:

 $$\frac{1}{q} = \frac{1}{C(z)}-i\frac{\lambda}{\pi w(z)^2}$$

with $$C(z)$$ the radius of curvature of the phase front

 $$C(z) = z[1+(\frac{z_R}{z})^2]$$

and

 $$z_R=\frac{\pi w_0^2}{\lambda}$$

However, I don't think this approach is fully valid as it assumes an ideal scenario, e.g. beam quality or turbulence effects on the divergence are not taken into account. Could this be tackled by using the definition of beam divergence $$\theta=\frac{\lambda}{\pi w_0}$$ to calculate a 'representative' value for $$w_0$$ (at any given location) to be used to calculate $$z_R$$ and then $$C(z)$$ and then $$\frac{1}{q}$$?

Still, this would mean I have to propagate in small increments to any propagation range I like in order to stay accurate - I can't just use a single equation such as equation .

Splitting equation 

Rather than trying to cram everything in equation , can I split it in both terms? The first part ($$R_0^2(1-\frac{z}{f})^2$$ corresponds to purely geometrical optics, whereas the second part $$\theta^2 z^2$$ corresponds to 'defects' that hinder an infinitely small spot size and add some divergence. Can I separately calculate both terms, e.g. calculate spotsize according to geometric optics for the first term (by using ABCD-matrices) and then later 'adding' the (atmospheric) divergence components later using $$\theta^2 z^2+\theta_{refl}^2 z_{refl}^2$$? Is that allowed?

All in all, I feel like there should be a simple(r) solution, but I'm afraid I've been staring at this for too long. Any help is greatly appreciated!

Concrete questions

1. Is equation  indeed valid when only considering the additional divergence components?
2. Can the ABCD-matrix approach be used when (ab)using the current divergence to calculate a representative w_0? (If yes, I have a solution, but it wouldn't be my favorite)
3. Would it be allowed to split equation  so that the first part can be tackled using the ABCD approach and the second part can be added later?
4. Are there any simpler ways to doing this?