LightPipes: beam propagation toolbox

Gleb Vdovin
OKO Technologies
Rontgenweg 1, 2624 BD Delft, The Netherlands
Phone: 31-15-2851547, email:
©1993-1999, Gleb Vdovin


Command line interface


LightPipes is a portable set of Unix-style tools written in C. It is designed to model the propagation of light in coherent optical devices in scalar approximation. The toolbox is based on efficient propagation algorithms, providing extended possibilities of beam manipulation. Simulations of interferometers, holographic setups, laser resonators, lasers, Fourier optics, waveguides are possible. The source code is portable, it can be compiled on any low-end workstation or a decent PC under Unix, MSDOS or OS2. The MSDOS binary is reported to run also in DOS box under OS2 and W95.

The toolbox consists of a number of filters. Each filter represents an optical element or a step in the light propagation. The components of a model can be connected via Unix system pipes in any order in (any) shell environment to form the model of an optical system. The syntax depends on the shell used.

All filters share a huge data structure, containing sampled two-dimensional distribution of the scalar complex amplitude of light. These structures can be saved and retrieved to/from the disk at any step of computation. The number of light beams treated simultaneously is limited only by the disk space available. Tools to generate presentation quality Postscript density plots and to export portable gray-map files (PGM graphic format with magic number P2 or P5), surface wireframes (Xplot or Gnuplot required) and contour plots (Gnuplot required) are included.

C shell syntax is used for all examples in this manual. A good reference for the syntax of C shell is given in  [1]. Any other shell environment or script language can be used with LightPipes as far as binary system pipes are supported.

The package is far from being perfect, there are more powerful and flexible (also much more expensive) programs around. LightPipes is written with the idea of maximum functionality with minimum of coding effort. The development was based on the tradeoff concept ``first 10% of coding effort give 90% of the functionality''.


There is no warranty for the program, to the extent permitted by applicable law. Except when otherwise stated in writing, the copyright holders and/or other parties provide the program "as is" without warranty of any kind, either expressed or implied, including, but not limited to, the implied warranties of merchantability and fitness for a particular purpose. The entire risk as to the quality and performance of the program is with you. Should the program prove defective, you assume the cost of all necessary servicing, repair or correction.


LightPipes beam propagation toolbox is copyright © (1993-1999) Gleb Vdovin. This manual (also © Gleb Vdovin) may be freely distributed and copied as a whole.

Unix/DOS source of LightPipes is distributed under GNU public license. Contact: Dr. Gleb Vdovin, OKO technologies, Reinier de Graafweg 300, 2625 DJ Delft, The Netherlands,,, +31-62-2794647 voice, +31-15-2574233 fax.

New features

LightPipes.1.2 introduces finite-difference propagation routine steps and a better (cubic two-dimensional) interpolation procedure interpol. The old linear interpolation procedure remains with a new name interp1. Steps can be used to model waveguides with complex three-dimensional distribution of refraction/absorption. The new distribution also includes more tools operating on intensity/phase distributions. These service filters are not considered to be the part of LightPipes, though they extend the package possibilities. Tools are distributed ``as is'' free with source. These sources can be used as templates for customization.

More bugs were fixed, for example now lens_forvard can correctly propagate behind the focus.


Install the executables into a directory included into your search path, or link them to /bin. For MSDOS distribution use pkunzip -d to create the directory tree.

There is a set of tools, to deal with intensity and phase distributions exported by LightPipes. These tools are freeware distributed with sources. There is no manual for tools. To build tools, type make tools at the command line. You may write your own tools by modifying existing sources.

To build the full package from the source, put all source files into an empty directory. Be sure the name of the directory is included into your search path. Edit the makefile: define the C compiler installed in your system. Gnu C compiler (gcc -O2 -Wall) is the development compiler, chosen for its portability, high quality and availability on the net. Other C compilers will also work but there is no guarantee. Very probably the package compiled with the native compiler will work a bit faster. Type make. To build MSDOS binaries with DJGPP (version 1.12) type make -f makedos. DJGPP is available from numerous Simtel mirrors, directory vendors. In the case of trouble first consult your system administrator, then send email to (valid in 1996 and perhaps later). To use LightPipes from other directory, include the directory containing LightPipes executables into your search path, or link all executables to a directory which is included into your search path, for example to /home/joe/bin.

First steps

Prompt help

Any LightPipes filter invoked with no command line arguments prints short syntax description to stderr.

The package is under development, the behavior and names of commands can be different from described here. Some features can be missed in this description. In all these cases please refer to the prompt help printed by the command itself. This message is supposed to reflect the very last modifications.

Incidentally some commands of your system may have the same names as the LightPipes filters (like, for example, begin). LightPipes filter causing the problem may be renamed. There are also other ways:

  1. Redefine the PATH environment to make shell search the LightPipes directory first.
  2. Use the full name of a command: like ./begin or /home/joe/LightPipes/begin.

Simple models

The simplest way to invoke LightPipes is to type the commands at the system (Unix/MSDOS/OS2) prompt. All examples from this manual may be executed in interactive mode, i.e. you type the command at the system prompt and the computer executes it. It's a serious advantage, allowing for building simple models in usual shell environment. Sometimes the chain of piped commands is too long to be printed in one line, or doesn't fit into the line buffer. Long lines can be broken, like here:

program_a | program_b | program_c > /dev/null
which is equivalent to:
program_a > foo
program_b < foo  > foo1
program_c < foo1 > /dev/null
rm foo foo1
In the last example the temporary files foo and foo1 were created to store the field distribution between executions of commands. After executions, the temporary files should be deleted (with Unix command rm). Fifos (named pipes) may be used instead of files. Before being used, the fifo must be created with a command mkfifo. Fifos are faster, consume more memory, they are not supported by some systems and, if you write to the fifo, the next command must read it immediately. There are too many rules about fifos, writing to a file is much more flexible and reliable operation.

Other way to use LightPipes is to put commands in a shell script. Simple (C-shell) script has the following structure:

#!/bin/csh -f 
LightPipes program ....
File containing a script must be readable and executable under Unix.

To begin with, the data structure for LightPipes must be created: we type command begin at the shell prompt. Invoked with no arguments begin prints to the screen the help message:

begin forms the initial data structure for use by all following programs

USAGE: begin B C [A],  where:

A is the grid dimension, grid of AxA points will be formed (if omitted, A=256)
B is the grid size, the real size of the grid will be put as BxB [units you use]
C is the wavelength in [units you use]
Three parameters: must be defined.

Obligatory for begin are only the size and the wavelength, the grid is 256x256 points by default (or smaller in the demo version). For example:

begin 0.01 1e-6 128 > /dev/null
defines a plane wave with a unity amplitude and wavelength of tex2html_wrap_inline1703 in 128x128 square grid with a size of 1x1cm.

The standard input and the standard output channels are used to pass huge binary data structures, containing the field complex amplitude, between filters. The redirection to dev/null is used in a Unix environment to prevent the binary output to the screen. Redirection to files also works:

begin 0.01 1e-6  > first_pipes_file
In the last example the field structure generated by the program begin is written to the file
first_pipes_file. Begin is the first block in any LightPipes program, it generates the basic structure and accepts no piped input. All other LightPipes programs accept piped input and can be connected by system pipes in any order.

The field structures in LightPipes may be defined only in square arrays. It is not always optimal, sometimes much more memory is used than necessary. The default structure produced by begin (256x256 points) is 1MB large, and it is so, even if only in a one-dimensional array of of 256x1 points is of the interest. This approach simplifies coding and will not be changed in future.

Begin allocates memory for a square array of complex numbers, fills it with values corresponding to a uniform distribution of intensity and phase (plane wave) and redirects the header (containing information about sampling/wavelength) and the formed structure to the standard output in binary, computer-dependent format. Binary files written by LightPipes are not portable in general, because two different types of data presentation (big/little endian) exist. Binary files can be moved between computers with similar endians, thus PC/Linux -> HP is forbidden. If the output is not redirected to a file or to /dev/null, it fills the screen with garbage and may cause other problems.

LightPipes can be started with different grid samplings, for example

begin 0.01 1e-6  500 >second_pipes_file
will create 4Mb file containing 500x500 structure.

These examples illustrate the amount of RAM needed to run LightPipes. More than one megabyte is needed for 256x256 grid, 4MB for 512x512 , 16MB for 1024x1024 and a really large system is necessary for 2048x2048. These numbers may be slightly different for different computers. Some filters allocate additional memory (up to 8 times of the basic data structure in the worst case), rising further the memory requirements.

If RAM is not sufficient, the system will swap to the hard-disk, slowing down not only your job, but the whole computer. The 512x512 size is safe for the majority of Unix systems.

The data structure of the demo version is restricted to less than 1MB, thus the demo is very fast and runs on any computer.

Dimensions of data structures

LightPipes works with any even sampling, i.e. 300, 196, 578...

>begin 1 1 99 | rect_ap 1 > /dev/null
>Sorry, number of points must be even, stopping
illustrates this feature. Program does not work with odd samplings.

Apertures and screens

The simplest component to model is an aperture. There are three types:

  1. circular aperture: circ_ap
  2. rectangular aperture rect_ap
  3. Gaussian diaphragm gauss.

circ_ap filters the field through a circular aperture of radius R given as the first command line argument. Two optional arguments define the decentering of the aperture in X and Y directions.

rect_ap filters the field through a square or rectangular aperture. The first command line argument defines the side length of the (square) aperture, the second is the length of the second side of a rectangle, third and fourth are shifts in x and y directions and the last argument is the rotation of the aperture. Only the first argument is obligatory (centered square aperture).

gauss filters the field through a Gaussian intensity filter. The transmitted intensity is given by:


where R is the radius of the Gaussian aperture and A is the maximum transmission. This filter can be used to model gauss beams: begin followed by Gauss produces a finite grid approximation of a gauss beam.

There are three commands describing screens: circ_screen for a circular screen (inversion of the circular aperture), rect_screen is the inversion of the rectangular aperture and gauss_screen is the inversion of the gauss. The syntax for screens is exactly the same as for apertures. The gauss_screen is complementary to gauss, so that the sum of transmitted ( gauss and ``reflected'' gauss_screen intensity equals to the incident intensity. These two filters are useful for modeling of gauss apodized output couplers. gauss should be used to model the beam reflected back into the resonator by an apodized mirror while gauss_screen will produce the output beam which is filtered through the apodized mirror.

Examples for apertures and screens

Here are examples of simple LightPipes scripts demonstrating the usage of apertures and screens:

  1. begin 0.01 1e-6  | circ_ap 0.0025 | file_ps  > /dev/null
  2. begin 0.01 1e-6  > foo
    circ_ap 0.0025 0.001 0.002 < foo | file_ps > /dev/null
  3. begin 0.01 1e-6   | circ_ap 0.0035 > foo
    rect_screen 0.002<foo | file_ps > /dev/null
  4. begin 0.01 1e-6  | gauss 0.004 >foo
    rect_screen 0.001 0.001 -0.0015 -0.002 < foo > foo1
    rect_screen 0.001 0.0035 -0.002 0.0025 45 < foo1 > foo2
    circ_screen 0.0005  0.0025 -0.003 < foo2 > foo3
    circ_screen 0.0007   0.001 0.0015 < foo3 | file_ps 128 > /dev/null
    rm foo foo1 foo2 foo3

The resulting intensity distributions are shown in Fig. 1.

Figure 1: Intensity distributions produced by four examples (from left to right) and saved to postscript files out[1...4].ps

When a script runs at the background - and it is frequently the case, the following trick shows to the screen the intensity distribution of your beam:

#!/bin/csh -f
cat $1 | file_ps temp_ps > /dev/null
ghostview temp_ps
# uncomment this if ghostview executable is called gv:
# gv temp_ps
rm temp_ps

It's named ``show'' and that's exactly what it does. You can use it in your scripts to view the intensity distributions. You need ghostview installed to run this script.

Graphic output

In the previous example a new operator file_ps was introduced. file_ps is one of two LightPipes programs providing the possibility to output graphics. It writes the intensity distribution (density plot) into a postscript file. The program accepts three command-line arguments: the filename, the grid dimension and the value of contrast correction constant tex2html_wrap_inline1721 . The dimension and tex2html_wrap_inline1721 are optional, defaults are 128x128 (different in the demo) and tex2html_wrap_inline1725 (it increases slightly the brightness in low-intensity regions, showing the distribution in a more natural way). The postscript output of the previous example (Fig. 1) is 128x128 and occupies approximately 32K on disk (different in demo). When tex2html_wrap_inline1727 sampling is used in calculations, the sampling of postscript graph must be also tex2html_wrap_inline1727 i.e. 32,64,128,256,512.... The filter does not report a wrong sampling, instead it adjusts automathically the sampling to nearest sensible value.

The second program which outputs graphics is file_pgm.

file_pgm writes intensity into *.pgm file F

USAGE: file_pgm F [N, G, MAX], where F is the output filename,
optional parameter N is the grid size, N=128 if omitted in command line
G is the $\gamma$ parameter, [0.1...10], higher G gives better
contrast in low intensities, default G=2.0
MAX is the number of gray levels, default MAX=255

Output file F can be processed with netpbm package

Files (intensity distributions) produced with file_pgm can be imported by a popular X graphical program XV, which saves to many other formats and allows for different conversions. Pgm stays for ``portable graymap'' and can be converted to other formats with netpbm or pbmplus packages. Files (intensity distributions) produced with file_pgm can also be imported by LightPipes using fil_ter. The tex2html_wrap_inline1721 parameter for files that are saved for further use by LightPipes has to be set to 1 (as mentioned before, the default is 2).

There are possibilities of graphical interface to gnuplot:

file_int writes the intensity distribution into a file which may be plotted with gnuplot graphical package. One command line argument is necessary: the output filename. Default grid size is 64x64, you may change the default dimension with the second command line argument. file_pha is similar to file_int, it writes wrapped phase distribution. Surface plots are not always acceptable - sometimes two-dimensional plots of the cross section of the light beam are of interest. It is especially actual for wrapped phase plots, which look awful in three dimensions. The command cros_out writes into a file the X and Y cross sections of the field distributions. For example the program:

begin 0.04 1e-6 | gauss 0.01 0.003 -0.002 > foo
cat foo | file_int in | file_pha pha 128 | cros_out out >/dev/null
rm foo
will form decentered Gaussian intensity distribution in a 256x256 grid and save the intensity distribution in a 64x64 file in, the phase distribution in file pha and cross sections in file cros_out. File in is 64x64 (by default), file pha is saved in a grid having 128x128 points. The number of points in the output file is governed by rules, described already for the postscript output.

Gnuplot has to be installed on your system to plot the distributions: If it's installed, type


and you will get on your screen something like:

        G N U P L O T
        Linux version 3.5
        patchlevel, 27 Aug 93
        last modified Fri Aug 27 05:21:33 GMT 1993 
        Copyright(C) 1986 - 1993   Colin Kelley, Thomas Williams 
        Send comments and requests for help to
        Send bugs, suggestions and mods to
Terminal type set to 'x11'

At the gnuplot prompt type:

gnuplot> set hidden
gnuplot> splot 'in' with lines
gnuplot> splot 'pha' with lines
gnuplot> plot  'out' using 1:2 with li,  'out' using 1:3 with li
gnuplot> plot  'out' using 1:4 with li,  'out' using 1:5 with li

The first operator plots the intensity plot, second- the phase plot, the third and the fourth plot the Y intensity cross-sections and the Y phase cross-section. The graphic output is shown in Fig. 2.

Figure 2: Graphs produced with gnuplot: Surface plot intensity, surface plot phase and intensity cross sections

Files saved with file_pha and file_int can also be imported with fil_ter, which automatically recognizes pgm and gnuplot formats.


There are six (four different algorithms, two of these four allow for coordinate conversions) different possibilities for modeling the light propagation in LightPipes.

FFT propagation (spectral method).

Let us consider the wave function U in two planes: U(x,y,0) and U(x,y,z). Suppose then that U(x,y,z) is the result of propagation of U(x,y,0) to the distance z, with the Fourier transforms of these two (initial and propagated ) wave functions given by tex2html_wrap_inline1745 and tex2html_wrap_inline1747 correspondently. In the Fresnel approximation, the Fourier transform of the diffracted wave function is related to the Fourier transform of the initial function via the frequency transfer characteristic of the free space tex2html_wrap_inline1749 , given by [4, 5]:





Expressions (2, 3, 4) provide a symmetrical relation between the initial and diffracted wave functions in the Fresnel approximation. Applied in the order (3tex2html_wrap_inline1751  (2) tex2html_wrap_inline1753  (4) they result in the diffracted wave function, while being applied in the reversed order they allow for reconstruction of the initial wave function from the result of diffraction. We shall denote the forward and the reversed propagation operations defined by expressions (2, 3 and 4) with operators tex2html_wrap_inline1755 and tex2html_wrap_inline1757 respectively.

The described algorithm can be implemented numerically using Fast Fourier Transform (FFT) [5, 6] on a finite rectangular grid with periodic border conditions. It results in a model of beam propagation inside a rectangular waveguide with reflective walls. To approximate a free-space propagation, wide empty guard bands have to be formed around the wave function defined on a grid. To eliminate the influence of the finite rectangular data window, Gaussian amplitude windowing in the frequency domain should be applied - see [5, 6] for extensive analysis of these computational aspects.

The simplest and fastest LightPipes command for propagation is forvard. It implements the spectral method described by (234). The syntax is simple, for example if you want to filter your field through 1cm circular aperture and then propagate the beam 1m forward, you type:

begin 0.02 1e-6 |circ_ap 0.005 |forvard 1 > foo 
file_int in < foo |cros_out out |file_ps 128 >/dev/null 
rm foo

Figure 3: The result of the propagation: density, surface and cross section intensity plots

We see the diffraction effects, intensity distribution is not uniform anymore.

The algorithm is very fast in comparison with direct calculation of diffraction integrals. Features to be taken into account:

Due to these two reasons the operator
begin 0.02 1e-6 | rect_ap 0.02 | forvard 1 > out
makes no sense. The cross section of the beam (argument of rect_ap) equals to the section of the grid (the first argument of begin), so we have a model of light propagation in a waveguide but not in a free space. One have to put
 begin 0.04 1e-6 | rect_ap 0.02 | forvard 1 > out
for propagation in the near field, and may be
begin 0.2 1e-6 512  | rect_ap 0.02 | forvard 400 > out
for far field propagation. Comparing the result of previous example with the result of
begin 0.06  1e-6 512  | rect_ap 0.02 | forvard 400 > out
we see the difference.

Described FFT algorithm has its drawbacks. The positive message is that it is very fast, works pretty well if properly used, is simple in implementation and does not require any allocation of extra memory. In LightPipes.1.1 and later a negative argument may be supplied to forvard. It means that the program will perform ``propagation back'' or in other words it will reconstruct the initial field from the one diffracted. For example :

#!/bin/csh -f
# initial uniform distribution:
begin .1 1e-6 |circ_ap 0.025|file_ps> foo
#diffraction 30 m
forvard 30 <foo |file_ps >foo1
#reconstruction of the initial field: 
forvard -30 < foo1 | file_ps >/dev/null
rm foo foo1
The result produced by this script is shown in Fig. 4.

Figure 4: The initial field, the field propagated to the near field and the field propagated back (from left to right)

Direct integration as a convolution: FFT approach

Another possibility of a fast computer implementation of the operator tex2html_wrap_inline1755 is free from many of the drawbacks of the described spectral algorithm. The operator tex2html_wrap_inline1755 may be numerically implemented with direct summation of the Fresnel-Kirchoff diffraction integral:


with functions U(x,y,0) and U(x,y,z) defined on rectangular grids. This integral may be converted into a convolution form which can be efficiently computed using FFT [7, 8]. This method is free from many drawbacks of the spectral method given by the sequence (3tex2html_wrap_inline1751  (2) tex2html_wrap_inline1753  (4), although it is still very fast due to its use of FFT for computing of the integral sums.

We'll explain this using two-dimensional example, following [7], p.100. Let the integral is defined in a finite interval tex2html_wrap_inline1775 :


Replacing functions U(x) and tex2html_wrap_inline1779 with step functions tex2html_wrap_inline1781 and tex2html_wrap_inline1783 , defined in the sampling points of the grid with tex2html_wrap_inline1785 , and tex2html_wrap_inline1787 we convert the integral 5 to the form:


Taking integrals in 5 we obtain:


where: tex2html_wrap_inline1789 are analytically expressed with the help of Fresnel integrals, depending only onto the difference of indices. Sums tex2html_wrap_inline1791 can be easily calculated for all indices m as one convolution with the help of FFT.

Filter Fresnel, defined starting from version 1.1, implements this algorithm for two-dimensional diffraction integrals. It is almost as fast as forvard (still 2 to 5 times slower), it uses 8 times more memory (to perform the two-dimensional convolution) than forvard and it allows for ``more honest'' calculation of Fresnel diffraction. As it does not require any protection bands at the edges of the region, the model may be built in a smaller grid, therefore the resources consumed and time of execution are comparable or even better than that of forvard. Fresnel does not accept negative propagation distance. When possible Fresnel has to be used as the main propagation engine within LightPipes.

Warning: Fresnel does not produce valid results if the distance of propagation is comparable with (or less than) the characteristic size of the aperture, at which the field is diffracted. In this case forvard or steps should be used.

Direct integration

Direct calculation of the Fresnel-Kirchoff integrals is very inefficient in two-dimensional grids. The number of operations is proportional to tex2html_wrap_inline1795 , where N is the grid sampling. With direct integration we do not have any reflection at the grid boundary, so the size of the grid can just match the cross section of the field distribution. LightPipes include a program forward realizing direct integration. Forward has following features:

The direct calculation of diffraction integral in the Fresnel approximation is used, thus forward can not be used to propagate the field to a short distances. Use forvard of steps instead.

Forward is a very slow program.

Finite difference method.

It can be shown that the propagation of the field U in a medium with complex refractive coefficient A, is described by the differential equation:


To solve this equation, we re-write it as a system of finite difference equations:


Collecting terms we obtain the standard three-diagonal system of linear equations, solution of which describes the complex amplitude of the light field in the layer tex2html_wrap_inline1803 as a function of the field defined in the layer Z:


where (we put tex2html_wrap_inline1807 )




The three-diagonal system of linear equations (9) is solved by standard elimination (double sweep) method, described for example in [11]. This scheme is absolutely stable (this variant is explicit with respect to the index i and implicit with respect to the index j). One step of propagation is divided into two sub-steps: the first sub-step applies the described procedure to all rows of the matrix, the second sub-step changes the direction of elimination and the procedure is applied to all columns of the matrix.

The main advantage of this approach is the possibility to take into account uniformly diffraction, absorption (amplification) and refraction. For example, the model of a waveguide with complex three-dimensional distribution of refraction index and absorption coefficient (both are defined as real and imaginary components of the (three-dimensional in general) matrix tex2html_wrap_inline1813 ) can be built easily.

It works also much faster than all described previously algorithms on one step of propagation, though to obtain a good result at a considerable distance, many steps should be done. As the scheme is absolutely stable (at least for free-space propagation), there is no stability limitation on the step size in the direction Z. Large steps cause high-frequency errors, therefore the number of steps should be determined by trial (increase the number of steps in a probe model till the result stabilizes), especially for strong variations of refraction and absorption inside the propagation path.

Zero amplitude boundary conditions are commonly used for the described system. This, again, creates the problem of the wave reflection at the grid boundary. The influence of these reflections in many cases can be reduced by introducing an additional absorbing layer in the proximity of the boundary, with the absorption smoothly (to reduce the reflection at the absorption gradient) increasing towards the boundary.

In LightPipes version 1.2 the refraction term is not included into the propagation formulas, instead the phase of the field is modified at each step according to the distribution of the refractive coefficient. This ``zero-order'' approximation happened to be much more stable numerically than the direct inclusion of refraction terms into propagation formulas. It does not take into account the change of the wavelength in the medium, it does not model backscattering and reflections back on interfaces between different media. Perhaps there are other details to be mentioned.

The described algorithm is implemented in a filter steps. Steps accepts 6 arguments, first - the step size, second- the number of steps, third and fourth - files containing absorption and refraction coefficients to be imported and the last two arguments are the filename where field cross-sections are saved after each N steps, where N is the last argument. All filenames are optional.

The format of imported datafiles is the same as exported for using in gnuplot. The format of the datafile into which cross sections of the field are saved at every N-th step is similar to the file generated by cros_out. It has one more column, correspondent to Z coordinate. Cross sections correspondent to different values of Z are divided by an empty line.

For example:

begin 0.004 0.63e-6  | circ_ap 0.0015 | lens 0.5 |\
steps 0.01 75 void void steps_out > /dev/null 

steps: void refractive coefficient file, skipping
steps: void absorption coefficient file, skipping

produces file steps_out which can be used to plot the radial intensity distribution in propagating beam, for example with gnuplot:

set parametric
splot 'steps_out' using 1:6:2 with lines

Figure 5: Steps output for a simple lens.

Steps has built-in absorption layer along the grid boundaries (to prevent reflections), occupying 10% of grid from each side. Steps can import absorption/refraction coefficients defined in the format similar to exported by file_int. Steps is the only filter in LightPipes allowing for modeling of (three-dimensional) waveguide devices.

Like forvard, Steps can inversely propagate the field, for example the sequence ...|steps 0.1 1 | steps -0.1 1 |... doesn't change anything in the field distribution. Author has tested this reversibility also for propagation in absorptive/refractive media, examples will follow.

Steps implements scalar approximation, it is not applicable for modeling of waveguide devices in the vector approximation, where two components of the field should be taken into account.

Splitting and mixing beams

There are two commands in LightPipes, useful to model interferometers. With b_split we split the field structure (amplitude division) - one part will be written to the disk file, another will be passed to stdout. There is similar Unix command tee, writing two copies of any input - one to a file and the second to stdout. B_split works only with LightPipes data structures and provides the possibility of amplitude division, preserving the sum of intensities of the two splitted beams. The two obtained files could be processed separately and then mixed again with a program b_mix. Consider a simple model of an interferometer consisting of two apertures:

Figure 6: A simple interferometer

#/bin/csh -f
#model of a screen with two round pinholes
       begin 0.005 0.55e-6 512 > beam
       circ_ap 0.00012 -0.0005 < beam > beam1
       circ_ap 0.00012 0.0005  < beam > beam2
#Mixing two beams, rewriting initial beam:
        b_mix beam1 < beam2 > beam
#intensity in the plane of pinholes:
        cat beam | file_ps 128 >/dev/null
#propagation 0.75 m forward and second postscript file
        forvard 0.75 < beam  | file_ps 128  >/dev/null
#removing all temp  files
        rm beam beam1 beam2
#all done

In this script we have formed two beams (plane waves) each containing one ``shifted'' hole. After mixing these two beams we have a screen with two holes: a Young's interferometer. After execution of this script, two intensity distributions are stored in files and The result of the interference is shown in (Fig 7):

Figure 7: Intensity distributions in the plane of the screen and 75cm behind the screen

We can ``play'' with the interferometer model, moving the pinholes and changing their sizes and shapes. The following script models the interference of a plane wave diffracted at three round apertures:

#!/bin/csh -f
#model of a screen with three round holes

        begin 0.005 0.55e-6 512 > beam
#Three circular apertures:
       circ_ap 0.00012 -0.0005 -0.00025< beam > beam1
       circ_ap 0.00012 0.0005  -0.00025< beam > beam2
       circ_ap 0.00012 0.       0.00025< beam > beam3
#Mixing three beams, rewriting initial beam:
       b_mix beam1 < beam2 | b_mix beam3 > beam
#writing  Postscript Picture 
        cat beam | file_ps 128 >/dev/null
#propagation 0.75 m and second file 
        forvard 0.75 < beam  | file_ps 128  >/dev/null
#removing all beam files
        rm beam beam1 beam2 beam3
#all done

Figure 8: Intensity distributions in the plane of the screen and 75cm behind the screen

The next interferometer is more interesting:

#!/bin/csh -f

#model of a screen with two slits, one  tilted
        begin 0.005 0.55e-6 512 > beam

#Two slits one tilted :
       rect_ap 0.0001 0.0025  -0.0005 < beam > beam1
       rect_ap 0.0001 0.0025  0.0005 0 15 < beam > beam2
#Mixing two beams, rewriting the initial beam:
       b_mix beam1 < beam2  > beam
#writing  Postscript Picture of the screen with slits
        cat beam | file_ps 128 >/dev/null
#propagation 0.75 m and second file
        forvard 0.75 < beam  | file_ps 128  >/dev/null
#removing all beam files
        rm beam beam1 beam2 

#all done

Figure 9: Intensity distributions in the plane of the screen and 75cm behind the screen

In the last example the intensity distribution is modulated by the wave, reflected from the grid edge, though it gives a good impression about the general character of the interference pattern. To obtain a better result, the calculations can be conducted in a larger grid or with another numerical method.

This example uses direct integration algorithm (the input and output have different scales and samplings):

#!/bin/csh -f

#model of a screen with two slits, one of them tilted
#note  the number of points and the size of the grid 
        begin 0.0026 0.55e-6 64 > beam
#Two slits one tilted, the slits are as high as the grid:
       rect_ap 0.0001 0.0025  -0.0005 < beam > beam1
       rect_ap 0.0001 0.0025  0.0005 0 15 < beam > beam2
#Mixing two beams, rewriting the initial beam:
       b_mix beam1 < beam2  > beam
#writing  postscript file of the screen with slits
        cat beam | file_ps 64 >/dev/null
#propagation 0.75 m into different grid, writing  second file
        forward 0.75 0.005 64 < beam  | file_ps 64  >/dev/null
#removing all beam files
        rm beam beam1 beam2
This example uses approximately 25 times less memory than the previous FFT, calculation takes from minutes to tens of minutes, depending on the floating-point efficiency.

Figure 10: Intensity distributions in plane of the screen and 75 cm after the screen, note that input and output have different scales, input grid is 2.6x2.6mm, output is 5x5mm


Interpol is the tool to manipulate the size and the dimension of the grid and to change the shift, rotation and the scale of the field distribution. It accepts six command line arguments, the first is the new size of the grid, the second and others are optional. The second argument gives the new sampling, the third and fourth determine shifts in X and Y directions, the fifth gives the field rotation (after shifting). The last argument defines the magnification, its action is equivalent to passing the beam through afocal system with magnification M (without diffraction, but preserving the integral intensity). For example if the field was propagated with FFT algorithm forvard, then the grid contains empty borders which are not necessary for further processing. Other way around, after forward we have to add some empty borders to continue with forvard. Interpol is useful for interpolating into grids with different size and samplings. Of course, it is not wise to interpolate a grid of 512x512 points into a grid of 4x4, and then back to 512x512, because all information about the initial field will be lost. The same is true for interpolating a distribution defined in the square of 1x1m to 1x1mm and back. When interpolating into a grid with larger size, for example from 1x1 to 2x2, the program puts zeros into added regions.

In LightPipes1.2 a better (and slower) interpolation routine (12-point two-dimensional cubic interpolation) is introduced. Linear interpolation remains with the name interp1.

Examples of transformations:

The following example (Fig. 11) illustrates the usage of interpol for transition from a fine grid used by forvard (near field) to the coarse grid used by forward (far field).
#!/bin/csh -f
#starting with forvard (near field)
begin 0.01 1e-6 | rect_ap 0.005 | forvard 0.2 >foo
#continuing with forward (far field)
file_ps <foo |interpol 0.0075 64|forward 400 0.4>foo1
file_ps <foo1| file_int in |cros_out cros >/dev/null
rm foo foo1

Figure 11: Intensity distributions in the near field, in the far field, far field surface plot and far field cross section (from left to right)

Phase filters

There are four phase filters available in LightPipes - wavefront tilt, the quadratic phase corrector lens, a general aberration in the form of a Zernike polynomial, and a user defined filter. To illustrate the usage of these filters let's consider the following examples:

begin 0.04 1e-6   | rect_ap 0.02 | lens 8 | cros_out outpha  > out
file_ps  < out  | forvard 4 |file_ps >/dev/null

The first sequence of operators forms the initial structure, filters the field through the rectangular aperture and then filters the field through the positive lens with an optical power of 0.125D (the focal distance of 1 / 0.125=8m). With the second command we propagate the field 4 m forvard. As 4m is exactly the half of the focal distance, the cross section of the beam must be reduces twice.

Figure 12: The phase distribution after passing the lens, intensity in the plane of the lens and at the distance equal to the half of the focal distance. (from left to right)

We have to be very careful, propagating the field to the distance which is close to the focal distance of a positive lens- the near-focal intensity and phase distributions are localized in the central region of the grid occupying only few sampling intervals. This leads to the major loss of information about the field distribution. The problem is solved by applying the coordinate system which is tied to the divergent or convergent light beam. The tools to do this will be described later.

The lens may be decentered, lens 8 0.01 0.01 produces a lens with a focal distance of 1/0.125, shifted by 0.01 in X and Y directions. When the lens is shifted, only the phase mask correspondent to the lens is shifted, not the aperture of the lens or light beam itself.

The wavefront tilt is illustrated by following examples:

begin 0.04 1e-6   | rect_ap 0.02 | tilt 1e-4 -1e-4 |forvard 8  > foo
cros_out out <foo >/dev/null

The result is shown in Fig. 13.

Figure 13: The intensity and the phase after tilting the wavefront by 1e-4 rad and propagating it 8 m forward. Note the transversal shift of the intensity distribution and the phase tilt.

In this example the wavefront was tilted by tex2html_wrap_inline1827 rad in X and Y directions, then propagated to the distance Z=8m. In the output distribution we observe the transversal shift of the whole intensity distribution by tex2html_wrap_inline1831 .

Two arguments must be passed to tilt even if one of the two arguments of tilt equals to zero.

Zernike polynomials

Any aberration in a circle can be decomposed over a sum of Zernike polynomials. Formulas given in [2] have been directly implemented in LightPipes. The filter is named Zernike and accepts four command line arguments:

  1. The radial order n (first column in Table 13.2 [3] p. 465).
  2. The azimuthal order m, tex2html_wrap_inline1837 , polynomials with negative n are rotated tex2html_wrap_inline1697 relative to the polynomials with positive n. For example Zernike 5 3 1 1 gives the same aberration as Zernike 5 -3 1 1, but the last is rotated tex2html_wrap_inline1697 . This index corresponds to the n-2m given in the third column of Table 13.2 in [3], p. 465.
  3. The radius R
  4. The amplitude of aberration in radians, at R
We can uniformly introduce the lens and tilt using Zernike. The difference is that lens and tilt accept conventional meters and radians, which are widely in use for description of optical setups, while Zernike uses the amplitude of the aberration, which frequently has to be derived from the technical description.

Cylindric lens can be modeled as a combination of Zernike:

begin 0.01 1e-6 | circ_ap 0.0045 | file_ps\
| Zernike 2 2 0.0045 20 | Zernike 2 0 0.0045 -10\
| Fresnel 1.55 | interpol same same 0 0 45 |\
file_ps > /dev/null

Figure 14: The intensity in the input plane and after propagation through a cylindrical system, defined as a combination of Zernike polynomials and free space propagation.

The result of propagation to the focal plane of cylindric lens, defined as a combination of two Zernike polynomials, is shown in Fig. 14.

Spherical coordinates

The principle of beam propagation in the ``floating'' coordinate system (for the case of a lens waveguide) is shown in Fig. 15

Figure 15: Illustration for the light propagation in lens waveguides with fixed and floating coordinate systems.

The spherical coordinates follow the geometrical section of the divergent or convergent light beam. Propagation in spherical coordinates is implemented with programs lens_forvard and lens_fresn. Both filters accept two parameters: the focal distance of the lens placed in the input of the propagation distance, and the distance of propagation. When lens_forvard or lens_fresn is called, it ``bends'' the coordinate system so that it follows the divergent or convergent spherical wavefront, and then propagates the field to the distance Z in the transformed coordinates. Filter convert should be used to convert the field back to the rectangular coordinate system. Some LightPipes filters can not be applied to the field in spherical coordinates.

As the coordinates follow the geometrical section of the light beam, operator lens_forvard 10 10 will produce floating exception because the calculations just can not be conducted in a grid with zero size (that is so in the geometrical approximation of a focal point). On the other hand, diffraction to the focus is equivalent to the diffraction to the far field (infinity), thus FFT convolution algorithm will not work properly anyway. To model the diffraction into the focal point, a more complicated trick should be used:

begin 0.01 1e-6|rect_ap 0.005|lens 10|lens_fresn 1.111111 1 \
|convert y|cros_out cros>/dev/null

Here we calculate the diffraction to the focus of a lens with a focal distance of 1m. It is represented as a combination of a weak phase mask lens 10 and a ``strong'' geometrical coordinate transform lens_fresn 1.111111 1. The grid after propagation is tex2html_wrap_inline1855 10 times narrower than in the input plane. The focal intensity is 650 times higher than the input intensity and the wavefront is plain as expected. The same method is used in the script farfield (in the example's directory) to calculate the farfield distributions.

Diffraction behind the beam waist is more tricky: it can be done with lens_forvard but will not work with lens_fresn. Spherical coordinate variant of steps is not implemented, because the spherical coordinates also should be applied to absorption/refraction coefficients, which is far from trivial.

User defined phase and intensity filters

A filter called fil_ter is used to manipulate the phase and intensity of the light beam.

USAGE: fil_ter C1 C2 F [N], where
C1 is character constant with valid values 'int' and 'pha'
C2 is character constant with valid values 'mult' and 'subst'
F is the name of a file containing intensity/phase mask
It may be portable graymap or anymap (*.pgm, *.pnm) with a grayscale data.
Matrix F must have the same dimensions as the grid 
N is (any) optional argument indicating that the intensity mask should
be normalized before applying

fil_ter int mult aa: filter THROUGH  the intensity mask from a file aa
fil_ter pha subst aa: substitute the phase with a phase taken from file aa
fil_ter pha mult aa: filter the field through phase filter aa
fil_ter int subst aa: substitute intensity with one taken from file aa
fil_ter int subst aa haha: substitute the intensity with a normalized
one taken from the file aa

The distributions for importing may be formed by an external program - here is a simple Fortran77 example (arb.f) of such an operation:

Figure 16: An arbitrary intensity distribution

c     Copyright Gleb Vdovin, 1995
c     This Fortran program forms complex
c     intensity and phase filters in files in and pha, 
c      to use with LightPipes
      open(1, file='in', form='formatted')
      open(2, file='pha', form='formatted')

      do i=1,256
c     ii corresponds to the X coordinate

         do j=1,256
c     jj corresponds to the Y coordinate

c     ai and af are  some arbitrary function 
c     of the coordinates ii and jj,
c     ai corresponds to intensity and af is the phase:


c     with this "if" a circular aperture is cut:
c     writing intensity:
            if (ii**2+jj**2 .le. 64**2) then
            write(1,*) ai
            write(1,*) 0.
            end if

c     writing phase:
            write(2,*) af
         end do

c     empty line is a row separator:
       end do

c     all done:

This program forms two files: in and pha. A beam can be filtered through these filters to produce the intensity distribution shown in Fig. 16:

begin 0.1 1e-6 | fil_ter int subs in|file_ter pha subs pha\
 | file_ps > /dev/null

The phase and intensity distributions may be produced with a LightPipes script itself. File_int, file_pha and file_pgm provide such a possibility. The dimensions of the exported and imported distributions must be the same (this can be fixed with interpol), otherwise the behavior is not defined. Any (square) file (gif,tiff, bmp) containing graphics, can be converted into pgm format - for example using xv, and imported by fil_ter as the intensity or the phase distribution.

Random filters

random is a random filter with a uniform distribution. It uses standard C function rand (many consider it as a bad one, but it exists everywhere). By default random takes the random seed from the system clock, but if the result should be repeated with the same seed, the random seed (unsigned int) may be supplied as the third command line argument. The first argument may be int for the intensity or pha for the phase, the second gives the amplitude of the random phase screen (uniform distribution) in radians, or the amplitude of the intensity noise (uniformly distributed additive) normalized to the maximum intensity.

FFT and spatial filters

LightPipes.1.1 and later versions provide a possibility to perform arbitrary filtering in the Fourier space. There is an operator, performing the Fourier transform of the whole data structure: pip_fft.

pip_fft  performs FFT of the field distribution

USAGE:  pip_fft IND , where IND is direction,
IND=1 for forward, or -1 for inverse

Here is the example of a spatial filter:

#!/bin/csh -f

# random screen and propagation
begin 0.01 1e-6 200|circ_ap 0.003| random pha 3 | Fresnel 1 > foo

# how it looks like:
cat foo|file_ps 100 >/dev/null

# spatial filter
cat foo |pip_fft 1|circ_ap 0.0005|pip_fft -1 >foo1

# how it looks after filtering
cat foo1 |file_ps 100 >/dev/null

# remove FOO's
rm foo foo1

Figure 17: Intensity distributions before and after applying the spatial filter.

Note, we still can apply all the intensity and phase filters in the Fourier-space. The whole grid size in the angular frequency domain corresponds to tex2html_wrap_inline1857 where tex2html_wrap_inline1859 is the grid step. One step in the frequency domain corresponds to tex2html_wrap_inline1861 where x is the total size of the grid. LightPipes filters do not know about all these transformations, so the user should take care about setting the proper size (using the relations mentioned) of the filter (still in linear units) in the frequency domain.

In many cases, the frequency spectrum occupies only few sampling intervals. In this case filtering mask has to be defined only on a few grid points, making precision definition of the mask impossible (too few points). In this case, interpolation into a grid with maximum possible dimensionality is advised, to define the filter mask with a higher precision.

Service routines

There are two programs for phase unwrapping unf3 and unf4. The phase distribution saved in the gnuplot format may be filtered to obtain unwrapped phase plot. The operation of unwrapping is sometimes (very frequently) somewhat unstable, therefore two different programs are given. There is no guarantee that these filters will always work. They work well on smooth distributions. For example:

#!/bin/csh -f
begin 0.01 1e-6 60|lens 5 |file_pha pha >/dev/null
cat pha | unf3 > pha1
cat pha | unf4 > pha2
You can easily see the difference between pha, pha1 and pha2 after plotting the files.

Filter normal normalizes the integral intensity to 1.

Filter convert is used to transform spherical coordinates into Decarte's ones.

Filter absorber introduces intensity attenuator.

Filter a_phase compares two phase distributions, prints some info to the stderr and outputs the phase difference to stdout. a_phase is optional and may be not included into this distribution.

There are more programs like scale, c_matlab or c_xplot, operating on only intensity or only phase distributions. All these programs are more or less self-documented. They are distributed with source, therefore no explicit description is given here.

Laser amplifiers

The l_amplif introduces a simple single-layer model of a laser amplifier.

l_amplif filters the field through a laser medium

USAGE: l_amplif G L I, where G is the gain,
L and I are the length and the saturation intensity of the medium.
It is useful for modeling of amplifiers and lasers with active medium.

Diagnostics: Strehl, the output power and the center of gravity

Strehl prints to the stderr the current parameters of beam. Using it in in combination with l_amplif you may calculate and print out the output power of a laser with active medium.


begin 1 1 | rect_ap 0.5 | normal y | l_amplif 0 1 1 | Strehl y > /dev/null
normal: Normalization coefficient= 2.539215e-01
Strehl: ratio= 1.000000e+00 energy= 1.000000e+00
Center_of_gravity: x= 9.926167e-24 y= 4.268252e-22
Grid size: 1.000000e+00, Grid sampling: 256  

begin 1 1 | rect_ap 0.5 | normal y | l_amplif 1 1 1 | Strehl y > /dev/null
normal: Normalization coefficient= 2.539215e-01
Strehl: ratio= 1.000000e+00 energy= 1.499308e+00
Center_of_gravity: x= 1.875808e-23 y= -3.237424e-21
Grid size: 1.000000e+00, Grid sampling: 256

Example models

Shearing interferometer

We can easily model a shearing interferometer using the operators introduced in the previous sections. Suppose we have to model a device shown in Fig. 18.

Figure 18: The shearing interferometer

Let tex2html_wrap_inline1865 , tex2html_wrap_inline1867 , tex2html_wrap_inline1869 , the focal length of the lens is -20m, the beam is divided 1:1 and one of two beams is shifted by D=3mm in the X direction and tex2html_wrap_inline1873 in the Y direction after reflection on the back surface of the plane-parallel plate.

#!/bin/csh -f
# starting, lens, Z1 
begin 0.04 5e-7 | circ_ap 0.01 |lens -20|forvard 0.5|b_split foo1 > foo
# shifting beam 3 mm without real interpolation, mixing it
# with centered beam
interpol 0.04 256 0.003 0.001 <foo1|b_mix foo |forvard 0.5>foo2
# postscript output
file_ps 128<foo2|cros_out out >/dev/null
rm foo foo1 foo2
In Fig 19 is shown the intensity distribution in the plane of the screen, obtained as a result of running this script:

Figure 19: The intensity distribution in the plane of the screen, obtained with the model of a shearing interferometer, the intensity cross-section is shown in the right graph

We can put a source of arbitrary aberration on the place of a lens, for example we can introduce a spherical aberration with amplitude of 10 Rad:

#!/bin/csh -f
# starting, lens, Z1 
begin 0.04 5e-7 | circ_ap 0.01 |Zernike 4 0 0.01 10|\
forvard 0.5|b_split foo1 > foo
# shifting beam 3 mm without real interpolation, mixing it
# with centered beam
interpol 0.04 256 0.003 0.00 <foo1|b_mix foo |forvard 0.5>foo2
# postscript output
file_ps 128<foo2 >/dev/null
rm foo foo1 foo2
Changing the amplitude and the order of aberration we can obtain all the shearing interferograms shown in the section about shearing interferometers of [3].

Figure 20: Shearing interferograms of the spherical aberration, coma and astigmatism

Rotational shear interferometer

In the rotational shear interferometer the beam interferes with a copy of itself rotated by the angle tex2html_wrap_inline1875 around optical axis. Interferometer is useful for detecting of asymmetrical aberrations. We shall consider only bare-bone (no propagation and diffraction) model of such an interferometer:

#!/bin/csh -f
# Rotational interferometer with  coma
begin 0.04 5e-7 | circ_ap 0.01 |Zernike 3 1 0.01 10|b_split foo1 > foo
# rotating and mixing
interpol 0.04 256 0. 0. 180 <foo1|b_mix foo >foo2
# postscript output
file_ps 128<foo2 >/dev/null
rm foo foo1 foo2
#!/bin/csh -f
# Rotational interferometer with  astigmatism
begin 0.04 5e-7 | circ_ap 0.01 |Zernike 2 2 0.01 10|b_split foo1 > foo
# rotating and mixing
interpol 0.04 256 0. 0. 90 <foo1|b_mix foo >foo2
# postscript output
file_ps 128<foo2 >/dev/null
rm foo foo1 foo2
#!/bin/csh -f
# Rotational interferometer with  high order aberration
begin 0.04 5e-7 | circ_ap 0.01 |Zernike 7 3 0.01 10|b_split foo1 > foo
# rotating and mixing
interpol 0.04 256 0. 0. 90 <foo1|b_mix foo >foo2
# postscript output
file_ps 128<foo2 >/dev/null
rm foo foo1 foo2

Figure 21: Rotational interferograms of coma tex2html_wrap_inline1695 , astigmatism tex2html_wrap_inline1697 and the high order Zernike aberration tex2html_wrap_inline1697

Radial shearing interferometer

In radial shearing interferometer the beam interferes with a copy of itself magnified by the factor M. Interferometer is useful for detecting of axisymetrical aberrations. As for the previous case, we shall consider only bare-bone (no propagation and diffraction) model of such an interferometer:

Figure 22: Radial shearing interferograms of defocus, spherical aberration and the higher-order aberration

#!/bin/csh -f
# Radial shearing  interferometers:

begin 0.04 5e-7 | circ_ap 0.01 |Zernike 2 0 0.01 10|b_split foo1 > foo
# Radial magnification
interpol 0.04 256 0. 0. 0 1.3 <foo1|b_mix foo >foo2
# postscript output
file_ps 128<foo2 >/dev/null
rm foo foo1 foo2

#spherical aberration
begin 0.04 5e-7 | circ_ap 0.01 |Zernike 4 0 0.01 10|b_split foo1 > foo
# Radial magnification
interpol 0.04 256 0. 0. 0 1.3 <foo1|b_mix foo >foo2
# postscript output
file_ps 128<foo2 >/dev/null
rm foo foo1 foo2

#High order  aberration without axial symmetry
begin 0.04 5e-7 | circ_ap 0.01 |Zernike 10 4 0.01 10|b_split foo1 > foo
# Radial magnification
interpol 0.04 256 0. 0. 0 1.3 <foo1|b_mix foo >foo2
# postscript output
file_ps 128<foo2 >/dev/null
rm foo foo1 foo2

Twyman-Green interferometer

Figure 23: Twyman-Green interferometer (left) and the output of the interferometer model

The scheme of the interferometer is shown in Fig. 23. In this example we put D=1cm, tex2html_wrap_inline1865 , tex2html_wrap_inline1889 , tex2html_wrap_inline1891 , tex2html_wrap_inline1893 , tex2html_wrap_inline1869 . The first mirror is plane and the second is aberrated (coma) with aberration amplitude of tex2html_wrap_inline1703 . The beam-splitter is ideal dividing the beam 3:7. We run the following script:

#!/bin/csh -f
# propagation and splitting
begin 0.03 5e-7 | circ_ap 0.005 | forvard 0.5 | b_split foo 0.3 >foo1
# Propagation to mirror1, reflection,
# propagation back and beam-splitter again
forvard 1 <foo|i_filter 0.7 >foo2
# Propagation to mirror2, reflection,
# propagation back and beam-splitter again
forvard 0.4 < foo1 | Zernike 3 1 0.005 25|forvard 0.4|i_filter 0.3 >foo3
# mixing the two beams and propagation to the screen
b_mix foo2 <foo3 |forvard 1| interpol 0.012 128 | file_ps 128>/dev/null
# deleting the temp files
rm foo foo1 foo2 foo3

Unstable laser resonator

We know enough to solve a real problem with LightPipes. Here is the script which calculates the output mode of an unstable resonator with square mirrors. The example is taken from [10]. The scheme of the resonator is shown in Fig. 24 together with the equivalent lens waveguide scheme used in calculations. The infinite sequence of reflections from resonator mirrors is replaced here by finite but long enough sequence of propagations through afocal magnifying lens system.

Figure 24: An unstable confocal resonator and a lens waveguide

Figure 25: Intensity and phase distributions at the output of the model of unstable resonator

#!/bin/csh -f

#	LightPipes script: a model of an unstable resonator
#	with N=10, M=2, 
#	See D.B. Rench, Applied  Optics 13, 2546...2561 (1974)
#	All dimensions are taken from the table on p. 2552
#	in this reference

#	theoretically the initial field distribution may be random 
#	but we'll use a good plane wave because it won't
#	converge fast with a random distribution
begin 20. 3e-4 512 >field1

#	20 field iteration inside the resonator
@ i=0
while( $i < 20)
#	reflection from the convex mirror and propagation 
	rect_ap 5.48 < field1 | lens -1e4 |forvard 1e4 > field2
#	reflection from the concave mirror and propagation
	rect_ap 10.96 <field2 | lens 2e-4  | forvard 1e4 | normal y > field1
	Strehl y < field1 >/dev/null
	@ i += 1
# output of the resonator is screened by the output mirror:
# the output 
rect_screen 5.48 <field1 | file_int in |\
 file_pha pha|file_ps |Strehl y>/dev/null
# removing the temp files
rm field1 field2

When executed, the script prints values of the feedback coefficient and other information for the output beam at each iteration.

From the output sequence we see that after 20 iterations the feedback coefficient is quite stable, so the calculations can be terminated. If the files field1 and field2 are not deleted, we can continue the calculations later. To do this, the operator begin should be commented, othervise it will start again the calculations from a plane wave.

We can compare the intensity and phase distribution in the output mode with original distributions obtained by Rench [10]:

In this example we calculated the mode profile of an ideal unstable resonator. A resonator with aberrated and misaligned mirrors can be also easily modeled. One problem with the considered model is the way, the mirror's curvature is defined. Curved mirrors are modeled as phase masks lens and for high curvatures the phase difference between the neighbour points of grid may be higer than tex2html_wrap_inline1899 , especially in the proximity of the mirror edges. If the phase difference is greater than tex2html_wrap_inline1899 , the results are not valid, so the nuber of points must be increased to fix the mentioned phase condition. It leads to enormous grids (we used 512x512 ) for a simple problem. This can be solved by applying spherical coordinate system to propagate spherical waves. In such a system the coordinate surfaces coincide approximately with the propagating wavefront.

Laser with unstable resonator in spherical coordinates.

The script for a laser with a confocal unstable resonator in spherical coordinates follows (including laser medium):

#!/bin/csh -f

#all dimensions are in mm here

begin 7. 3e-4 100 >field1

#	40 iterations inside the resonator
@ i=0
while( $i < 20)
#	reflection from the convex mirror and propagation 
	rect_ap 5.48 < field1|l_amplif 1e-4 1e4 1 | lens_fresn -1e4  1e4 > field2
#	reflection from the concave mirror and propagation
	rect_ap 10.96 <field2 |l_amplif 1e-4 1e4 1| lens_fresn 2e4  1e4\
|Strehl y > field1
	cp field1 field_out
	interpol 7. <field1 > field2 
	mv field2 field1
	@ i += 1
# output of the resonator is screened by the output mirror:
# the output 
echo 'Output beam:'
cat field_out |convert y |rect_screen 5.48 | file_ps\
 |Strehl y >/dev/null
# removing the temp files
rm field1 field2

Figure 26: The output (intensity and phase) of an unstable resonator, modeled in spherical coordinates

Here is the tail of the screen output:

Strehl: ratio= 8.415976e-01 energy= 3.978092e+02
Center_of_gravity: x= -1.330218e-03 y= -1.100146e-03
Grid size: 1.400000e+01, Grid sampling: 100
Output beam:
Strehl: ratio= 8.209055e-01 energy= 2.917869e+02
Center_of_gravity: x= -1.716706e-03 y= -1.600309e-03
Grid size: 1.400000e+01, Grid sampling: 100

The last script is much faster and the field occupies a larger part of the grid. It produces almost the same feedback coefficient, the Strehl number is also close to that obtained in the previous example, the difference is due to the better computational approach used in the latter case.

Propagation in lens-like/absorptive medium

Here the model of the propagation of a Gaussian beam in a lens-like waveguide is built. The profile of the refractive index is chosen so, that the beam preserves approximately its diameter in the waveguide (we use fundamental mode). We'll consider the propagation of an axial mode, tilted with respect to the waveguide axis and a non-axial mode. First, we form the distribution of a refractive coefficient in a file. We do that with a Fortran program (any other way, for example by hand ;-), can be used:

      program med
c     The program to form a quadratic 
c     distribution of the refractive index;
c     nc is the grid  sampling
c     an0 is the axis value of refractive index
c     an1 is n1 in the formula for refractive index: 
c     n=sqrt(an0^2-an0*an1*r^2) where 
c     r is the radial coordinate
c     xmax is the grid size

      read*, nc, an0, an1, xmax
      do i=1,nc
         do j=1,nc
            print*, sqrt(an0*an0 -an1*an0*r2)
            end do
            end do

Here we use the approximation for the profile of the refractive coefficient in th form: tex2html_wrap_inline1903 . It is a well-known fact [12] that the half-width of the fundamental Gaussian mode of a lens-like waveguide is defined as: tex2html_wrap_inline1905 . For a waveguide of 1x1mm, tex2html_wrap_inline1907 , tex2html_wrap_inline1909 and tex2html_wrap_inline1911 , the Gaussian mode has a diameter of 226 tex2html_wrap_inline1913 m. We shall use this value in the following script, which demonstrates three-dimensional propagation of tilted and off-axis beams in square absorptive medium:

#!/bin/csh -f

#	We form here the distribution of the refractive index
#	in the file refr.ind;
#       med is a Fortran program to do that:

f77 -o med  med.f -lm 
echo 100 1.5 400 1e-3  | med > refr.ind

#	The array of refractive indexes (1.5 ) is used here
#	to form an array of absorbtion coefficients.
#	No optics here, it's just a reuse of the table of numbers!

scale -1. < refr.ind >

# 	propagation of a  gauss beam in a quadratic medium: 
#	the steps back are implemented with the 
#	second operator to show the reversibility of the
#	operator:

begin 1e-3 1e-6 100 |  gauss 113e-6 | tilt 1e-3 0 |\
steps 1e-3 1000  refr.ind s.out 10 |\
steps -1e-3 1000  refr.ind si.out 10  > /dev/null

#	The graphs can be plotted with gnuplot:

#	set par
#	set con
#	splot 's.out' using 1:6:2 w l

# here is the  propagation of a non-axial gauss beam:

begin 1e-3 1e-6 100 |  gauss 113e-6 2e-4 |\
steps 1e-3 1000  refr.ind s_s.out 10 > /dev/null

# here we have two non-axial beams:

begin 1e-3 1e-6 100 |  gauss 113e-6 2e-4 >foo1
begin 1e-3 1e-6 100 |  gauss 113e-6 -2e-4 >foo2
b_mix foo2 < foo1 |\
steps 1e-3 1000  refr.ind s_ss.out 10 > /dev/null

rm foo*

The array of refractive coefficients is plotted in Fig. 27.

Figure 27: Distribution of the refractive index in a lens-like waveguide.

The result of the first example - centered Gaussian beam under small angle to the axis in a lens-like medium is shown in Fig. 28

Figure 28: Radial (x at the left, y at the right) energy distribution in a tilted Gaussian beam propagating in a lens-like absorptive medium.

The result of ``reversed'' propagation is shown in Fig. 29. Comparison with Fig. 28 proves the validity of algorithm for refractive/absorptive media.

Figure 29: Radial (x at the left, y at the right) energy distribution in a tilted Gaussian beam, reversely propagated in a lens-like absorptive medium.

The result of propagation of ``shifted'' Gaussian beam, incident parallel to the axis is shown in Fig. 30.

Figure 30: Radial (x at the left, y at the right) energy distribution in a shifted Gaussian beam, incident parallel to the waveguide axis.

The result of propagation of two ``shifted'' Gaussian beams, both incident parallel to the axis of the waveguide is shown in Fig. 31.

Figure 31: Radial (x at the left, y at the right) energy distribution in two Gaussian beams, incident parallel to the waveguide axis.

Steps is a ``rich'' filter. Whole classes of waveguide devices can be simulated. Simple examples in this manual illustrate only the very basic capabilities.

Inverse problem: reconstructing the phase from measured intensities

The idea is very simple: we have two near field intensity distributions, but we don't know the phase of the field. To reconstruct the phase we start from a field with uniform (random) distribution and iterate then the field ``forward'' and ``back'', substituting each time the measured intensity and leaving the phase untouched. This process converges to a field distribution with reconstructed phase, providing a good match to the measured intensities. This very specific algorithm may be easily implemented with LightPipes script:

Figure 32: Given intensity distributions (top), reconstructed intensities (middle) and reconstructed phase

#!/bin/csh -f # (C) Gleb Vdovin 1993-1995 
# This file is a part of
LightPipes.1.1 distribution 
# this script demonstrates the possibility of 
# phase reconstruction from two near-field intensity
# distributions i_1 and i_2
@ num_point1 = 128
@ num_point2=256
@ num_iter= 100
@ dist = 2

# here we start from a plane phase:
begin  0.011 6.33e-7 $num_point1| random 0 > f_far1
#begin  0.02 6.33e-7 $num_point| random 9 > f_near1

#	num_iter  iterations 
@ i=0
while( $i < $num_iter)
	echo $i
#	Here we substitute the known second intensity and interpolate
#	into a bigger grid - to make "forvard" happy
	fil_ter int subs i_2 < f_far1 |interpol 0.02 256>foo
#	propagation "back" to the plane of the first image
	cat foo| forvard -2|interpol 1.1e-2 128 > f_near1

#	substitution of the first intensity distribution
#	and propagation forward:
	fil_ter int subs  i_1 < f_near1 | Fresnel 2 > f_far1
	@ i += 1

Example script inverse solves similar problem for Young interferometer.

Optical information processing

Figure 33: Initial transparant (left), filtered with a high-frequency component cut (center), the same filter with low frequency components cut (right).

LightPipes can be used to model the classical setups of Fourier-optics information processing. For example, this script models coherent optical computer, performing Fourier-plane image filtering:

#!/bin/csh -f

# 	importing the intensity mask:
begin 1e-2 5e-7 200 | fil_ter int subs reg.pgm | forvard 1 > foo
#	lens, focal filter, free space...
lens 1 < foo | forvard 1 | circ_ap 0.0015 | forvard 1 > foo1
#	lens free space, output.
lens 1 < foo1 | forvard 1 | file_ps 200 4 > /dev/null
#	all done

compare the result with the result of:

lens 1 < foo | forvard 1 | circ_screen 0.0015 | forvard 1 > foo1
#	lens free space, output.
lens 1 < foo1 | forvard 1 | file_ps 200 4 > /dev/null
#	all done

The initial distribution (file reg.pgm) and two filtered images are shown in Fig. 33.

Generation and reconstruction of interferograms

In the simplest case the interferogram (hologram) is generated by mixing of an aberrated beam with a plane wave:

#!/bin/csh -f
#	Aberrated beam:
begin 1e-2 1e-6 200 | Zernike 7 1 0.72e-2 10 > foo
#	Saving and unfolding the phase:
file_pha pha 50 < foo >/dev/null
unf3 < pha > pha1  
#	Addition of tilt, to obtain the vertical fringes:
tilt 2e-3 0 < foo > foo1
#	Ideal beam, mixing -> generation of interferogram:
begin 1e-2 1e-6 200 | b_mix foo1 | file_int int_in 256 > foo2
#	Plot the interferogram:
file_ps 200 1 < foo2 > /dev/null
#	All done
The generated interferogram, saved in file may be directly plotted, as shown in Fig. 34.

Figure 34: The phase distribution and the interferogram, generated with LightPipes script.

To reconstruct the phase from the obtained interferogram, we shall build a model of the algorithm of phase reconstruction, using shift and filtering in the fourier domain. Let the fringe pattern has the form:




where the phase tex2html_wrap_inline1915 contains information about the mirror shape and the term tex2html_wrap_inline1917 describes the wavefront tilt. Fourier transform of the expression (13) gives:


where the capitals A and C denote the Fourier spectra and f is the spatial frequency.

Figure 35: Fourier spectra of a fringe pattern modulated by the spatial frequency tex2html_wrap_inline1701 (top) and a single spectrum translated to the origin (bottom).

Figure: Spectrum of the interferogram shown in Fig. 34, with a side lobe shifted to the origin (left), the same after filtering (middle) and finally, the reconstructed phase (right).

We may take one of two side-lobe spectra tex2html_wrap_inline1927 or tex2html_wrap_inline1929 and translate it to the origin with zero spatial frequency as shown in Fig. 35. Now we can take the inverse Fourier transform of the translated spectrum to obtain c(x,y) defined in the expression(13, 14). Calculating the complex logarithm of the expression (14) we obtain the phase tex2html_wrap_inline1915 :


The phase tex2html_wrap_inline1915 is obtained indeterminate to a factor of tex2html_wrap_inline1937 with its principal value localized in the range of tex2html_wrap_inline1939 . To obtain the continuous phase map, a special unwrapping algorithm, removing discontinuities with an amplitude close to tex2html_wrap_inline1937 is applied to the reconstructed phase map. The algorithm implementation includes the following steps:

  1. Calculation of the spectrum of the initial intensity distribution pip_fft 1. It consists of a central lobe and two side-lobes, containing information about the phase.
  2. Shift of the obtained spectrum to put one of the side lobes into the origin: interpol.
  3. Filtering, leaving only the side lobe, containing th phase information. This operation is performed as a product with rectangular spectral window rect_ap.
  4. Reversed Fourier transform pip_fft -1 and saving of the wrapped phase.
  5. Phase unwrapping.
Described procedure is performed by the following script:
#!/bin/csh -f
#Importing the interferogram intensity pattern:
begin 1e-2 1e-6 200 | fil_ter int subs int_in >foo
#Fourier transform and shift in the spectral domain: 
pip_fft 1 < foo | interpol 1e-2 200  10.e-4 0 > foo1
#Postscript of the intensity distribution in the Fourier domain:
interpol 4e-3 < foo1 | file_ps 200 14 |rect_ap 15e-4 1\
|file_ps 200 14> /dev/null

#Filtering in the Fourier domain, reverse FFT and saving of the phase:
rect_ap 15e-4 1 < foo1 |pip_fft -1|file_pha pha 50 > foo
#Unwrapping the phase:
unf3 -1 < pha > pha1
#All done
The initial interferogram is even (symmetric) with respect to the vertical axis, but the phase is odd, thus the phase is reconstructed with inverted sign (this depends on the direction of the shift in the frequency domain). The amount of shift in the frequency domain can be calculated using relations, given in the section containing description of pip_fft, it also can be determined experimentally by visualizing the spectrum and shifting it, to position one of its side lobes at the origin. This can be done with interpol. Visualization of the spectrum side lobes requires high tex2html_wrap_inline1721 (of the order of 10).

LightPipes: port to MSDOS

MSDOS port has the full functionality of the Unix distribution. All command names are cut to 8 letters, capitals are not supported. Thus, for example, instead of a Unix command LongCommandOfPipes you must use the abridged version: longcomm. Simple NULL have to be used instead of /dev/null The .EXE files are compiled with DJGPP compiler (version 1.12) - many thanks to DJ Delorie. DJGPP is the MSDOS port of Gnu C compiler, it is freely available from any Simtel mirror. If you don't have enough RAM, the program will swap to the harddisk. Under MSDOS it was approximately 2 times slower than under Linux on the same computer (Pentium90 of unknown origin with 24MB of RAM). I tested it on a 386 with 3MB of RAM- it runs all the example scripts and some of them take number of days to execute. It runs happily in a DOS box under MS Windows 3.1. It also runs all scripts under DOSEMU 0.60.1 with 8MB of DPMI (I tried it again under DOSEMU 0.60.4. - it works but it was not really stable).

Use the part 1 of this manual as the command reference and tutorial for MSDOS port. You just have to take into account some basic differences between DOS and Unix.

The syntax depends on the version of you use. Microsoft is supposed to be present. 4DOS is better. There is a problem with programming loops in script batch files, the simplest solution is just a mechanical repetition of the command sequence inside the batch file. The script files are getting longer, but the result is achieved. There is a full set of ported to MSDOS example scripts in the directory scripts.

Some rules to port a script from Unix to DOS (the list is not complete):

Command reference

The command reference is excluded from this manual. If you have a working program, just type the command you want to know about at the shell prompt, it will print short explanation to the screen.

If you want to get a list of all descriptions, the following command (under Unix Csh) will do the job:

rm refer; touch refer;
find . -perm -555 -exec {} >> & refer \;

The resulting file refer contains all descriptions of all executables.


G. Anderson, P. Anderson, The Unix C shell guide, Prentice-Hall (1986)

M. Born, E. Wolf, Principles of Optics, Pergamon, 464, 767-772 (1993)

D. Malacara, Optical shop testing, J.Wiley & Sons, (1992)

J.W. Goodman, Introduction to Fourier Optics, McGraw-Hill, 49-56 (1968)

W.H. Southwell, J. Opt. Soc. Am., 71, 7-14 (1981)

A.E. Siegman , E.A. Sziklas, Mode calculations in unstable resonators with flowing saturable gain: 2. Fast Fourier transform method, Applied Optics 14, 1874-1889 (1975)

N.N. Elkin, A.P. Napartovich, in "Applied optics of Lasers", Moscow TsniiAtomInform, 66 (1989) (in Russian)

R. Baracat, in "The Computer in Optical Research", ed. B.R. Frieden, Springer-Verlag, 72-75 (1980)

A. E. Siegman, Lasers, Oxford University Press (1986)

D.B. Rench, Three- dimensional Unstable Resonator Calculations with Laser Medium, Appl. Opt. 13, 2546-2561 (1974)

A.A. Samarskii, E.S. Nikolaev, Numerical Methods for Grid Equations, V.1, Direct methods, pp. 61-65, Birkhäuser Verlag 1989.

Dietrich Marcuse, Light Transmission Optics, Van Nostrand Reinhold, 267-280, 1972.

Gleb Vdovin
Mon Jan 13 16:55:15 MET 1999