Some Basics of Radio Telescope Optics
as implemented in the OptiCass program


Free-space Taper

Consider the energy distribution on a circular aperture of a paraboloidal reflector resulting from illumination by a prime focus feed with uniform radiation pattern. A ray emitted at the angle θ away from symmetry axis of the paraboloid reaches the reflector at a point, whose distance from the aperture centre is

r = 2f tan(θ/2),
where f is the focal length of the telescope considered. This is just the equation of parabola. The energy contained between cones defined by θ and θ + dθ is cast onto a ring on our aperture limited by the radii r and
r + dr = 2f tan[(θ + dθ)/2] = r + f/cos2(θ/2) dθ
whose area is equal to 2πr[(r + dr) – r], i.e
2π 2f tan(θ/2) f/cos2(θ/2) dθ = 4πf2 sin(θ/2)/cos3(θ/2) dθ.
This area corresponds to energy flux through the spherical strip of area
2π(ρ sinθ) ρdθ.
When divided by the sphere radius squared ρ2 it becomes a solid angle equal to 4π sin(θ/2) cos(θ/2) dθ, which is proportional to total energy emitted in this angle. Thus on the aperture ring we have the energy density proportional to this last expression divided by the ring area on the aperture, i.e.
(1/f2)cos4(θ/2).
It can be shown easily that the same ralation holds also for a Cassegrain system, with the understanding that in this case θ is the angle between a ray originating in the feed (i.e. in secondary focus) and the paraboloid axis and f is replaced by F, the effective focal length. Thus in any case, relative to axial ray (with θ = 0) the energy density falls as cos4(θ/2). For example, the 32-m Torun telescope has subreflector subtended angle of 2θ = 18.826°, which means that signals reflected near edge of the dish at the feeds are received weakened by the factor of cos4(4.7065°) = 0.98658 or about –0.06 dB due to this so called free-space taper.

It follows that to account for the free space taper in simulation of amplitude distribution on the aperture of Cassegrain telescope one effectively multiplies the amplitudes by

cos2(θ/2) = 1/[1+0.25(r/F)2]

in addition to the taper function of the feed itself. In practice frequently this free-space illumination function is assumed to be included in the feed illumination function (feed radiation pattern), however in general, and especially when dealing with assymetrical telescopes, the two functions play different roles.

Ray Tracing Subroutine

This procedure is used to trace ray paths all the way between a feed and aperture plane with reflections from a hiperboloidal subreflector and paraboloidal main mirror. It consists of the following steps:

(1) A ray taken from a regular array on the aperture at the distance rp from the z-axis is projected from corresponding point on the paraboloid onto subreflector along the direction towards primary focus, assuming perfect alignment (no offsets) of the secondary.
(2) For the coordinates on the secondary, where this ray intersects the hyperboloid, the normal to it is calculated and rotated according to given angular offsets (in two planes, yz and xz) of the secondary.
(3) Coordinates of the point are also similarly rotated according to the same angular offsets and, additionally, translated by linear offsets of the secondary in x, y and z.
(4) Pathlength, direction cosines cx, cy and cz, and angle with respect to rotated normal, from the feed at its position (in general away from the secondary focus) to the point on subreflector are calculated.
(5) A few equations are used to solve for pathlength and coordinates (xp,yp) of the point of intersection on the paraboloid of the same ray after reflection from the secondary.
(6) Using the direction cosines of normal to paraboloid, the final section of ray path (direction cosines, cax, cay and caz, and distance) from paraboloid to aperture plain is calculated.

The procedure returns initial (at feed, cax,cay,caz) and final (cx,cy,cz) direction cosines, the pathlength (pathl) and coordinates of the reflection point on paraboloid (xp,yp) and coordinates on the aperture plane (Xapert,Yapert). These coordinates on aperture are used to transform (DFT) amplitude distribution into the telescope power pattern.

	subroutine RayTrace(rp,fih ,xp,yp, cax,cay,caz, cx,cy,cz,
     *                    Rdish,Xapert,Yapert,pathl)

c rp,fih        - initial/reference coordinates of ray on aperture plane,
c                 radial distance and angle
c xp,yp         - final coordinates of reflection point on paraboloid
c cax,cay,caz   - direction cosines at feed
c cx,cy,cz      - direction cosines at aperture
c Rdish         - radial distance from dish centre to reflection point
c Xapert,Yapert - coordinates of ray on aperture plane
c pathl         - length of ray path from feed to aperture plane

	implicit real*8 (a-h,o-z)
	data deg2rad/0.0174532925199433d0/
	common /params/ D,r1,f, r2s,hf2, xsub,ysub,zsub ,subtiltX,
     * subtiltY,xoff,yoff,zoff,tiltfX,tiltfy,datNr,parNfu,parNFv, ! par 9-17
     * freqmx,freqmy,wavel,TaperdB,zApert,differ,pivot,xpiv,zpiv
     *,ddu,ddv,dNsh ! par 18
     *,f_D,fbyD,D2,Depth,xDepth,t_V,  t_H,t_0,t_1,t_2	! par 30 - 39
     *,Dhc,Dhm,Dh2,g,r1s,fiV,fic,phi0d,fi1,fi2		! par 40 - 49
     *,f1,f2,fs,Feff,amagn,  a,c,b,e,Aratio  ,beamu	!     50 - 60
     *,beamv,HPBWu,HPBWv,fHPBWd, Squint,dphiX,dphiY,dfX0,dfY0
     *, Aberrl,SpillL,AberSp,Gillum,TotL
     *, ComaL,ComaH,deta_difr,gainl,TotL0 


	tiltsubX=subtiltX*deg2rad ! needed also in Amplitude procedure
	sintX=dsin(tiltsubX)
	costX=dsqrt(1-sintX*sintX)
	tiltsubY=subtiltY*deg2rad
	sintY=dsin(tiltsubY)
	costY=dsqrt(1-sintY*sintY)

c	Theta=2*datan(rp/(2*f))	   ! direction to focus (inclination to -z)
c							tanTh_2=rp/(2*f)
	sinTheta = rp/(f+rp*rp/(4*f)) 	 != 2*tanTh_2/(1+tanTh_2**2)
	cosTheta=(4*f*f-rp*rp)/(4*f*f+rp*rp) !=(1-tanTh_2**2)/(1+tanTh_2**2)
	b2 = (c+a)*(c-a)
	b2_denom = b2/(a+c*cosTheta)
c corresponding coordinates of this ray on hyperboloid
	z = -b2_denom*cosTheta
	r =  b2_denom*sinTheta
	x = r*dcos(fih)
	y = r*dsin(fih)
c Components of normal to hyperboloid at the above point:
	aN = dsqrt((a**2*r)**2+(b2*(z+c))**2)
	cfxr0 = x*a**2/aN
	cfyr0 = y*a**2/aN
	cfzr0 = -(z+c)*b2/aN
c rotate this unit vector with respect to the pivot in yz plane (by tY)
	cfxr = cfxr0
	cfyr = cfyr0*costY - cfzr0*sintY
	cfzr = cfzr0*costY + cfyr0*sintY
c and in xz plane (by tX)
	cfx = cfxr*costX - cfzr*sintX
	cfy = cfyr
	cfz = cfzr*costX + cfxr*sintX
c rotate coordinates of point as well (first in yz plane)
c in this plane we shall use 
	zpivY=0	! i.e. rotation pivot fixed in the prime focus
	rpivotY = dsqrt(y**2+(z-zpivY)**2)		! since ypiv=0
	if(rpivotY.ne.0d0) dzetaY = datan2(z-zpivY,y)
	yhp  =   rpivotY*dcos(dzetaY+tiltsubY)
	zyChng =  zpivY + rpivotY*dsin(dzetaY+tiltsubY) - z
c ...and xz plane
	rpivotX = dsqrt((x-xpiv)**2+(z-zpiv)**2)
	if(rpivotX.ne.0d0) dzetaX = datan2(z-zpiv,x-xpiv)
	xhp = xpiv + rpivotX*dcos(dzetaX+tiltsubX)
	zhpX = zpiv + rpivotX*dsin(dzetaX+tiltsubX)
	zhp=zhpX + zyChng

c now translate the coordinates from hyp. frame to parabloid frame
	xh = xhp + xsub
	yh = yhp + ysub
	zh = zhp + zsub
c pathlength from the secondary focus at (xoff,yoff,zoff - 2 c)
	pathl = dsqrt((xh-xoff)**2+(yh-yoff)**2+(zh-(zoff - 2*c))**2)
c direction cosines (from subreflector towards the feed)
	cx = (xh-xoff)/pathl
	cy = (yh-yoff)/pathl
	cz = (zh-(zoff - 2*c))/pathl
c double incident angle cosine C.N (i.e the vector product of C and N)
	cosIa2 = (cx*cfx+cy*cfy+cz*cfz)*2
c reflected ray cosines; in vector notation Cr = C - 2 Cf(C.Cf)
c C.Cf = cx*cfx+cy*cfy+cz*cfz =|C|*|Cf|*cos(C-Cf)=cosIa2/2
	chx = cx - cfx*cosIa2
	chy = cy - cfy*cosIa2
	chz = cz - cfz*cosIa2
c Find the point of ray intersection with paraboloid
c Equations to solve are: 	zp=(xp**2+yp**2)/(4f)-f,  xp=xh+dist*chx,
c    yp=yh+dist*chy,   	zp=zh+dist*chz, 	 dist=(zp-zh)/chz,
c or dist**2= [(yp-yh)**2+(xp-xh)**2]/(chx+chy)**2
c 0=[xh+(z-zh)*chx/chz]**2+[yh+(z-zh)*chy/chz]**2-4f(z+f), or
c 0=[xh*chz+(z-zh)*chx]**2+[yh*chz+(z-zh)*chy]**2-4f(z+f)*chz**2
	x1 = xh*chz - zh*chx
	y1 = yh*chz - zh*chy
c 0=[x1+z*chx]**2+[y1+z*chy]**2-4*f*(z+f)*chz**2, or
c 0=[x1+z*chx]**2+[y1+z*chy]**2-4*f*z*chz**2-4*f**2*chz**2
	aa = chx**2 + chy**2
	bb = 2*(x1*chx+y1*chy) - 4*f*chz**2
	cc=(x1**2+y1**2 - (2*f*chz)**2)
	delta = bb**2 - 4*aa*cc
		if(aa.gt.0.5) then
 	zp = (-bb - dsqrt(delta))/(aa + aa)
		if((zp - zh)*chz.lt.0d0) 
     *zp=(-bb + dsqrt(delta))/(aa+aa) ! the other solution
		else
	zp=2*cc/(-bb + dsqrt(delta))
		if((zp - zh)*chz.lt.0d0) 
     *zp=2*cc/(-bb - dsqrt(delta)) ! the other solution
		endif

	rpar   = 2*sqrt(f*(zp + f))

		if(abs(chz).gt..30d0) then ! to account for chz ~= 0
	dist = (zp - zh)/chz
		else 
c	dist=(4*f*(f+zh)-xh**2-yh**2)/(-2*f*chz + chx*xh + chy*yh+
c     + dSqrt((chx*xh+chy*yh-2*chz*f)**2 +aa*(4*f*(f+zh) -xh**2-yh**2))) 
c or equivalently:	rhopar**2 = zp**2+rpar**2
c 	dist+- = rhoCos +/- sqrt(rhoCos**2 + rhopar**2 - rho**2), thus
	rhoCos = xh*chx+yh*chy+zh*chz
	rhoSq = xh**2 + yh**2 + zh**2
	dist=dsqrt(rhoCos**2+zp**2+rpar**2-rhoSq) - rhoCos
		endif

	xp = xh + dist*chx	! reflection point on dish
	yp = yh + dist*chy
c Distance from aperture center, required later to recognize
c rays missing the dish or falling into the subreflector shadow
	Rdish = dsqrt((xp-r1-D/2)**2 + yp**2) 

c Normal to paraboloid
	aNr = 1/dsqrt(rpar*rpar + 4*f*f)
	cpx = xp*aNr
	cpy = yp*aNr
	cpz = -(f + f)*aNr
c incident ray angle cosine
	cosIp2 = (chx*cpx + chy*cpy + chz*cpz)*2
c reflected ray cosines
	cax = chx - cpx*cosIp2
	cay = chy - cpy*cosIp2
	caz = chz - cpz*cosIp2
	pathl  = pathl + dist
c intersection point on aperture plane that contains (0,0,zApert) point
c	dist=-(xp*cox+yp*coy+(zp-zApert)*coz)/(cax*cox+cay*coy+caz*coz)
	dist   = (zApert - zp)/caz
	pathl  = pathl + dist
	Xapert = xp + dist*cax - r1 - D/2	! <=== final
	Yapert = yp + dist*cay			! <===    coordinates
	end

Beam Deviation due to Subreflector Movements

Formulae for these were presented with enough detail by J.W.M. Baars. He has however neglected to give explicit expression for the beam offset in case of arbitrary location of the secondary rotation center. An expression apparently devised for this purpose is given in J.W. Lamb's report, but comparison with results obtained from ray-tracing suggests that as printed there it is not general enough. The Baars folmulae are based on the so called beam deviation factors, Kp and Ks, computed for offsets in prime focus and secondary focus, respectively. Any subreflector movement of interest can be expressed as sum of the rotation about its vertex through angle α and lateral translation by x. The corresponding small angle and displacement approximation formulae are given as:
    α(2c/f)(Ks + Kp)/(M + 1) = α(f1/f)(Ks + Kp), and
    (x/f)(Ks/M – Kp)
where 2c is the interfocal distance, f — the primary focal length, f1 — the secondary vertex to prime focus distance, and M — the Cassegrain system magnification. Assuming the center of subreflector rotation has coordinates (xc,zc) it is easy to show that a rotation of such subreflector by α (a small quantity) is associated with lateral shift of
    x = α (f1 + zc).
Note that this x is independent of xc. Now putting these three equations together one arrives at the desired formula, which reduces to this simple expression:

(α/f)[(2c + zc)Ks/M – zcKp].

It can be noted, that Lamb's formula obtains from the above when both the K quantities are removed from it (or equated to 1).