# Diffraction and Computer Generated Holography Calculations

I've tried this through Mathematica, and hit my own limit in math ability trying to do this, both to no avail. I'm assuming there is no way to do so, as a simple solution to this problem would be a breakthrough for Ray Tracing, but WHY is it so difficult to integrate a euclidean distance function. (and by extension, the surrounding trigonometric functions)

$\int {\cos {\sqrt {x^2+y^2+z^2}}dx}$

I find myself REALLY wanting to accomplish this for the sake of a project I'm working on, and only being able to do this via Riemann sums makes for a VASTLY inferior calculation. The full equation is much uglier than this, but the stumbling block is this part.

Here's the full equation I'm working with at the moment. (There is a corresponding sin equation as well, the combination of which is used to determine the phase offset)

$\int_{\text{vs}}^{\text{ve}} \left(\int_{\text{us}}^{\text{ue}} \text{Cos}\left[\frac{\sqrt{(-\text{px}+\text{t4}+\text{t1} u+\text{t2} v)^2+(-\text{py}+\text{t8}+\text{t5} u+\text{t6} v)^2+(\text{t12}+\text{t9} u+\text{t10} v)^2}}{w}\right.] \, du\right) \, dv$

There are a lot of variables here because I'm a little too much computer background, and too little math background. u,v are defined over the range us-ue,vs-ve, (hence the integral ranges) and t[1-16] are components of a 4x4 translation matrix that pushes the 2d surface u,v into 3d coordinates. p[x,y] is the location of the holographic plate "pixel" being computed. w is the wavelength of light under consideration. The idea is to push the same equation to multiple cores on the GPU, each with a different pixel co-ordinate.

For the purposes of a single computation, all t,p and w variables are constants. us and ue cannot be taken for granted unfortulately as there is the possibility that us and ue are functions of v, however, vs and ve would be constants.

The idea here is to create a holographic plate representation of a surface instead of a point. I don't want to do Riemann sums because it creates holes, and I have the distinct impression that it's terribly inefficient.

Yes, this is what I originally posted to MathOverflow.

I simplified the equation above using the Fourier Aproximation, and then further reduced it's complexity using Euler's Formula. I can't speak to the efficiency of calculation of Euler's formula vs the integral of a Cos/Sin pair, but, it is technically more correct since I'm talking in electromagnetic waves anyway. The resulting Complex Amplitude actually seems to give the correct answer, as I can demonstrate through pictures more easily than words.

The speed of calculation slows down to a crawl the moment I try to do any rotations of the u,v coordinate surface, and I'm guessing that's because a whole bunch of terms zero out when it's parallel to the z=0 plane, and don't when it's not. Making the z coordinate of the surface a function of u,v also has a similar computational hit, which is unfortunate, because that's one of the easier ways to make a continuous wavefront pattern from a z-buffer.

In any event, here are the graphics I generated to test whether I was on the right track or not with my calculations, and they actually match the real wavefront pattern surprisingly well. Take a gander if you're interested:

-
The question is, as where you asked on MO, are you indeed asking for an indefinite integral in the single variable $x,$ because that is what you typed. In that case, you might as well rewrite $y^2 + z^2$ as a single symbol $k,$ since they are treated as constants in the single-variable integral. – Will Jagy Jun 10 '11 at 20:24
Full equation posted. I was trying to limit the scope of the question, but as you can see, the full equation is a double integral in u,v that's translated into x,y,z by a 4,4 matrix t[1-16] in my notation. As I say, computer background more than math background, though I do have a non-trivial amount of both. – SplinterReality Jun 11 '11 at 2:01
Why don't you have $-pz$ in the $t_{12} + t_9 u + t_{10}v$ term? – Will Jagy Jun 11 '11 at 3:23
@Will Jagy Because the holographic pixel surface (The surface the fringe pattern will be re-created on) is flat, so to make the equation simpler, albeit only slightly, the plate is aligned along the x,y axis and arbitrary placed at pz = 0 so the term drops out. Conceptually it also makes talking about the pixel coordinates easier, as it's a two term coordinate instead of three. – SplinterReality Jun 11 '11 at 4:29
So... I've spent some time working on this over the past couple of months. (It's a lot of reading for me, I'm just starting to grasp Differential EQs now) The Fresnel Integrals were on the right track, but won't get me all the way there. I've updated my math using the Fresnel Approximation, and am able to do the math quite rapidly now. Euler's Formula gave me an interesting way to generate both the real and imaginary sides of this equation as well. I'll post how that goes as far as generalizing the equation, and any interesting results I'm able to generate. Thanks for your help everyone! – SplinterReality Aug 19 '11 at 7:42

Contrary to what you typically see in a calculus course, most elementary functions don't have closed-form antiderivatives. However, I doubt that Riemann sums are the only alternative. There are much better ways of approximating integrals.

-
I'm getting a suggestion of Fresnel integrals after $k = y^2 + z^2$ and $u = \sqrt{x^2 + k},$ substitution, then integration by parts. – Will Jagy Jun 10 '11 at 18:11
I don't quite see how you can get a Fresnel integral here. If you write $u = \sqrt{x^2 + k}$, the integral becomes $\int \frac{u \cos(u)}{\sqrt{u^2-k}}\ du$. Now you could try a power series in k: $\frac{u}{\sqrt{u^2-k}} = \sum_{n=0}^\infty \frac{(2n)! k^n}{4^n (n!)^2 u^{2n}}$, and then each term can be integrated using the Si function. – Robert Israel Jun 10 '11 at 19:35
I just did parts, $t = \cos u, ds = \frac{u du}{\sqrt{u^2 - k}},$ I was left with integrating $\sqrt{u^2 - k} \sin u du,$ which wasn't great either. It might be better to say that the closest thing in Gradshteyn and Ryzhik involved Fresnel integrals but still did not do the job. Your way looks better. – Will Jagy Jun 10 '11 at 20:11
@Robert Israel: I've noticed that attempting to do math in the real world, including optics and orbital mechanics in the past. Full formula posted. I'm struggling with the math you mentioned here, but somewhat following. (I'll keep researching) Trying to translate this back to the full formula is what will likely give me the most trouble. – SplinterReality Jun 11 '11 at 1:46
Given values of the constants, numerical methods are the way to go. See the list of free implementations here: en.wikipedia.org/wiki/… – Robert Israel Jun 12 '11 at 5:55

The answer to the above is to first simplify massively through substitution.

The equation after substitution becomes simply:

$\int_{\text{vs}}^{\text{ve}} \left(\int_{\text{us}}^{\text{ue}} \text{Cos}\left[\frac{r}{w}\right] \, du\right) \, dv$

Where

r = $\sqrt{z^2+\left(-x+x'\right)^2+\left(-y+y'\right)^2}$

Run this through the Fresnel Aproximation, and:

$r \approx z + 1/2 z ((-x + x')^2 + (-y + y')^2)$

Furthermore, Euler's Formula states that:

$e^{\text{ix}}=\text{Cos}[x]+i \text{Sin}[x]$

Which is amazingly close to what I want, as I need a Cos/Sin pair for my differential. This is used to describe electromagnetic waves, of which, luckily, light is one, so I can get both numbers I want to taking the derivative of $e^{\text{ix}}$ and using the real and imaginary components. Normalizing even gets me amplitude, which is what I need to calculate anyway. Yay shortcuts.

So what I want is now: $\int\int e^{\frac{i \sqrt{z^2+\left(-x+x'\right)^2+\left(-y+y'\right)^2}}{w}} \, du \, dv$

Substitute the original values back in for x',x,y',y and z, and that's a solvable, if not nightmarish equation. Evaluate for each {px, py}, and you get the graphics I have added to this question.

Thank you for all the pointers, I couldn't have mottled through quite this much math without the suggestions from all of you.

-