Thanks for pointing out the article that discusses use of phase space plots to help with the design of an illumination system.  It’s interesting work.

In looking over this paper, it appears that the authors use phase space plots to help better understand the fundamental behavior of the starting optical path (in this case a reflective ellipsoidal collector).  By looking at the initial phase space plot in the illumination plane, they realize what is needed to improve irradiance uniformity is a two-surface freeform corrector optic that helps suppress variation of the magnification over the pupil.  It doesn’t appear that they actually use the entire phase space data directly for optimization purposes (although perhaps a subset may be utilized), but after optimization they show that the new phase space plot is consistent with the improved irradiance uniformity they achieved.

In any event, you asked how one can go about generating phase space plots programmatically.  Since, in general, phase space is 4D, it is not possible to plot the complete version.  However, for axially symmetric systems we can just use a 2D version (like the authors of the paper did).  To get this data, we can use ZOS-API to trace a set of rays emanating from a line in the object (source) plane.  For each starting x location we can define a fan of rays in angle (a sequence of L values, i.e., a set of x-direction cosine values); the complete set of source rays can then be used for batch ray tracing. 

Here is a sequential model that replicates the starting system in the paper:

Next, we need to write a script for batch ray tracing (if needed, I can provide more detail, but see example #22 in the ZOS-API help documentation).  I did this in Matlab, and here is what I find by simply plotting L vs. x for the rays in the object plane and the image plane.

This yields good agreement with the plots shown in the paper:

Hope this helps…

Regards,

Jeff