OptiCass — a Cassegrain Radio Telescope Optimization and Design Tool

General introduction
Parameter description
    Input parameters
    Output parameters        
    Block diagram
    Cassegrain geometry
    Ray tracing
    Amplitude distribution
    Power pattern and losses
    Closed formulae

Tests and first results
Notes on new versions

Abstract: A software for analysis of optical properties of Cassegrain radio telescopes with offset primary is detaily described. An executable for Windows (9x, XP) environment is made available for download. First results of its application to symmetrical (non-offset) designs suggest that the widely practiced Ruze approximation for gain losses may be less universal than thought of till now.

General introduction

OptiCass is a program developed in Fortran as a tool to study properties of two-mirror Cassegrain type radio telescopes, especially these with offset primary mirror for which well known geometrical optics approximations may prove inadequate. Its name stands for OPTIcs or OPTImization of CASSegrain telescopes. Here the word optimization is used to mean mainly (but not exclusively) minimization of aperture efficiency losses through varying certain parameters.

OptiCass basic working relies on the ray tracing method and the Fourier transformation (DFT) of so obtained amplitude distribution on the telescope aperture. Unlike the sophisticated telescope design programs found on the market, OptiCass never uses paraxial technique — its ray tracing is exact for all rays also for axially symmetric configurations. As presented here it was developed during the first quarter of 2005 as completely original program. No borrowed procedures are involved except a graphics library for writing output plots in PostScript, the PSPLOT, which is free.

The program accepts some two dozens of input parameters and returns a number of derived parameters describing geometry and optical properties of a telescope corresponding to the current setup. User is given possibility to modify input parameters until the desired properties are arrived at. During this process the user may save intermediate results in text or PostScript files or copy some of them directly from the screen. Initial input parameters can be read from a parameter file or taken from an earlier saved session.


Program uses number of its 'standard' files. These have the following extensions:
— .par input file; contains settable parameters in the form displayed on OptiCass screen,
— .txt output textual summary written automatically just before session conclusion,
— .ps PostSript output plots of power pattern, its distortions and outline of mirrors with marks of traced rays where they intersect the aperture plane; also written automatically just before session conclusion,
— .dat table of coordinates, amplitudes and phases of all rays; written upon user request,
— .pat table of power pattern; written upon request,
— .res table of residual power pattern (after reference pattern subtraction); written upon request.
Besides these, the program upon request writes a temporary file with all the current parameter values to be possibly read-in at start of another session (the $ option).

Filenames of .par files read at the start are limited to 8 characters, but those read through the command line argument (drag'n'dropped) haven't any restriction as to the name length or extension type. Other files of the above list are given the .par file name. In case of the $ and ENTER input, or of error in the attempt to open a .par file requested (whereupon default parameters are assumed) OptiCass gives names for the outputs in the form RTddd or RTdd_, where dd(d) are two (or three) digits of the telescope diameter, D, in metres. E.g., RT05_ is for D = 5 m and RT200 for D = 200 m.


The program is invoked in typical Windows style, including the drag and drop manner. Unless a parameter file is dropped (in which case the first prompt is skipped), OptiCas first screen looks like this:

Fig. 1: OptiCass opening window. It contains a brief summary of commands and prompts
for specifying the source of input parameters. One may just press ENTER
and modify the defaults assumed later, on next prompt.

The user is given a few options to chose from. One of six preset files can be read in (it must be placed in the parent directory) or a user prepared file name typed-in following a * character (the name in this case can be 8 character long at most, and must be typed without extension which is assumed to be .par). Parameters can be read also from most recently saved parameters (the $ option) or they may be made to assume values preset internally by the program (this happens also if a requested par file does not exist).

Of the six .par files known to OptiCass the last one may be useful for checking whether OptiCass works properly in a new environment by comparing TestO.txt and TestO.ps files, generated through simple run of the program, typing-in "6" on the first prompt and then "q" on the second (or just dropping TestO.par onto OptiCass.exe and typing "q"), with existing files named TestOptiCass.txt and TestOptiCass.ps.

Having successfully read 25 to 30 input parameters the program computes the remaining parameters (about 50 in number) and displays all of them in a compact table. The second screen that will appear may look like this:

Fig. 2: Display of read in (SETTABLE) and derived (COMPUTED) parameters. Such screen can be copied in its entirety into a new .par file. On input OptiCass will search for a line with exactly "       SETTABLE" (7 initial spaces are significant) at the beginning and then will read parameter values from the following 6 lines (ignoring all the content until the column marked #, inclusive). Note that here on first prompt at the table bottom a user has already typed number 110 and got corresponding response, 6 significant digit value of the parameter located in row 11 and column marked #0, and short explanation of the parameter meaning and its units.

At this point the user may change any of the settable parameters, one at a time or a few consecutive parameters at once, each time the updated SETTABLE table is displayed. Upon just ENTER, if any of the parameters has earlier been modified, the COMPUTED table is recomputed and displayed. An example follows.

Fig. 3: An example of parameter modification. Here a user replaced current value of parameter No 21 (feed offset in y coordinate) with 0.05 (metres are assumed in this case), then pressed sole ENTER key to get new values for the COMPUTED parameters.

Parameter description

All the displayed parameters are organized in the table of 16 rows and 5 columns. Each row is preceded by textual matter meant to indicate in a very abbreviated form the content of following five columns of decimal numbers. The five columns are numbered on their left. This row number (# = 0 to 15) in composition with colum number (#0 to #5) serves to refer by the user to his chosen parameter while working with OptiCass.

The parameter units are only metres, degrees, decibells, percents or none for dimensionless quantities. They are also indicated in the textual body, although not uniquely, but this aspect should not cause any problem for one who knows the parameter meaning. In case of doubt the user may issue the parameter number alone followed by ENTER key to get 6 significant digits of the value together with brief info and units (see bottom part of Fig. 2).

Input parameters
There are 5 rows of fully SETTABLE parameters in the displayed table. The 6th row (unnumbered) contains two conditionally settable parameters. Brief explanation of the content of successive rows follows.

Cassegrain: D,X_min,f,Xs_max,h2   [m]  0 
These 5 parameters are most basic and specify fully the optical geometry of the Cassegrain telescope without any imperfections or offsets. The parameters are:
D — diameter of circular aperture of the primary parabolic mirror.
X_min — dish rim (inner tip of offset paraboloid) separation from paraboloid axis (x-coordinate); in case of symmetrical telescope it assumes the value of –D/2.
f — focal distance of parent paraboloidal dish.
Xs_max — distance of outermost point of the subreflector edge from the axis. For symmetrical telescope it is just the radius of subreflector.
h2 — height of secondary focus above paraboloid vertex. It is defined to be positive towards the subreflector.

The next two rows contain 10 offsets, 5 of the subreflector and 5 of feeds. The first three in each case are linear translations in the x, y, and z direction, in this order. These Cartesian coordinates have the origin in the prime focus. The z-coordinate is directed along telescope axis (the common axis of paraboloid and hyperboloid) away from the mirrors, x- and y-coordinates lie in plane perpendicular to the z-axis. For the offset paraboloid, x is directed toward the point of dish that is furthest from the z-axis and y forms the right-hand triple. For symmetrical telescope the direction of x-axis is undefined and can be chosen freely by the user. With horizontal telescope mount however it may be practical to think of it as lying in the vertical plane perpendicular to the telescope altitude axis.

Subreflector: x,y,z,  tx, ty  [m,deg]  1
Feed offsets: x,y,z, Dtx,Dty  [m,deg]  2
x,y,z — offsets (translation or just coordinates) of the secondary mirror first focus (the one which in perfect alignment coincides with the prime focus); in case of the feeds these are its displacements from nominal place, the secondary focus (z measured positive toward the subreflector).
tx, ty — tilts of the secondary in the x-z plane and y-z plane, respectively. These angles are defined positive when a rotation moves the mirror vertex (or the normal, if the pivot is defined to lie at the vertex) towards positive x and y, respectively. The rotation by the ty angle is ALWAYS made about the prime focus, and by the tx angle about a center defined by the pivot parameter.
Dtx,Dty — similarly defined rotations of the feed radiation pattern with respect to its nominal direction (i.e. toward the secondary vertex).
Ray-tracing: N/R, Nu,Nv, u_max,v_max   3
N/R — integer number of equidistant points on the circular aperture radius. For unblocked aperture it cannot exceed 68 (corresponding to about 15000 points or rays to be traced). The larger is this number the slower is computing of power patterns.
Nu,Nv — integer number of nonzero spacial frequencies in u and v dimensions, respectively; |Nu| ≤ 200 and |Nv| ≤ 100. Total number of channels (points in the spacial frequency domain) is (|Nu| + 1)(|Nv| + 1). The larger is this number the slower is the execution. If Nu < 0 a reference power pattern is computed using the approximate (faster) DFT; if Nv < 0 the telescope power pattern is computed using the approximate DFT.
u_max,v_max — patterns are computed in the range –u_max to +u_max in u dimension and analogously in v.
Wavel,Taper,zAprt,pltD,pivot [m,dB,-]  4
Wavel — observing wavelength.
Taper — edge taper resulting from aperture illumination function. If Taper < 0 the Gaussian function is assumed; otherwise amplitudes will be weighted with the parabolic on pedestal function. Note that setting Taper = 0. reduces the parabolic distribution to the uniform one.
zAprt — z-coordinate of assumed aperture plane for ray tracing. On the program start it is computed as the mean of z of main dish center and that corresponding to X_min, but the user may change it as he wishes. With nonzero offsets the change will somewhat affect results due to divergence of the beam of rays, so it is better to keep it close to certain mean height reflection of rays on the paraboloid.
pltD — index that, if set to value less than 0., directs the program to compute and plot distortions of power pattern with respect to aberration free pattern (i.e. for just phaseless rays), otherwise the reference is made to phaseless amplitude computed for offset-free illumination function (i.e. symmetrical wrt to central ray) and for non-offset dish (i.e. with free-space taper computed for symmetrical antenna). If the absolute value of pltD equals to 1 power pattern is calculated for equivalent symmetrical telescope and no distortions looked for (effectively one obtains just the reference pattern used with pltD = 2.). The case pltD = 3. gives the same results as pltD = 2., except that full phaseless pattern is also computed (with pltD = 2. only the peak power is computed; see the partial output tables in .txt file).
pivot — defines rotation center of secondary mirror in the x-z plane. If integer part of pivot equals to 1, coordinates of the center (x_piv and z_piv) are calculated at the vertex of subreflector; if it equals to 0, they are zeroed (which means prime focus location of the center) and other values for pivot direct the program to accept the existing values. In the latter case (e.g. pivot = –1.) the next (6th) line in the .par file is also read and the first two values are taken for x- and z-coordinate of subreflector rotation center. Later they will not be modified by the program unless the user changes the pivot parameter value to one with 0 or 1 as its integer part.
       COMPUTED param's  & x,z_piv,du,v,N
The first two parameters in this row are x and z coordinates of the secondary center of rotation (in x-z plane), then there are working parameters of little use (they represent shift in u and v of peak of power pattern wrt direction of central ray — the one reflected from the secondary at the point corresponding to center of aperture in absence of offsets; the last parameter represents the same shift but in numbering of pattern points as organized in the internal array; its nonzero value may be indicative of excessive distortions of the main beam in u (N is a small integer) or v (|N| is close to a multiple of |Nu|).

Output parameters

Main dish: f/D,f/Dp, D2, G,coordG [m]  6
f/D — focal distance to diameter ratio.
f/Dp — same for parent (symmetrical, nonoffset) paraboloid.
D2 — the other (longer, 'vertical') physical diameter of main dish.
G — maximum depth of the dish.
coordG — distance from inner (closer to z-axis) edge of offset paraboloid to the point of maximum deph measured on the plane of dish rim.
 angles: t_V,t_H, t_0, t_1, t_2 [deg]  7
t_V — dish vertical opening (subtended) angle.
t_H — dish lateral ('horizontal') opening angle.
t_0 — center of aperture angle as seen from prime focus.
t_1 — similar angle to the nearest point on dish rim.
t_2 — same to the most distant point on dish rim.
The last three angles are inclinations to –z axis. For symmetrical telescopes t_V = t_H, t_0 = 0 and t_1 = –t_2.
Subreflector: Dc,Dm,Ds2,Gs,Xs_min [m]  8
Dc — secondary mirror lateral diameter at center (below mid point of longest diameter).
Dm — lateral diameter corresponding to maximum lateral diameter of main dish.
Ds2 — longest (vertical) diameter of secondary.
Gs — depth of subreflector at its center (mid of Dc).
Xs_min — distance from z-axis of closest point on the rim (for symmetrical telescope Xs_min = –Xs_max).
 feed angles: t_VH,t_c,_0,_1,_2 [deg]  9
t_VH — secondary opening cone (subtended) angle.
t_c — inclination angle of opening cone axis to z-axis.
_0 — similar inclination of imagery ray reaching center of aperture.
_1 — inclination of a ray through the innermost (closest to z-axis) point on dish rim.
_2 — same for the outermost point on dish rim.
Cassegrain: f1, f2, fs, F,   M  [m,-] 10
f1 — first focal distance (parent hyperboloid vertex to prime focus).
f2 — second focal distance (the vertex to secondary focus).
fs — interfocal distance (f1 + f2).
F — effective focal length (Mf).
M — magnification of the telescope (f2/f1).
Hyperbola: a,c,b,e,  FSAmpRatio [m,-] 11
a — constant of canonical (in Cartesian coordinates) form of parent hyperbola equation.
c — c2 = a2 + b2. Also: 2c = f1 + f2 = fs.
b — another constant of canonical form of parent hyperbola equation.
e — eccentricity of the parent hyperbola (c/a).
FSAmpRatio — amplitude ratio (rim to maximum) due to free-space taper.
Beam: t_u,t_v,     HPBW_u,_v,_0 [deg] 12
t_u — main beam direction offset downward in vertical plane (u-dimention).
t_v — same laterally leftward (corresponding to v-dimension).
HPBW_u — half power beam width in vertical plane.
_v — same in perpendicular plane (lateral offset).
_0 — beamwidth calculated traditionally.
  Squint, phi_X,Y, apert_X,Y  [deg,m] 13
Squint — absolute value of main beam offset (square root of t_u2 + t_v2).
phi_X — feed pattern inclination to YZ plane (+ve toward increasing x).
Y — same toward y (wrt XZ plane).
apert_X — x-coordinate on aperture of feed pattern maximum.
Y — same in y-direction.
Loss: Aberr,Spill,A&S,IllDec,Totl [%] 14
Aberr — aberration (coma and astigmatism) loss of antenna efficiency.
Spill — loss of gain due to increased spillover (percentage of power contained in rays reflected from the secondary but missing the main dish off the outer rim or into the inner hole).
A&S — the above two together.
IllDec — certain measure of loss due to decentered illumination (this envelopes both the loss resulting from asymmetry of the free-space taper present naturally in telescopes with offset aperture, and from adjusting the direction of feed radiation pattern through angles Dtx and Dty).
Totl — 'total' loss (useful for optimization purposes only).
  pathR, ComaLow,-Hi,Aber0,Totl'[m,%] 15
pathR — range of ray path (focus to aperture) errors (max – min).
ComaLow — lowest value of differential pattern (distortions due to offsets).
–Hi — highest value of the same.
Aber0 — aberration loss calculated according to approximate formulae.
Totl' — another measure of total gain or aperture loss.


When the following line appears:
 Press ENTER,q,$,c,m,o,a,p #i, or #i value (e.g. 40    0.060):
the user may issue one of the indicated commands, where:
Enter — evaluate (on change of any parameter) and display derived parameters for current set of input data.
q — causes writing the .txt and .ps files; if prior to q any parameter has been modified, the program after re-evaluation of derived parameters and writing these two files returns to the discussed command line prompt, otherwise the execution is stopped. This allows for viewing current pattern distortions with an auxiliary program (e.g. GSView) while continue running OptiCass. A 'qv' directs the program to plot (into the .ps file) a 'voltage' antenna radiation pattern instead of the power pattern (the former represents just the square-rooted power). Immediate quitting results after '' or after the standard '<Ctrl>C'.
$ — save all current parameters into a temporary file (for future use on start of another session with OptiCass). If the $ sign is followed by nonzero integer its nonzero units cause the program to write out to the .res file (residua), the nonzero tenths — to the .pat (power pattern) file and if the integer divided by 100 is nonzero then amplitude distribution array is saved to the .dat file.
c — clear all the 10 offsets, i.e. set 0. to parameters 1i and 2i, where i assumes values from 0 to 4.
m [n] — secondary mirror angular shape as seen by feed and its lateral diameter [in n sections; n (an odd number) defaults to 5]. The subreflector diameters (Diam. column in the command output table) are calculated at places (diameter mid-point distances from the telescope axis, listed in r(c) column) corresponding to n equidistant diameters of the circular aperture. This command just calls the SubShape routine.
o — optimize (for maximum gain) offsets in feed's z-coordinate and its pattern direction (Dtx and Dty parameters). See also the Remark in the –#i #j [step] command description below. Each o command is equivalent to typing first –22 140 (i.e. minimize the aberration), and then –23 144 (if the feed x is nonzero) and/or –24 144 (if the feed y is nonzero).
a — show classical approximations for feed offset gain losses and beam rates (see Fig. 4).
p #i — show pointing and loss (A&S) change due to parameter No #i. E.g. p 21 (21 stands for feed y-offset) might result in such response:
 Pointing errors per m or deg in u (elev) & v; dEff: 0.0' & -35.2'; 1.03%
which tells that feed lateral displacement of each 1 m results in about 35 arcmin main beam squint in opposite direction and about 1% loss in this antenna aperture efficiency. Caution: these numbers are obtained for 0.01 metre or degree change and then multiplied by a 100, thus they represent derivatives at current set of all the parameters rather than actual change. The assumed step (.01) may however happen to be inadequate for certain configurations.
#i [value] — display the current value of parameter No #i to 6 significant digits with short explanation of the parameter meaning. If the field 'value' is not empty the parameter No #i will be ascribed the value given, provided #i is less than 55 (i.e. belongs to the SETTABLE group); in this case all the settable parameters will be re-displayed. The user may set any value, not only those displayable in the format adopted, i.e. with 3 decimal digits.
n#i v1,v2,...,vn — set n new values, v1 to vn, to n consecutive parameters: No #i to #(i+n–1). The index (i+n–1) may assume values exceeding 4 to indicate parameters from next line (row following #).
–#i #j [step] — attempt to minimize parameter No #j with respect to #i [using new 'step' value in search for minimum (default is 0.1 m or deg)]; the newly set 'step' value will stay effective until next change; e.g. –21,144 .01 causes search for minimum of 'Total' loss by varying offset in y-coordinate of feed with step of 1 cm. Remark: The minimization is always done with following temporary presets for parameters Ray-tracing: N/R, Nu,Nv, u_max,v_max: 45. –14. –14. 1.2 1.2, i.e. using the faster DFT version and fixed parameters for pattern computations. Thus it is likely to fail to work for large offsets.

Fig. 4: An example of the 'a' command output. None of the results in this display is based on the ray-tracing technique — used are approximate formulae taken exclusively from literature with sole exception of the beam deviation for arbitrary 'z' of center of subreflector rotation (which is a combination of other formulae).


Fig. 5: Example of PostScript output for the TestO.par parameters. The black horizontal lines represent the ray tracing HPBW (with numerical value given in the upper figure). The HPBW0 calculated accordinag to a classical formula is delimited by upward arrowheads. The rays distributed on aperture plane (the right upper inset) are colored according to their amplitude (center is somewhat brighter) and phase (green corresponds to 0 and yellow to 360° of phase relative to that of central ray). Note a few rays marked by red dots that fell outside the dish (increasing the spillover).

The three files (.dat, .pat, .res) generated with e.g. $ 111 command are written in the form of 3- or 4-column tables, where the first two columns contain coordinates (x and y in .dat file, in metres, and u and v in the other two files). The first 3 rows of the .dat file may look like this:
  -15.6094   -3.1981  0.25704  0.38305
  -15.6095   -2.4173  0.26946  0.13258
  -15.6095   -1.6366  0.27833  6.16520
where the third column contains the ray amplitude and the last column — its phase in radians. This example was obtained for the 32OCRA parameters of Fig. 3.

The .pat and .res files both have the same structure. The third column contains just the pattern or differential pattern value at indicated spatial frequencies. Here are the first three rows of the 32OCRA.pat file, again generated with parameters of Fig. 3.

   -2.5000   -4.1403       350.4
   -2.2917   -4.1403       111.9
   -2.0833   -4.1403         3.1
The 32OCRA.res looks the same except that the third column contains in this case values 4.2, 4.4 and 0.9.

The .txt file consists of a list of various parameters followed by excerpts of the .dat, .pat and .res files. Ray amplitudes and phases are given here only for points originally lying (i.e. corresponding to offset-free positions) along x-axis and perpendicularly to it through the aperture centre. Similarly, the power patterns are represented by their sections along u and v through the mid-point of the array (corresponding to beam direction). The three rightmost columns represent the actual pattern, aberration-free pattern (computed from the same amplitudes but with phases zeroed), and reference pattern (for amplitudes weighted with illumination functions calculated as if there weren't any offsets), respectively.

The .ps file is a plot of the power pattern or 'voltage' pattern (i.e. the square root of power), residua and ray distribution on the aperture plane (see Fig. 5 for an example). It has also an outline of the two mirrors and some numerical results.


All the routines of OptiCas are original, written specifically for this purpose and usually with an eye on possibly shortest execution time. The most time-consuming are Fourier transforms. Because of the flexibility here implemented are direct algorithms (DFT as opposed to FFT). However, when largest arrays are being transformed the computing becomes somewhat too slow — with full precision it may take several seconds to once evaluate all the parameters. Therefore, for commands that involve many calls of DFT subroutine (such as search for minimum loss) OptiCass uses somewhat simplified DFT (simplification consists in using an array, presently every about 0.5°, of the sinus function values instead of a call to the function itself), which is about three times as fast but slightly less accurate. The user may choose to work with the faster versions by setting negative values to the Nfu and/or Nfv parameters.

Functional diagram
The figure that follows (Fig. 6) presents an OptiCass internal structure. Most of interactions with the user take place in the subroutine ReadParam called by the the main program (the latter reads initial .par files and writes .txt and .ps files). Essentially all the parameters are computed through call to the RTPattern subprogram. RTPattern calls the Amplitudes subprogram to get an amplitude distribution on the telescope aperture and transforms it to power pattern via DFT. Now the location of pattern peak is determined by a 3-point fit of parabola (subroutine FindMax finds an exact solution to this problem; independently for each of spatial frequencies) around a maximum value of the transform and then the pattern is computed anew for shifted (usually slightly) spatial frequencies so that now the peak exactly concides with the center of user defined u and v range. So found location of maximum is represented by the parameters t_u and t_v; the third parameter describing the beam deviation, Squint, is equal just to (t_u2 + t_v2)1/2.

Fig. 6: Building blocks of OptiCass. (This flow graph is based on an output from the Fujitsu Visual Analyzer that comes with the Lahey Fortran 90 package.)

Cassegrain geometry
Basic geometrical parameters are calculated from well known formulae (see e.g. this set of parameters) and are implemented in various parts of the program. Therefore adopting other than Cassegrain configurations, e.g. of the Gregorian telescope, although certainly possible, would not be a straightforward task.

Ray tracing
The RayTrace subroutine takes as the input coordinates of a point on the aperture (distance from paraboloid axis Ri, and angle seen from the axis between direction to the point and x-axis) to determine coordinates of the point of incidence on the subreflector assuming a perfect collimation. Then it computes (this time accounting for all the offsets) full path of a ray emitted from actual position of feed, and reflected from this physical point on the secondary mirror, up to the aperture plane as fixed by the user. As geometrical optics requires it, the rays are assumed to propagate along straight lines and follow strictly the law of reflection. The subroutine returns coordinates of the ray on the aperture plane, pathlength (distance from feed to aperture), direction cosines at the feed and at the aperture, and coordinates of point of reflection on the primary mirror. Also determined is the distance of reflection point on paraboloid from aperture center, required later to recognize rays missing the dish or those blocked by the subreflector.

The procedure is fast because care was taken to avoid trigonometric functions, wherever practical. Computation of all the amplitudes turned out to be about two orders of magnitude faster than to DFT them.

Amplitude distribution
To determine the power pattern, amplitudes are computed at number of points on an aperture plane and later are Fourier transformed to the spatial frequencies domain. The amplitudes, in OptiCass generated in the Amplitude subroutine, are weighted by one of the following two illumination functions:

[1 – τ (2 ri/D)2]/[1 + (0.5 Ri/F)2], used if Taper ≥ 0, or
exp[α 2 ri/D)2]/[1 + (0.5 Ri/F)2], used if Taper < 0
where the numerator represents the feed voltage pattern (parabolic on pedestal function or Gaussian function) and the denominator that corresponding to the so called free-space taper. In this expression D is the dish diameter, F — effective focal length, τ = 1 – 10–Taper/20, α = ln(10) Taper/20, Taper being one of the input parameters (assumed illumination function taper at the dish edge in dB), ri is distance on the aperture counted from a point corresponding to the direction of feed radiation pattern maximum (thus affected by Dtx and Dty parameters), and finally Ri is the distance of the i-th ray initial position on the aperture from the paraboloid axis. The last statement implies that the free-space taper function used here is completely independent of any offsets (save the paraboloid offset of course, if at all applicable). However it will be noted that the effect is accounted for anyway, owing to ray tracing itself, because the energy distribution defined by the applied free-space taper gets automatically modified according to changes of initially uniform distribution of rays on the aperture plane. Another way to say the same is that if we traced rays uniformly distributed on sphere centered at the feed, we would need not worry at all about any extra weighting to account for the free-space taper.

Thus weighted amplitudes are supplemented with signal phases (computed through division of ray pathlength by the wavelength) to form complex quantities.

Parallelly to the above described main process of computing the complex amplitudes the same procedure (Amplitudes) creates also similar array for reference, in which amplitudes are differently weighted and are phaseless (i.e. real). The weighting function is here of the same form but computed for completely offset-free telescope.

The aggregate power, PL (sum of squared amplitudes), associated with rays that miss the dish is used to calculate the parameter named Spill. For offset paraboloid telescope (more precisely, if X_min > 0) it is simply equal to
      Spill = 100 PL/PT,
where PT is the sum of squared amplitudes of all the rays traced. In the case of symmetrical telescope a correction to this quantity is made for rays that due to skewed main beam happen to fall within a subreflector shadow (are blocked by the secondary). Since in this small central part of the dish there are too few rays for sensibly continuous correction, use is made of a mathematical expression for the additionally blocked area. It is based on the fact that a fraction of area sticking out on one side (a crescent) of two equal circles with centers σ (here expressed in the circle diameters) apart is equl to

2[arcsin(σ) + σ (1 – σ2)1/2]/π.
This fraction is suitably (according to ray density) rescaled, approximately weighted and then added to PL before evaluating Spill. This method was found very inappropriate for similar estimation of power loss due to the rays missing the dish outside its outer edge because of considerably deformed ray density there (it is the same effect which is responsible for the earlier discussed free-space taper).

Power pattern and losses
Computation of power patterns and characteristics of a radio telescope performance based on ray tracing is all done in the RTPattern procedure. As mentioned, such pattern is obtained with the help of a Direct Fourier Transformation of amplitudes to spatial frequencies domain. The sum of squares of the real and imagery components of this transform is taken for the power pattern measure.

The spatial frequencies, u and v, as originally defined are number of (sinusoidal) cycles per the aperture radius (D/2). Given the observing wavelength λ, they may be conveniently expressed in terms of angular distance from the telescope symmetry axis:

u = (D/λ)sinΘx,   and   v = (D/λ)sinΘy,
where Θx/y are angles counted from the direction of parent paraboloid axis in the vertical plane (x-index) and one perpendicular to it (y-index). Thus, e.g. for D/λ = 32/0.06 (observing frequency of about 5 GHz and 32 m telescope) u or v equal to 1 corresponds to arcsin(0.001875) radians or about 6.5' (worth to note is that this unit is usually of the order of, but some 20% smaller than, the half power beam width, HPBW).

Up to three patterns are computed: the one from the complex amplitude distribution, P(u,v), one from the same amplitudes but phases zeroed, P0(u,v), and a reference pattern obtained from the auxiliary real 'offset-free' amplitudes, Po(u,v). One of the most important parameters, the aberration loss of antenna efficiency is calculated as
      Aberr = 100 (1 – Pmax/P0max),
where the subscript max indicates the peak of respective pattern. It is understood that this parameter includes the astigmatic and comatic fractional losses, which in classical approach are calculated separately. Other parameters related to the gain loss are calculated as follows:
      A&S = 100 – (100 – Aberr) (1 – Spill/100)
      IllDec = 100 (1 – P0max/Po'max)
      TotL = 100 – (100 – Aberr) (1 – Spill/100) (1 – IllDec/100)
      TotL' = 100 [1 – (Pmax/P0max) PT/Pomax]
      ComaLow = 100 PDmin/Pmax
      –Hi = 100 PDmax/Pmax
where Po'max is same as Pomax but computed with zeroed amplitudes of lost rays (rays spilled over), PDmin/max refer to the extreme values of a differential pattern. The differential pattern is obtained from the 'offset' pattern, P, by subtraction from it of the reference power pattern, Po, normalized to Pmax, i.e. PD = P – Po (Pmax/Pomax). The .res file and the distortions plot in .ps file will usually contain just these differences. Should, however, the user set pltD = –2., the reference would be made instead to the pattern obtained from 'offset' amplitudes but with phases zeroed (P0) and the 'offset-free' amplitudes would be discarded altogether (Po not computed).

The power pattern P(u,v) is further analysed in the vicinity of its maximum for points of half maximum, separately in u and v dimensions on either side of the maximum. Each of the 4 points is determined by linear interpolation between nearest two discrete values of P found. These then are used to determine HPBW_u and _v parameters.

Closed formulae
To have an 'on line' check on results obtained from the ray tracing analysis, the OptiCass computes also some parameters based on closed formulae taken from literature. All of them are grouped in the subroutine Formulae, although a few are used also elsewhere (for the subroutine is executed only when the user issues the 'a' command and to write to the .txt file). The sources are as follows:

J.W.M. Baars, 2003: Characteristics of reflector antenna, ALMA Memo 456.
J.W. Lamb, 1999: Optimized Optical Layout for MMA 12-m Antenna, MMA Memo 246.
R. Padman, 1995: Optical fundamentals for array feeds, Multi-Feed Systems for Radio Telescopes, (D.T. Emerson and J.M. Payne, Eds), ASP Conference Series, vol. 7.
P.-S. Kildal, 1983: The effects of subreflector diffraction on the aperture efficiency of a conventional Cassegrain antenna—An analytical approach, IEEE Trans. Antennas Propagat., vol. AP-31, no. 6, pp. 903–909, Nov. 1983.
Useful information in this respect is also contained in:
B. Shillue, 1997: Gain Degradation in a Symmetrical Cassegrain Antenna Due to Laterally Offset Feeds, MMA Memo 175.
B.J.Butler, 2003, Requirements for Subreflector and Feed Positioning for ALMA Antennas, ALMA Memo 479.

The beam deviation factor (BDF), used in a few formulae, is calculated according to the prescription given by Baars, however reduced to the following simpler equation:
    K = fD{2(τ fD + 1)[1 – fD log(1/fD + 1)] – τ}/(1 – 2 τ/3),
where τ = 1 – 10–|Taper|/20, and fD = (4 M f/D)2 with M being the magnification of Cassegrain system (M is set to 1 while computing the prime focus BDF, Kp). This assumes parabolic on pedestal illumination function and is used in OptiCass also for the Gaussian illumination function (therefore the absolute value in expression for τ). The HPBW_0, first side-lobe level, all the other beam deviations and feed offsets listed in the discussed table are also calculated according to Baars except for the case of secondary mirror rotation about center at an arbitrary z-coordinate. These formulae were checked to nicely agree with numerical results from OptiCass' ray tracing. The mentioned case of arbitrary z is calculated from a combination of Baars' expressions for translation and rotation about vertex, namely:
    [(2c + z_piv) K/M – z_piv Kp]/f.
(the derivation can be found in a separate document).
Formulae applied for calculation of gain loss due to diffraction on secondary and due to feed offsets as well as that for the Petzval radius are all adapted directly from Lamb.

Tests and first results

Number of checks are being carried to see how well the described program performs and to uncover possible bugs. It was found e.g. that the optimized aberration losses in antenna gain as output by OptiCass are frequently smaller than those found with the classical approximate formulae based on widely used J. Ruze expansion (see e.g. Shillue). OptiCass finds the losses according to unapproximated definition as the ratio of the peaks in power spectra computed from amplitude distribution complete with original phases (not even corrected for the slope due to beam squint) and the one devoid of phases.

For example, the cited Shillue's Memo states that the 8-m MMA telescope at λ = 0.35 mm and with the 12-dB parabolic illumination function would have 0.75% gain loss for 152.4 mm laterally offset feed (this is essentially, within 0.01% of, what obtains with classical approach formulae adopted in OptiCass). With the OptiCass optimization of z-coordinate of so offset feed (optimal z-offset turns out to be equal to 42.5 mm) the aberration loss is found to be 0.61% for the uniform illumination of the aperture due to the feed radiation pattern, and only 0.48% for parabolic 12 dB taper. This is about 2/3 of what Shillue has found — somewhat disturbing observation. The difference cannot be explained out by the free-space tapering (included in OptiCass), because in this example it amounts to only about 0.01 dB at the dish edge. With OptiCass it is possible to optimize also the feed radiation pattern direction. In this case, however, such optimization reduces the aberration losses only insignificantly.
We note on passing that with these feed offsets (lateral and axial) the range of pathlength errors (discussed in the following paragraph) is 24 mm i.e. ±34 λ, which corresponds to the beam squint of about 69 u-units (t_u = –0.173° or some 60 beamwidths).

In order to check that the OptiCass ray tracing algorithm does not contain any serious flaw the program output phases were compared with five examples based on ALMA telescope studies (which included ray tracing) performed by J.W. Lamb (in Butler, 2003) and published in the form of drawings. OptiCass produced results that were in very good agreement in all the 5 Lamb's cases as to the range of the path errors, shape of curve and the sign convention (except for the last case, where OptiCass uses opposite signing of the tilt angle).

1.925623 E–5 m — axial movement of feed of 1 cm
1.218378 E–3 m — lateral feed shift of 1 cm
5.565040 E–3 m — axial movement of secondary mirror of 1 cm
1.656128 E–2 m — lateral movement of secondary mirror of 1 cm
2.163380 E–2 m — tilt of secondary mirror of +1 or –1 deg

In the above table the first column is an eye estimate of what was drawn in the Lamb's Memo and the second column gives a full precision output of OptiCass parameter pathR (obtained by typing '150' on the command prompt). These many digits here are meant for other prospective users of the program as yet another check of their installation.

Fig. 7 presents more complete comparison for the 5-th of the discussed cases corresponding to most involved term (subreflector tilt about the mirror vertex) and quite an extreme situation (the tilt of 1° at 300 GHz would result in nearly 70% loss of gain!). Nevertheless, as seen in this figure, OptiCass numerical phases agree very nicely (although not perfectly) with Ruze's approximation formula. For this comparison the Ruze term was calculated from the expression given in Lamb's Tab. 1. The gap in the figure centre corresponds to the ALMA subreflector shadow on the aperture. This sort of data is normally contained in the OptiCass .txt file (amplitudes and phases along sections through the centre of aperture), even though the phase there is given in radians and in presence of large offsets may exhibit discontinuities of 2π (to avoid these jumps it suffices for the user to set a suitably longer wavelength, λ = Wavel). The phases multiplied by λ/(2π) correspond to the pathlengths modulo the wavelength.

Fig. 7: Comparison of path length errors in aperture distribution resulting from subreflector tilt of 1° computed for the ALMA 12 m telescope using OptiCass and one of J. Ruze's approximate formulas.

Certain practical analyses were carried with OptiCass in relation to the existing Torun 32-m radio telescope. These also demonstrate good agreement of ray tracing beam parameters (deviations and width) with classical calculations. An English reader may wish to verify this statement by having a glance at Tab. 1a therein. Fig. 1 in the same report gives an example of another useful application of OptiCass.

The described tests seem to demonstrate that OptiCass does not contain any fundamental fault in implementation of the ray tracing method, which would lead to underestimates of aberratic losses of antenna gain. If the program erratically estimated the position of DFT peak value, it would result in rather overestimated losses. Thus it is difficult to avoid the conclusion that the approximate solution of Ruze with an expression for the Seidel aberrations relying on smallness of phase errors may not be as straightforwardly applicable as presently assumed. An inadequate elimination of phase drifts across the aperture (that correspond to beam squint) from phase errors represented by the aberration expression (see Shillue's Memo) would lead to overestimates of losses since the losses would then represent just the antenna gain somewhere off the main lobe peak. Butler mentions of fitting a plane to the errors derived from approximate expressions while computing his losses. He doesn't give full details, but his estimates for ALMA feed offsets seem to be much overestimated. E.g. from his Fig. 2 we see about 1% gain loss for 10 λ of the lateral offset while both the Lamb's formulae and OptiCass indicate the loss to be a few orders of magnitude smaller (assuming λ = 0.46 mm). Similar situation is with his Fig. 1 while the other three estimates for subreflector are in good agreement with what OptiCass gave. Obviously enough, subtracting a (somehow) fitted plane does not guarantee that the rest of analysis is carried at exactly the beam peak, where the least losses occur. OptiCass does it by explicitly locating the peak position on the u-v plane.

Finally, it should be mentioned that an early (March 2005) version of OptiCass has been used for a preliminary performance analysis of an 80-m Cassegrain telescope with offset paraboloid (tackling with such a design has called for a special software and this program is the outcome). A short report is available here.


The idea of this program belongs essentially to prof. Andrzej Kus. Detailed algorithms and their implementation in FORTRAN are in fact mine, but virtually every step of development was preceded by discussions with him, so I can firmly state that the success of it would not be possible without his continued support and suggestions.

Notes on new versions

The above description refers to the OptiCass version 1.1 (of June the 1st, 2005).

New in version 1.2 (September the 1st, 2005)
— Added experimentally a multiple feed (array of up to 31×31 feeds) mode in which complex amplitudes corresponding to all feeds are summed for each ray computed to get the aperture amplitude and phase. Feeds are represented by 5 parameters each: x,y,z offsets (relative to offsets defined for the single-feed mode, i.e. as displayed in the line Feed offsets: x,y,z ... in Fig. 2 and 3), amplitude and phase. These can be read from console or a file with a structure similar to that given in the example of 'feeds.dat' file. Anything in it before the line 'Array size Nx,Ny' is assumed to be a comment. The two numbers (Nx,Ny) tell OptiCass how many feeds are in each column and row to be read subsequently. This structure may or may not correspond to physical dimensions of the array; e.g. the user may choose to arrange data of Nx×Ny = N feeds in one column specifying 'Array size N,1'. The 5 parameters of each feed may be modified after reading the file (upon OptiCass query). To initiate this multi-feed mode the user sets any nonzero value to the displayed fractional part of the parameter No 30 (old N/R changed now to N/R.NxNy)), then (on the following evaluation) he will be given option to enter filename with the feed data.
— If parameter h2 > 1000, it is replaced (in RTPattern subroutine) by f + Xs_max[(X_min + D)/(2f) – 2f/(X_min + D)], which corresponds to a flat subreflector.
— Added 1st side lobe level (along +ve u-axis) to parameter set.
— N = 200 changed to 1223 as used in SinArray part (entry) of the dftxyAProx subroutine (in order to better approximate DFT).
— Added analytic (somewhat involved) expressions of BDF for Gaussian illumination function (besides parabolic and uniform functions).
— Added analytic computation of gain loss due to tapered illumination for Gaussian, uniform & parabolic function.
— Added plot of power pattern cut-off at 1% level (g command).
— Also a number of bug fixes and cosmetic changes were made. Let the user pay special attention to new meaning of q and g commands as explained in the 'help' display. To be sure, check parameter meanings by typing their number alone.

New in versions 1.3 to 1.5 (September the 13th, 2005 to 19 April 2011)
— Improved and corrected partial blockage by subreflector (for D/2 > r2s > r1).
— Plotting of beam cut at 1 to 99 % level (e.g g500/g501 power/voltage at 50%).
— Corrected bug in f/D as printed to the .txt file.
— Input files 1 to 6 plus any other (no need to precede them with *).
— Modifications to 'help' outputs. — If Xs_max < 0 it is understood as –f2.
— added calculation of axis-to-outer-dish-edge angles (DishAngle) and of vertical extent of blind area on subreflector (exSub_blind).
— Changed format for Wavel on display (to 4 decimal places).
— Changed labels and info in graphs.

New in version 1.6 (May the 24th, 2011)
— One of 'if' instructions in RTPattern (for ddu) modified.
— Inscriptions on graph modified.
— Coma bar has now a sign (it is drawn up or down).
— Weighting of outer edge rays added (power lost is now calculated as squared sum of amplitudes).
— Precision of numbers on graphs changed (for small angles only).

New in version 1.7 & 1.8 (July 18th, 2011)
— Corrected drawing of subreflector outline on graph.
— New procedure option: "+par1,par2,step1,par3,step3,Ntimes,outfilename" — it tries to minimize par1 with respect to par2 using step1 and Ntimes repeating this, each time advancing par3 by step3; the results are written to outfilename file.
— New procedure option to scan in frequency with optimization at 100λ feed offset: "*40 fstep Nsteps outfilename". At each of Nsteps, wavelength (Wavel) is multiplied by fstep, par 20 (feed x-offset) is calculated as 100*Wavel, optimization in z-offset performed and results written to outfilename.
— Ray phase on graphs is now absolute (not relative to the centre).
— New option to shift pattern graph by fraction of the window: '54 +/–fract' (it spoils some measurements, e.g. aberration loss; see comments on usage in this document).
— Waighting of edge rays changed; now they extend beyond the edge up to 1/2 of the distance between rays; spillover corrected for this extra loss (in perfect alignment rays are situated partly outside dish).
— Improved sidelobe detection and done now on both sides of the main beam.

— KMB (kb@astro.uni.torun.pl)
May 2, 2005. (last modified 2011.06.24)