Plot recipes; Docs; RTM of a system; waistradiusfunc

[ghost]

This PR achieves many things:

1. Documentation is generated via Documenter.jl via a separate project in directory \docs
2. The plot meta package Plots.jl is supported: Plotting works and is documented with 2 examples
3. A method to calculate the Ray Transfer Matrix of an entire system is added to the function RTM
4. The function waistradiusfunc is added to provide access to a function returning a system and beam's waist radius as a function of the location along the beam axis

Opionion on the sensibility of including this may vary. There are disadvantages:

1. Several external dependencies are required, making the installation of this package less "light"
2. Documentation generation even adds an entire new project (to avoid making package Plots a requirement for the basic package itself)
3. Perhaps the interface to waistradiusfunc is unwanted?

[ngedwin98]

Thank you for the pull request. I like the first step towards better documentation, as well as the idea of support for plotting of beam profiles. I would like to start a conversation regarding waistradiusfunc:

My current (very basic) approach to plotting has been to generate a desired optical system, call discretize on it to create a more finely meshed system, evaluate it, and then plot it. For example, with some system::Vector{<:Element}, an initial beam0::Beam, and using PyPlot

system = discretize(system, 100)
beam_trace = beamtrace(system, beam0)
plot(location.(beam_trace), spotsize.(beam_trace))

I have found this to be more or less sufficient for visualizing a particular optical design.

However, I suppose it might be desirable to evaluate an optical system at arbitrary locations which are not known to the user a priori (as could be imagined for an optimization routine). Is this the use case you have in mind?

In this case, I think the implementation could be improved by taking in a non-discretized system, and then using logic to find the element of FreeSpace which contains the desired location. We could then split that element exactly at the location desired into two elements freespaceleft and freespaceright, and then call beamtrace up to freespaceleft. This has the advantage of not calling beamtrace on elements that are not relevant, as well as avoiding re-implementing parts of beamtrace which appear in the current function. It would also be a good opportunity to redesign the "discretize" functionality into something more sensible, such a "freespace partitioner" which might live only within the context of this "evaluation" model -- after all, it seems like its only purpose is for visualization and evaluation of a beamtrace.

Next, I am somewhat wary of writing function-generating code as in this case. In addition to issues with closures, I find that they tend to be more application-specific in nature. Perhaps we can devise an interface function which does not produce a specific function with an application-specific function signature, but which is still close enough to the desired use case that the calling code can do something like:

f = let system=system, beam0=beam0
    z -> evaluate_spot_radius(system, beam0, z)
end

(The name, of course, is only for demonstration purposes.) If closure performance is not an issue, it's a simple one-line replacement to waistradiusfunc, as in z -> evaluate_spot_radius(system, beam0, z).

Additionally, I think in this situation, there is a need to implement both a method for evaluating a beamtrace at a single location and for evaluating at an array of locations. This is mostly for performance, since (as you've already anticipated in your implementation) both the elements and the list of desired locations (in principle) are sorted. However, I am not sure which is in the majority: situations in which evaluation occurs at a single location (though perhaps called iteratively), or situations in which evaluation occurs at multiple locations at once.

Finally, in my experience, the term "waist" refers specifically to focus of the beam (where the q parameter is strictly imaginary). Elsewhere, the 1/e^2 radius is usually just called the spot size (or rather, the spot radius, as you've correctly pointed out).

I'm curious to hear your thoughts on this matter, as I hope we can support this functionality with the above points addressed. (Unfortunately, it seems like the documentation and plotting functionalities are tied up with this at the moment so I'm not sure about pulling them in separately.) Let me know if you have the interest and time to look into these points. Otherwise, I'll try to do some implementation within the week, subject to further points on this discussion.

[ghost]

Thank you for the helpful critique, and for inviting discussion.

BTW, I ended up not using waistradiusfunc for the plotting code after all, although I had intended to: The original idea was that through the plotting libraries function evaluation code, it would better represent the beam in the neighborhood of beam waists (using your correct definition). I also wonder if (and perhaps incorrectly considered that) such a function would be useful when trying to do some symbolic evaluation, e.g. for cavity design, which I haven't tried yet but would make for a neat example for the documentation. I suppose, however, that everything I wanted to do with symbolic evaluation could even benefit from avoiding the if-then logic associated with finding the right boundary elements and instead using those directly.

I do not favor the idea of simply calling discretize (at least in its current form) to prepare plot data because I anticipate that some optical systems will have both very short and very long [Free]Space segments that do not necessarily contain a waist and hence may not need to be subdivided using the same number of subdivisions. In particular, I'm toying with the idea of adding support for optical plates (or crystals) and it seems stupendous to say, subdivide a 1mm thick plate the same way as a 1m free space segment. One approach that I have considered but not found perfect, either, would be to give a keyword option to discretize that either gives a target space between elements or somehow relates the spacing to the change in wavefront curvature or distance to a beam waist (to be more sensitive to the neighborhood of beam waists).

Thank you for pointing out the performance issue with closures which I had once read but forgotten about. Note that this issue can be addressed with local copies of the variables, as you already incorporated in your example.

In general, I may have to rethink my approach to using this code, as I anticipate needing a function such as waistradiusfunc and would prefer the brevity of a specifically generated function f(z) to calling a calculateXXX(system_or_element, beam, z). However, my shorthand does not need to live in this package, although I would prefer to have it there (and likely even a function that can also handle 2d or 3d coordinates to return an off-axis, relative intensity).

Regarding the name you make a great point; I think I have somehow mangled the concept of 1/e^2 intensity boundary with the name waist which obviously, as you point out, should be reserved to the actual beam waist. Of course you convinced me of renaming waistradius (and waistradiusfunc, if that survives at all); what do you think of spotradius?

[ngedwin98]

Thanks for the detailed responses.

I agree that discretize in its present form is not the best way forward, for many of the reasons you gave. (Incidentally, in the initial planning stages, I was deciding between a fixed discretization length vs. a fixed number of elements to split into, and I decided on the latter for not-particularly-compelling reasons.) I think that once a better way of determining where and how to evaluate the beamtrace, we may even be able to potentially retire discretize.

For the time being, let me merge this PR as-is. We can decide what to do with waistradiusfunc separately, but since it's a functional piece of code, I don't mind introducing it into the package for the time being (even if we might remove it at a later time).

After merging this PR, I should be current with your fork up to 393fea3. Since I anticipate continuing to pull your changes, it may be convenient for us to sync up. Would you like to have a discussion in finding a good workflow for us to do so?

[ghost]

Since I anticipate continuing to pull your changes, it may be convenient for us to sync up. Would you like to have a discussion in finding a good workflow for us to do so?

Yes, I very much welcome that. What do you suggest?

[ghost]

BTW, I'm halfway through type parametrization and stumble on the field x in Beam. What is its intended purpose? I'm afraid I can't quite guess it.

[ngedwin98]

BTW, I'm halfway through type parametrization and stumble on the field x in Beam. What is its intended purpose? I'm afraid I can't quite guess it.

The field x is supposed to represent the transverse (with respect to the optical axis) component of the paraxial ray. It's the first element of the ray vector in Siegman's book (the other being the slope).

How do you feel about the idea of using Gitter to continue chatting about points of clarification and other workflow topics?

[ghost]

I've never used Gitter, but sure, I can sign up for our discussion. I'll report back when that's done.