C     Last change:  KB    April 19, 2006
c19 Apr 2006: corrected bug in f/D printed to the .txt file
c14 Sep: plot beam cut at 1 to 99 % level (e.g g500/g501 power/voltage at 50%)
c13 Sep: partial blockage by subreflector improved/corrected (for D/2>r2s>r1)

c28 Aug:  added analytic expression of BDF for Gaussian illumination
c24 June: feed array enabled to be rectangular
c21 June: added gaussian/parabolic taper to the analytic losses
c19 June: gain loss due to tapered illum. added (Gaussian, uniform, parab.)
c17 June: 1% level pattern graph; SinArray N=200 changed to 1223
c16 June: Added 1st side lobe level (along +ve u-axis) to param's
c10 June: Now up to 31x31 array of feeds can be played with
c 8 June: Flat mirror amended; i.e when h2 = f+rs2((r1+D)/(2f)-2f/(r1+D))
c         if hf2.ge.1000 the above value is substituted (in RTPattern)
c 7 June: Multi-feed (5 feeds) option introduced
c 1 June: Added phase and amplitude to rays aperture distribution
c lf90 "<name>" -winconsole -bind -nwrap -s 32768 -pause -lib psplot plotsurf
	program OptiCass
c GEOMETRY AND OPTICAL PROPERTIES OF CASSEGRAIN RADIO TELESCOPE
c Analysis of geometry of Cassegrain type radio telescopes
c intended particularly for offset aperture telescopes

	implicit real*8 (a-h,o-z)
	real*4 xlen,ylen,fac,w,red,green,blue
	real*4 xd(15000),yd(15000),datAmp(15000),datPha(15000)
     *,xs(0:20300),ys(0:20300),Sre(0:20300),Sim(0:20300)
     *,Sre0(0:20300),Sim0(0:20300),ploty(0:300)
	character*80 fname,text,parRT(10)*8
	real*8 par(0:29)
	logical portr
	data pi/3.1415926535898d0/
	data parRT/'32OCRA','RT15','ALMA','p8MMA','p8015','TestO',
     *'p8010','p32-',2*' '/
	equivalence (par(0),D),(par(1),r1),(par(2),f),(par(3),r2s)
     *,(par(4),hf2),(par(5),xsub),(par(6),ysub),(par(7),zsub)
     *,(par(8),subtiltX),(par(9),subtiltY)
     *,(par(10),xoff),(par(11),yoff),(par(12),zoff)
     *,(par(13),tiltfx),(par(14),tiltfy),(par(15),datNr)
     *,(par(16),parNfu),(par(17),parNfv),(par(18),freqmx)
     *,(par(19),freqmy),(par(20),wavel),(par(21),TaperdB)
     *,(par(22),zAprt),(par(23),differ),(par(24),pivot)
c     *,(par(25),xpiv),(par(25),zpiv)
	common /arrays/ xd,yd,datAmp,datPha,xs,ys,Sre,Sim,Sre0,Sim0
	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
     *, pathRange,ComaL,ComaH,FstLob,TotL0
	common /output/ ie,in,ip,fname 

	parameter (MaxSide=31, Maxfeed=MaxSide*MaxSide)
	common /feeds/ feedx(Maxfeed),feedy(Maxfeed),feedz(Maxfeed),
     *feedAmp(Maxfeed),feedPh(Maxfeed)

	zhyp(theta)=a*(c+a*dcos(theta))/(a+c*dcos(theta))-c
c = -(c-a)*(c+a)*cos(theta)/(a+c*cos(theta))
	rhyp(theta)=(c-a)*(c+a)*sin(theta)/(a+c*cos(theta))
	rhohyp(theta)=(c-a)*(c+a)/(a+c*cos(theta))

c	read(*,*)

	print'(///T11,a/)','GEOMETRY AND OPTICAL PROPERTIES'
     *//' OF CASSEGRAIN RADIO TELESCOPE'
	print'(T21,a/)','OptiCass v. 1.2 (last revised Sep 14, 2005)'
	call HelpComm
	print*

	NrRT=10
	ip=0
	call getcl(text)
		if(text.ne.' ') then
	do 29 i=80, 1,-1
	if(text(i:i).eq.'\') go to 31
29	continue
c	print*,text(1:20)
	print*,' Dropped '//text(1:20)//
     *'... unrecognized. Ignoring it.'
	go to 4
31	do 32 j=i+1,76	!i+8
	if(text(j:j).eq.' '.or.text(j:j).eq.'.') go to 33
32	continue
	j=j-1
33	ie=j
	ip=i
	in=ie-ip-1
	iRT=in
	parRT(NrRT)=text(1+ip:ie-1)
	fname=text
	  if(text(i+1:i+12).eq.'working$.par'
     *.or.text(i+1:i+12).eq.'WORKING$.PAR') go to 35	
	go to 5
		endif

4	print'(8a)','  Enter 1 to 6 (',parRT(1)(1:7),(parRT(i)(1:5),
     &i=2,3),(parRT(i)(1:6),i=4,5),parRT(6)(1:5),'),'//
     &" *filename (.par extension    assumed!), $ (last saved setup)"
     &//', or ENTER for defaults: '
c     &//''
	read(*,'(a)') text
	if(text.eq.' ') go to 5
		if(text(1:1).ne.'$') go to 36
	fname='working$.par'
	inquire(file=fname,exist=portr)
			if(.not.portr) then
	Print*,'  Cannot find the saved parameter file! '
	go to 4
			endif
35	open(19,fname)
	read(19,*) (par(i),i=0,29)
	close(19)
	iRT=5
	in=iRT
	ie=ip+in+1
	iD=nint(D)
	if(iD.ge.100) write(fname(1+ip:80),'(2HRT,i3.3)') mod(iD,1000)
	if(iD.lt.100) write(fname(1+ip:80),'(a,i2.2,1h_)') 'RT',iD
c	write(fname(4:12),'(i2.2)') mod(idint(D),100)
	go to 10
c		endif	! end of $ option

36		if(text(1:1).eq.'*') then
	NrRT=10
	read(text(2:9),'(a)') parRT(NrRT)
		else	!go to 5
	read(text,*) NrRT
	NrRT=max(1,min(10,NrRT))
		endif

5	wavel=0.06d0!	wavelength in m

	datNR=20 !68 ! number of points on the aperture radius (max 68)
c	MaxNs=20300	! maximum number of points/bins in u-v pattern
c What to compute:  0 compute and draw only offset paraboloid pattern
	differ=2	! 1 compute (non-offset) normal Cassegrain pattern
			! 2 compute also difference offset - normal
			! 3 compute difference offset - normal and phaseless
			! differ < 0 - reference are just phaseless amplitudes

c============== Feed and subreflector mounting parameters========
c Offsets of feed relative to nominal secondary focus
	xoff= 0	! positive toward main dish
	yoff=0	! towards East when RT looks South
	zoff= 0 	! zoff=-hf2 corresponds to the parab. vertex level
	tiltfX= 0 	! extra feed tilt toward dish [deg]
	tiltfY= 0 	! lateral extra feed tilt [deg]

c tilt angle of subreflector (if pivot=1 center of rotation defined later
c as a point on subreflector surface where the central ray should fall)
	subtiltX=0 ! from x towards z direction, in degrees (up/down)
	subtiltY=0 ! from y (in yz plane; left/right)
	pivot=1 ! pivot location 1 - subrefl. center, 0 - prime focus, -1 set
c Position of (unrotated) subreflector relative to prime focus
c to be added to subreflector surface coordinates after rotation
	xsub=0 ! for offset paraboloid positive toward main dish
	ysub=0 ! leftward +ve
	zsub=0 ! upward +ve
c  ================== Power pattern parameters =================
c max. u (v) = (0.5 period)/(least gap between data) = +/- 0.5/1
c = +/- 0.5, but some tests show that max u may be even a multiple 
c of this value!
	freqmx= 2.5 ! maximum spatial frequency [in cycles per D/2];
	freqmy= 2.5 ! maximum spatial frequency [in cycles per D/2];
c these frequencies scale to angles using freqx = (D/lambda)sin(anglex) 
	parNfu= 24 ! number of channels (up to 200, apart of 0 frequency) along u
	parNfv= 24 ! <=100  these two numbers should be even to center at maximum
	TaperdB= 12 ! taper at dish edge [in dB], 12dB->taperP=.75

c Initial parameters of telescope geometry. Length units are meters.
	D = 80	! aperture diameter (the smaller diameter of dish rim)
	f = 63	! focal distance
	r1= 10 !-D/2! dish rim (inner tip) separation from paraboloid axis
	r2s = r1 	! distance of outer edge of subreflector from the axis
	hf2 = 0 	! height of secondary focus above paraboloid vertex

	do 7 iRT=8,2,-1
	if(parRT(NrRT)(iRT:iRT).ne.' ') go to 8
7	continue
8		if(iRT.lt.2) then
        parRT(NrRT)='RT'
        iRT=2
                endif




		if(ip.eq.0) then
	fname=parRT(NrRT)(1:iRT)//'.par'
	in=iRT
	ie=in+1
		endif
	inquire(file=fname,exist=portr)
			if(.not.portr) then
	in=5	!iRT=5
	ie=ip+in+1
	if(D.ge.99.5d0) write(fname(ip+1:80),'(a,i3.3)')
     & 'RT',mod(nint(D),1000)
	if(D.lt.99.5d0) write(fname(ip+1:80),'(a,i2.2,1h_)')
     & 'RT',nint(D)
		if(text.ne.' ') then
	Print*,'  Cannot find '//parRT(NrRT)(1:in)//'.par par'//
     *'ameter file!  Defaults assumed. Name: '//fname(1:ie-1)
		else
	Print*,'  Defaults assumed. Name: '//fname(1:ie-1)
		endif
	go to 10
			endif
	open(1,fname)

87	read(1,'(a)',end=34) text
	if(text(8:15).ne.'SETTABLE') go to 87

	do 9 i=0,4
88	read(1,'(a)',end=34) text
	if(text(1:1).eq.'c'.or.text(1:1).eq.'*') go to 88	!comment lines
9	read(text(41:80),'(5f8.3)') (par(i*5+j),j=0,4)
	if(idint(par(24)).eq.0.or.idint(par(24)).eq.1) go to 11
	read(1,'(a)',end=11) text
	if(text.ne.' ') read(text(42:80),'(f7.3,4f8.3)')
     * (par(5*5+j),j=0,1)
	go to 11
34	print*,' Failed par reading, so defaults. Name remains.'

11	close(1)
c	print'(3x,5f8.3)', ((par(i*5+j),j=0,4),i=0,4)

10	write(fname(ie:80),'(a)') '.ps'



	call SinArray	! generate table of sines for dftxyAProx

c Aperture plane height, upward positive (0 -> prime focus level)
c zApert=(z1+z2)/2, z-coord. of aperture plane, z1 = r1**2/(4*f)-f
	zApert=(r1**2+(r1+D/2)**2)/(8*f)-f 	! z-coord. of aperture plane

	
1000	call ReadPar(iGraph,phifx,phify,du,dv,u1,u2,Nmax,Smx,Ndat)

	Smax=dabs(Smx)

	open(6,fname(1:ie)//'txt')
	

	FSTaperdB =-20*dlog10(Aratio)
	fGHz=299.792458d-3/wavel
	rhoh2 =r2s*(f/(r1+D)+.25*(r1+D)/f)	!rhoh2 =r2s*(f/r2+.25*r2/f)
c	NRdat=datNr
	Nfu=dabs(parNfu)
	Nfv=dabs(parNFv)
	Ns=(Nfu+1)*(Nfv+1) 	! total number of channels
	u_vsize=dabs(parNfu)+1+(dabs(parNfv)+1)/1000.

	Write(6,'(a/(1x,a,T55,f9.3,1x,a))') 
     *'       Cassegrain radio telescope parameters'
     *,'Circular aperture diameter, D',D,'m'
     *,'Dish size (''vertical'' diameter), D2',D2,'m'
     *,'Focal length of primary mirror, f',f,'m'
     *,'Focal ratio of primary mirror, f/D',f_D,' '
     *,'Focal ratio of parent paraboloid, f/D_parent',fbyD,' '
     *,'Telescope F-number (or F/#): F_eff/D',amagn*f/D,' '
     *,'Deepest dish point is off mid point by',D2/2-xDepth,'m'
     *,'Dish depth, g_max',depth,'m'
     *,'Dish subtended angle - vertical plane',t_V,'deg'
     *,'Dish subtended angle - laterally',t_H,'deg'
     *,'Dish separation from axis, x1',r1,'m'
     *,'Axis to outermost dish point, x2',r1+D,'m'
     *,'Height of secondary focus above paraboloid vertex, h',hf2,'m'
     *,'Axis – subreflector outer edge distance',r2s,'m'
     *,'Prime focus – subreflector outer edge distance',rhoh2,'m'
     *,'Prime focus to hyperboloid vertex, f1=c-a',c-a,'m'
     *,'Secondary focus to hyperboloid vertex, f2=c+a',c+a,'m'
     *,'Secondary interfocal length, f1+f2=2c',2*c,'m'
     *,'Eccentricity of the hyperboloid, e=c/a', c/a,' '
     *,'Ray path difference to the foci, (f-h)a/c',(f-hf2)*a/c,'m'
     *,'Cassegrain telescope magnification, M=f2/f1',amagn,' '
     *,'Effective focal length, F=M*f',Feff,'m'
     *,'Aperture center angle as seen from sec. focus',phi0d,'deg'
     *,'Subreflector physical center angle',fic,'deg'
c     *(atan(r1/(2*Feff))+atan((r1+D)/(2*Feff)))*180/pi
     *,'Subreflector subtended angle',
     *(datan((r1+D)/(2*Feff))-datan(r1/(2*Feff)))*360/pi,'deg'
     *,'Subreflector size in vertical plane',Dh2,'m'
     *,'Subreflector lateral size at center',Dhc,'m'
     *,'Subreflector depth at center',g,'m'
c Radius of Petzval surface  (MMA memo 246/Born-Wolf)     
c Appr. radius of curvature of the Petzval surface is D*f*2/(Dhc+Dh2) 

	Write(6,'((1x,a,T40,3f8.3,1x,a),/1x,a,T48,2f8.3,1x,a,
     *2(/1x,a,T40,3f8.3,1x,a)/(1x,a,T48,2f8.3,1x,a))') 
     * 'Subreflector pivot coordinates, x,y,z',xpiv,0.d0,zpiv,'m'
     *,'Subreflector tilt towards +ve X and Y',subtiltX,subtiltY,'deg'
     *,'Subreflector offsets in x, y and z',xsub,ysub,zsub,'m'
     *,'Feed offsets in x, y and z',xoff,yoff,zoff,'m'
     *,'Feed pattern nominal angles wrt YZ,XZ planes',phifX,phifY,'deg'
     *,' and respective offsets',tiltfX,tiltfY,'deg'
     *,'Observing frequncy [GHz] and wavelength',fGHz,wavel*100,'cm'

	write(6,'(1(a,T55,f9.3,a/),1(a,T55,i9/),(a,T55,f9.3,a))')
     * ' Assumed z-coordinate for aperture plane',zApert,' m',
     * ' Number of rays traced',Ndat,
     * ' Number of pattern bins in u.v dimensions',u_vsize,' ',!Nfu+1,
     * ' Span of u coordinate [spacial cycles]',2*freqmx,' ',
     * ' Span of v coordinate',2*freqmy,' ',
     * ' Beam offset, vertically (downward +ve) ',beamu,' deg',
     * ' Beam offset, laterally (leftward +ve)',beamv,' deg',
     * ' HPBW, beamwidth in vertical plane',HPBWu*60,' ''',
     * '   same in orthogonal plane (laterally)',HPBWv*60,' ''',
     * ' Free space illumination taper or just asymmetry',FSTaperdB,
     *' dB', ' Feed illumination taper at the edge',TaperdB,' dB',
     * ' Ray tracing aberration loss',Aberrl,' %',
     * ' Ray tracing spillover loss',SpillL,' %',
     * ' Loss due to decentered illumination',Gillum,' %'

	if(dabs(differ).ge.2) write(6,'(a,T46,2f9.3,a)')
     * ' Pattern distortions range (min and max):',ComaL,ComaH,' %'
     *, ' '


c	Nside=min(MaxSide,idint(.5d0+dmod(dabs(datNr),1d0)*1d2))
	Nsidex=min(MaxSide,idint(dmod(.1d0+dabs(datNr)*1d2,1d2)))
	Nsidey=min(MaxSide,idint(dmod(.1d0+dabs(datNr)*1d4,1d2)))
	if(Nsidey.eq.0) Nsidey=Nsidex
		if(Nsidex*Nsidey.gt.1) then
c	Nfeeds=Nside*Nside
	write(6,'(a/a,2i3/1x,3hFNo,i7,34i8)') '   Feed array offsets [m],'
     *//' amplitudes and phases [rad]',
     *' Array size',Nsidex,Nsidey,(i,i=1,Nsidey)
	do 60 i=1,Nsidex
	iNs=(i-1)*Nsidey
	write(6,'(1x,i3,1hx,35f8.3)') i,(feedx(j+iNs),j=1,Nsidey)
	write(6,'(4x,1hy,35f8.3)') (feedy(j+iNs),j=1,Nsidey)
	write(6,'(4x,1hz,35f8.3)') (feedz(j+iNs),j=1,Nsidey)
	write(6,'(4x,1ha,35f8.3)') (feedAmp(j+iNs),j=1,Nsidey)
60	write(6,'(4x,1hp,35f8.3)') (feedPh(j+iNs),j=1,Nsidey)
	write(6,*)
		endif

	call Formulae

c*****
c To extract rays on the x & y sections through center of aperture 
c repeat the process done in Amplitudes routine.
	dR=D/(2*idint(datNr)+1)	! array step for initial aperture amplitude distr.
	NRdat=datNr

	write(6, '(a/a,i5,a,i3,a,f10.7,a/a)') 
     *' Ray data - sample of amplitude distribution on aperture', 
     *' No of rays:',Ndat,'.  Rays per diameter:', 2*Nrdat+1, 
     *'.  Step:',dR,' m','    (x [m],    y [m],  Ampl,  Phase [rad])'

c First extract points originally lying exactly along x-axis:
	write(6,*) ' Section along x-axis'
	R0=r1+D/2
	Ndi=0	! ray counter (index of the arrays)
	do 100 i=-NRdat,NRdat
	Xi=i*dR 				! original x coordinates (in m)
	Nxi=dsqrt((D/2)**2-Xi**2)/dR	! 2*Nxi+1 is the number of y-points at Xi
	if(i.eq.0) Ndi0=Ndi	! remember this for othogonal section
	do 100 j=-Nxi,Nxi
	Yi=j*dR 				! original y coordinates (in m)
c these (Xi,Yi) are starting coordinates on circular aperture
	rp=dsqrt((R0+Xi)**2+Yi**2)	! aperture point distance from z axis
		if(rp.lt.r2s) then	! region blocked by subreflector
	if(j.eq.0) write(6,'(2f10.4,2f9.5)') Xi,Yi,0.,0.
		else
	Ndi=Ndi+1
	if(j.eq.0) 
     * write(6,'(2f10.4,2f9.5)')xd(Ndi),yd(Ndi),datAmp(Ndi),datPha(Ndi)
c####     & *wavel*1000/(2*pi)	! in mm
		endif
100	continue
c Now in the direction of y-axis
	write(6,*) ' Central section perpendicular to x-axis'
	Ndi=Ndi0
	Xi=0
	do 101 j=-NRdat,NRdat
	Yi=j*dR 				! original y coordinates (in m)
	rp=dsqrt((R0+Xi)**2+Yi**2)	! aperture point distance from z axis
		if(rp.lt.r2s) then	! skip region blocked by subreflector
	write(6,'(2f10.4,2f9.5)') Xi,Yi,0.,0.
		else
	Ndi=Ndi+1
	write(6,'(2f10.4,2f9.5)')xd(Ndi),yd(Ndi),datAmp(Ndi),datPha(Ndi)
		endif
101	continue
c******

	ideg=0
c	beamu=dasin(du*wavel/(2*pi))*180/pi
c	beamv=dasin(dv*wavel/(2*pi))*180/pi
		if(ideg.eq.1) then	! to zle dziala !!!
	do 44 i=0,Ns-1
	ys(i)=dasin(ys(i)*wavel/D)*180/pi
44	xs(i)=dasin(xs(i)*wavel/D)*180/pi
	du=dasin(du*wavel/D)*180/pi
	dv=dasin(dv*wavel/D)*180/pi
		endif

	write(6, '(/a/a,i5/(1x,2f10.4,3f12.1))') '  Power spectra', 
     *' Central section (u,v0,Power,PowerNoPhase,RefPattern); Nu ='
     *,Nfu+1,((xs(i+(Nfu+1)*j),ys(i+(Nfu+1)*j),Sre(i+(Nfu+1)*j)
     *,Sim(i+(Nfu+1)*j),Sre0(i+(Nfu+1)*j),i=0,Nfu),j=Nfv/2,Nfv/2)
	write(6, '(/a,i5/(1x,2f10.4,3f12.1))') 
     *' Central section (u0,v,Power,PowerNoPhase,RefPattern); Nv ='
     *,Nfv+1,((xs(i+(Nfu+1)*j),ys(i+(Nfu+1)*j),Sre(i+(Nfu+1)*j)
     *,Sim(i+(Nfu+1)*j),Sre0(i+(Nfu+1)*j),i=Nfu/2,Nfu/2),j=0,Nfv)


	inquire(file=fname,exist=portr)
	if(portr) call system('del '//fname)

	call newdev(fname,ie+2)
	portr=.true.
	call psinit(portr)
	call setcolr(1.,0.,0.)

	call plot(+.3,4.5,-3) !u gory strony
	xlen=7.0
	ylen=7.8*.75
	w=.1
	CALL RECTFILC( 0., 0., xlen, 0.,ylen,.9,1.,.9)	! RED, GREEN BLUE)
c	call SYMBol(.5,ylen+w/2,w*1.2,carg(fname),0.,11)
	fac=1. 
	call factor(fac)
	call setcolr(0.,0.,1.)
c	CALL SETLW (0.015)
c	call setfnt(19)


	xmin=xs(0)
	ymin=0
	xsc=xlen/amax1(1.,xs(Nfu)-xs(0))
	ysc=ylen*(1-(Nmax/(Nfu+1))/(0.+Nfv)*.5)/Smax
	half=.5d0
	iVolt=mod(iGraph,10)

c To plot absolute voltage pattern take sqrt of power 
	if(iVolt.eq.0) go to 201
	ysc=ylen*(1-(Nmax/(Nfu+1))/(0.+Nfv)*.5)/dsqrt(Smax)
	half=dsqrt(.5d0)

201	dely=0.5*ylen/Nfv/ysc
	CALL SETLW(0.01)

c	ityp=-2	!-/+2 = plot with/without shadowing, -1 simple shadowing
	ityp=-2+iGraph/1000
	i1perc=mod(iGraph,1000)/10

	call graph(xs(0),Sre(0),Sim(0),0,xmin,xsc,ymin,ysc,ylen,ityp)

	N0=(Nfu+1)*(Nfv/2)
	xs0=xs(N0+Nfu/2)
	ys0=ys(N0)

c40.00   70.     50.      2.500   2.5
	scPc=1
		if(i1perc.ne.0) scPc=.47

	do 20 j=0,Nfv
	do 200 i=0,Nfu
	ploty(i)=Sre(j*(Nfu+1)+i)
		if(i1perc.ne.0) then
	ploty(i) = 100./i1perc*ploty(i)		!	*.4!!
	rs=sqrt((xs(j*(Nfu+1)+i)-xs0)**2 + (ys(j*(Nfu+1)+i)-ys0)**2)
	if(rs.lt.1.4d0) ploty(i) = dmin1(Smax,0d0+ploty(i))		!Smax*.6
		endif
	if(iVolt.ne.0) ploty(i) = sqrt(ploty(i))
200	continue
	call graph(xs(j*(Nfu+1)),ploty(0),Sim(j*(Nfu+1)),
     *Nfu,xmin,xsc,ymin,ysc*sngl(scPc),ylen,ityp)!
20	ymin=ymin-dely/scPc

	call graph(xs(1),Sre(1),Sim(1),0,xmin,xsc,ymin,ysc,ylen,ityp)
	call setcolr(0.5,0.,0.)
	CALL SETLW (0.03)
	do 202 i=0,Nfu
	ploty(i)=sqrt(Sre(N0+i))
		if(i1perc.ne.0) then
	rs=sqrt((xs(N0+i)-xs0)**2 + (ys(N0+i)-ys0)**2)
	ploty(i) = sqrt(100./i1perc)*ploty(i)
	if(rs.lt.1.4d0) ploty(i) = dmin1(dsqrt(Smax),0d0+ploty(i))	!*.6
		endif
202	if(iVolt.eq.0) ploty(i)=ploty(i)**2
	call graph(xs(N0),ploty(0),Sim(N0),Nfu,xmin,xsc,
     *ymin+(Nfv/2+1)*dely/sngl(scPc),ysc*sngl(scPc),ylen,ityp)

	call setcolr(0.,0.,0.)
	CALL SETLW (0.04)
	call setfnt(30)
		if(iVolt.eq.0) then
	write(text,'(a,i4,a,f5.1,a)') "Power pattern of",idint(D),
     * " m Cassegrain antenna at",fGHz,' GHz'
	call kekSYMc(xlen/2,ylen+3*w,w*1.7,text,0.,54,1)
		else
	write(text,'(a,i4,a,f5.1,a)') "Voltage pattern of",idint(D),
     * " m Cassegrain antenna at",fGHz,' GHz'
	call kekSYMc(xlen/2,ylen+3*w,w*1.7,text,0.,56,1)
		endif
	call plot(xlen,0.0,3)
	call plot(-0.02,0.0,2)
	call plot(0.0,0.0,3)
	call plot(0.0,ylen,2)
	call plot(xlen/2,0.0,3)
	call plot(xlen/2,2*w,2)
	call plot(xlen,-0.02,3)
	call plot(xlen,2*w,2)
	call plot(0.0,ylen*.25,3)
	call plot(2*w,ylen*.25,2)
	call kekSYMc(2*w,-2.75*w,w*2,"u",0.,1,0)
	call kekNUM(xlen/2,-3.0*w,w*1.8,sngl(du),0.,3,1)
	call kekNUM(xlen,-3.0*w,w*1.8,sngl(freqmx+du),0.,3,2)
	call kekSYMc(-2*w,ylen*.5,w*1.8,"v",90.,1,1)
	call kekNUM (-2*w,ylen/4,w*1.8,sngl(dv),90.,3,1)
	call kekNUM (-2*w,-3*w,w*1.8,sngl(-freqmy+dv),90.,3,0)

c	fHPBWd=2*dasin(fHPBW*wavel/D)*180/pi  -> sin(fHP.../360*pi)*D/wavel
c	hpu=xlen/(2*freqmx)*fHPBW	! u value at half peak power
	hpu=xlen/(2*freqmx)*fHPBWd/360*pi *D/wavel
	ypos=.75*ylen*half+ylen/4
	if(i1perc.ne.0) ypos=0.94*ylen
c	call plot(xlen/2-sngl(hpu),sngl(ypos),3)
c	call plot(xlen/2+sngl(hpu),sngl(ypos),2)
	x1=(u1-xmin)*xsc
	x2=(u2-xmin)*xsc
	call plot(sngl(x1),sngl(ypos),3)
	call plot(sngl(x2),sngl(ypos),2)
	CALL AROHED(xlen/2-sngl(hpu),sngl(ypos)-.07,0.,.08,18.,0)	!270.
	CALL AROHED(xlen/2+sngl(hpu),sngl(ypos)-.07,0.,.08,18.,0)	! 90.
c	call setfnt(30)
	write(text,'(f6.2,1h'')') HPBWu*60
	call kekSYMc(sngl(x1+x2)/2.,sngl(ypos)+w,w*1.4,text,0.,7,1)
	call kekSYMc(xlen/2,sngl(ypos)-w*2.8,w*1.3,"HPBW0",0.,5,1)

203	if(dabs(differ).lt.2) go to 16! to skip plotting difference
	ysc=ylen*(1-(Nmax/(Nfu+1))/(0.+Nfv)*.5)/Smax
	dely=0.5*ylen/Nfv/ysc

	call plot(-.0,-4.3,-3)
	CALL SETLW (0.007)
	call setcolr(1.,0.,0.)

	SminD=ComaL/100*Smax
	SmaxD=ComaH/100*Smax
	scaleD1p=dmax1(SmaxD,dabs(SminD),1d-4)
c		if(scaleD1p.le.1d-8) go to 290
	scaleD=.3*Smax/scaleD1p


	write(6, '(/a,i5/(1x,3f11.2))') 
     *' Central section (u,v0,Difference); Nu =',Nfu+1,((xs(i+(Nfu+1)*j)
     *,ys(i+(Nfu+1)*j),Sim0(i+(Nfu+1)*j),i=0,Nfu),j=Nfv/2,Nfv/2)
	write(6, '(/a,i5/(1x,3f11.2))') 
     *' Central section (u0,v,Difference); Nv =',Nfv+1,((xs(i+(Nfu+1)*j)
     *,ys(i+(Nfu+1)*j),Sim0(i+(Nfu+1)*j),i=Nfu/2,Nfu/2),j=0,Nfv)

	do 17 i=0,Ns-1
17	Sim0(i)=(Sim0(i)-SminD)*scaleD


	xmin=xs(0)
	ymin=Sim0(0)

c	ityp=-2	!-2 = z cieniem, +2 = bez
	ityp=-2+iGraph/1000
	call graph(xs(1),Sim0(1),Sim(1),0,xmin,xsc, ymin,ysc,ylen,ityp)
     
	do 30 j=0,Nfv! Nfv/2 !
	call graph(xs(j*(Nfu+1)),Sim0(j*(Nfu+1)),Sim(j*(Nfu+1)),
     * Nfu,xmin,xsc,ymin,ysc,ylen,ityp)
30	ymin=ymin-dely
	call setcolr(0.,0.,0.5)
	CALL SETLW (0.03)
	call graph(xs(1),Sim0(1),Sim(1),0,xmin,xsc,ymin,ysc,ylen,ityp)
	call graph(xs(N0),Sim0(N0),Sim(N0),Nfu,xmin,xsc,
     *ymin+(Nfv/2+1)*dely,ysc,ylen,ityp)

	do 18 i=0,Ns-1
c	Sim0(i)=(Sim0(i)-SminD)*scaleD
18	Sim0(i)=Sim0(i)/scaleD+SminD

	call setcolr(0.,0.,0.)
	Foffset=xoff
	if(abs(yoff).gt.abs(xoff)) Foffset=yoff
	if(yoff*xoff.ne.0.) Foffset=dsqrt(yoff**2+xoff**2)
		if(dabs(Foffset).gt.0.0005d0) then
	write(text,'(a,f7.3,2h m)')
     *"Power pattern distortions with feed offset of",Foffset
	call kekSYMc(xlen/2,ylen/2+3.*w,w*1.7,text,0.,54,1)
		else
	call kekSYMc(xlen/2,ylen/2+4.5*w,w*1.7,
     *"Power pattern distortions due to offset(s)",0.,42,1)
		endif
	CALL SETLW (0.04)

	percbar=scaleD1p/Smax*100	!dmax1(,.01d0)

	if(percbar/100*ysc*scaleD*Smax.lt.ylen/12) percbar=percbar*2
	ybar=percbar/100*ysc*scaleD*Smax
	write(text,'(f7.3,2h %)') percbar
	if(percbar.lt..01d0) write(text,'(1p,g7.1e1,2h %)') percbar
	call plot(xlen/10-w,ylen/4+sngl(ybar),3)
	call plot(xlen/10+w,ylen/4+sngl(ybar),2)
	call plot(xlen/10+0,ylen/4+sngl(ybar),2)
	call plot(xlen/10+0,ylen/4,2)
	call plot(xlen/10+w,ylen/4,2)
	call plot(xlen/10-w,ylen/4,2)
c	call setfnt(30)
	call kekSYMc(xlen/10+w,ylen/4+sngl(ybar)/2-w,
     * w*1.6,text,0.,9,0)
	write(text,'(a,f6.3,2h %)') 'Spillover loss =', SpillL!PowerL/PowerT*100
	call kekSYMc(xlen-w,-3*w,w*1.6,text,0.,24,2)
	write(text,'(a,f6.3,2h %)') 'Aberration loss =', AberrL
	call kekSYMc(w,-3*w,w*1.6,text,0.,25,0)
	
c	hpu=xlen/(2*freqmx)*fHPBW	! position of u value at half peak power
	ypos=ylen/2-5.5*w
	hpu=(x2-x1)/2
	call plot(xlen/2-sngl(hpu),3.5*w+sngl(ypos),3)
	call plot(xlen/2-sngl(hpu),1.5*w+sngl(ypos),2)
	call plot(xlen/2-sngl(hpu),2.5*w+sngl(ypos),2)
	call plot(xlen/2+sngl(hpu),2.5*w+sngl(ypos),2)
	call plot(xlen/2+sngl(hpu),3.5*w+sngl(ypos),2)
	call plot(xlen/2+sngl(hpu),1.5*w+sngl(ypos),2)
c	call setfnt(30)
	call kekSYMc(xlen/2,3.5*w+sngl(ypos),w*1.4,"HPBW",0.,4,1)

16	if(i1perc.ne.0) go to 19
c Plotting distribution of points as an inset
290	fac=1.2*2/D!
	call plot(xlen-1.2,4.4+4.4,-3)

	call setcolr(0.8,.5,0.2)
	CALL SETLW (0.01)	!.006)
	CALL CIRCLE (0., 0., sngl(D/2*fac), .false.)
	dotsize=0.006*sqrt(15000./Ndat)


c ===========Trace the central ray=============
	call RayTrace(r1+D/2,0d0, xp,dfY0, cox,coy,coz,
     * cphifx,cphify,cphifz, Rdish,Xapert,Yapert,pathl)

c cox/y/z are direction cosines of central ray, dish to aperture plane
c cphifx/y/z are feed-subreflector_center direction cosines

	phase0=dmod(pathl,wavel)/wavel	!*2*pi
	do 301 i=1,Ndat
		Amp=amax1(0.,amin1(datAmp(i),1.))
	Pha=datPha(i)/(2*pi)-phase0
	Pha=dmod(1+Pha,1d0)
c green - phase = 0, yellow - phase = 2pi
	blue=.3+(1-Amp)*.7
	red=Pha
	green=.6+.4*Amp
	
	if(Amp.eq.0d0) then	! rays lost outside the dish
		red=1
		green=0	!.2
		blue=0.2
	endif
		call setcolr(red,green,blue)!*.9
301	CALL CIRCLE (Xd(i)*fac, Yd(i)*fac,sngl(dotsize), .true.)

c Mark the parent paraboloid vertex and x & y axes 
	call setcolr(0.6,0.,0.6)
c	CALL SETLW (0.01)
	sdiam=r2s-r1s	! subreflector diameter
	sradius=sdiam/2*fac	! or: Dhc/2*fac
c Draw subreflector projection
	call circle(sngl((r2s+r1s)/2-R0)*fac,0.,sngl(sradius),.false.)
	call setcolr(0.75,0.,0.0)
	xarr=r2s+dmax1(sdiam/2,D/8)-R0
	yarr=sdiam/2+dmax1(sdiam/2,D/8)
	CALL CIRCLE (-sngl(R0)*fac, 0., .025, .true.)	! focus
	CALL ARROW (-sngl(R0*fac),0.,sngl(xarr*fac),0.,
     * .13,20.,0)
	CALL ARROW (-sngl(R0)*fac,0.,-sngl(R0)*fac
     *,sngl(yarr*fac), .13,20.,0)
	call kekSYMc(sngl(xarr*fac)+w/2,-w/2,w*1.8,"x/u",0.,3,0)
	call kekSYMc(-sngl(R0)*fac,sngl(yarr*fac)+w,w*1.8,"y/v",0.,3,1)

19	call plotnd
	print*,'  See files  ===>  '//fname(1:ie+2)//
c     *' and '//fname(ip+1:ie)//'txt  (may be erased next time!!!)'
     *' and .txt  (may be erased next time!!!)'

	close(6)
	go to 1000
	end	




	subroutine graph(x,y,e,N,xm,xsc,ym,ysc,ylen,ityp)
c ityp < 0 error bars only
c ityp = 0 punkty
c      = 1 punkty z error bars
c       -1 broken line with simple shadowing
c        2 broken line
c       -2 broken line with shadowing and interpolation
c       >1 [octal] number of symbol to be plotted (155 - circle)
	implicit real*8 (a-h,o-z)
	real*4 xp,yp,xpr,ypr,radius,ylen, x,y,e
	parameter (Npt=300)
	dimension x(0:N),y(0:N),e(0:N),cien(-1:Npt),cien2(-1:Npt),
     * cien3(-1:Npt),nup(-1:Npt+1)
	data cien,cien2,cien3,nup/Npt*-1.d3,Npt*-1.d3,
     * Npt*-1.d3,6*-1d3,Npt*0,3*0/,radius/.02/
	if(N.gt.0) go to 100
	do 90 i=0,Npt
	cien(i)=-1000
	cien2(i)=-1000
	cien3(i)=-1000
90	nup(i)=0
	return
100	if(ityp.eq.-1) go to 200
	do 3 i=0,N
	xp=(x(i)-xm)*xsc
	yp=(y(i)-ym)*ysc
		if(ityp.eq.-3) then
	yp=(y(i)*100.-ym)*ysc
	yp=amin1(yp,ylen)
		endif
	if(i.eq.0) call plot(xp,yp,3)
	if(ityp.eq.-2.or.ityp.eq.-3) go to 12	! linia lamana z cieniem
	if(ityp.lt.0) go to 4
	if(ityp.eq.2) go to 2	! linia lamana
	if(ityp.le.1) call circle(xp,yp,radius,.true.)
	if(ityp.gt.2) call keksym(Xp,Yp,radius,ityp,0.0,-999,1)
	if(e(i)*ysc.le.radius.or.ityp.eq.0) go to 3
4	call plot(xp,amin1(yp+sngl(e(i)*ysc),ylen),3)
	yp=dmax1(0.d0,yp-e(i)*ysc)
2	call plot(xp,yp,2)
	go to 3

12	if(yp.gt.cien(i)) then
	 	if(nup(i-1).eq.0) then
	 call plot(xp,yp,2)	! plot
	 	else
c	cienP=dmax1(cien(i-1),cien2(i-1),cien3(i-1))
	den=(cien3(i)-cien(i-1)-yp+ypp)
	if(dabs(den).lt.1d-8) den=dsign(1d-1,den)
	ratio=(yp-cien3(i))/den
cc	print'(i5,10f7.4)',i,ratio,xpp,xpr,xp,cien(i-1),cien3(i)!,cienP
cc     *,ypp,yp
c*
	ratio=dmax1(-1d0,dmin1(0d0,ratio))
	xpr=xp+ratio*(xp-xpp)
	ypr=yp+ratio*(yp-ypp)
	 call plot(xpr,ypr,3)	! move (pen up)
	 call plot(xp,yp,2)	! plot (pen down)
		endif
	 nup(i)=0
c
	 if(nup(i+1).ne.1) 	cien2(i)=dmax1(cien(i),cien2(i),DBLE(yp))
c
	 if(nup(i-1).ne.1) 	cien3(i)=dmax1(cien(i),cien3(i),DBLE(yp))
c
	 if(nup(i-1).ne.1)	cien(i)=dmax1(cien(i),DBLE(yp))
	else
	 	if(nup(i-1).eq.1) then
	 call plot(xp,yp,3)
		else
	den=(cien(i)-cien2(i-1)-yp+ypp)
	if(dabs(den).lt.1d-8) den=dsign(1d-1,den)
	ratio=(ypp-cien2(i-1))/den
c*
	ratio=dmax1(0d0,dmin1(1d0,ratio))
	xpr=xpp+ratio*(xp-xpp)
	ypr=ypp+ratio*(yp-ypp)
	 call plot(xpr,ypr,2)	! plot
	 call plot(xp,yp,3)	! move (pen up)
		endif
	 nup(i)=1
	endif
	xpp=xp
	ypp=yp
3	continue
	return

200	do 30 i=0,N
	xp=(x(i)-xm)*xsc
	yp=(y(i)-ym)*ysc
	if(i.eq.0) call plot(xp,yp,3)
	if(ityp.eq.-1) go to 120	! linia lamana z cieniem
	if(ityp.lt.0) go to 40
	if(ityp.eq.2) go to 20	! linia lamana
	if(ityp.le.1) call circle(xp,yp,radius,.true.)
	if(ityp.gt.2) call keksym(Xp,Yp,radius,ityp,0.0,-999,1)
	if(e(i)*ysc.le.radius.or.ityp.eq.0) go to 30
40	call plot(xp,amin1(yp+sngl(e(i)*ysc),ylen),3)
	yp=dmax1(0d0,yp-e(i)*ysc)
20	call plot(xp,yp,2)
	go to 30
120	if(yp.gt.cien(i)) then
	 	if(nup(i-1).eq.0) then
	 call plot(xp,yp,2)
	 	else
	 call plot(xp,yp,3)
		endif
	 nup(i)=0
	 cien(i)=yp
	else
	 call plot(xp,yp,3)
	 nup(i)=1
	endif
30	continue

	end



	subroutine dftxyAP(x,y,Ampl,Phase,N,u,v,Are,Aim)
c   Simple two-dimensional discrete Fourier transform. Same as dftxy
c exept the input data are amplitude Ampl(x,y) and phase Phase(x,y)
c of the complex data to be transformed. Computed is only one point
c of complex spectrum (Are,Aim) in frequency domain at (u,v).
c The spatial frequencies are assumed to be scaled numbers of cycles
c (fu,fv) per data span (Dx,Dy): u = 2*pi*fu/Dx and v = 2*pi*fv/Dy.
c Output is :   Are + i*Aim = Sum{Ampl*exp[Phase - (u*x + v*y)]}

c Procedura oblicza jeden punkt zespolonego dwuwymiarowego widma 
c Fouriera (Are,Aim) funkcji zespolonej reprezentowanej przez tablice 
c x, y, Ampl i Phase, o wymiarach N. Czestosci kolowe u i v sa
c wyrazone w liczbie okresow (cykli, fu i fv) na przedziale danych
c (Dx i Dy) podzielonej przez te przedzialy i pomnozonej przez 2pi:
c (u=2*pi*fu/Dx i v=2*pi*fv/Dy):
c       Sum{Ampl*exp[Phase - (u*x + v*y)]} 	ewent. *(1/N)
c Procedure mozna uzyc do liczenia jednowymiarowej transformaty tak:
c	call dftxyAP(x,y,Ampl,Phase,N,u,0.d0,Are,Aim)
c gdzie tablica y(N) moze zawierac dowolne wartosci, np. same zera.
c Transformate odwrotna dostaniemy zmieniajac znak argumentow u i v:
c	call dftxyAP(x,y,AmplS,PhaseS,N,-u,-v,Are,Aim)

	implicit real*8 (a-h,o-z)
	real*4 x(*),y(*),Ampl(*),Phase(*)

	Are = 0
	Aim = 0
		if(N.gt.0) then
	do 1 i = 1,N
	faza = u*x(i) + v*y(i) - Phase(i)
	Are = Are + Ampl(i)*dcos(faza)
1	Aim = Aim - Ampl(i)*dsin(faza)
		else
	do 2 i = 1,-N
	faza = u*x(i) + v*y(i)
	Are = Are + Ampl(i)*dcos(faza)
2	Aim = Aim - Ampl(i)*dsin(faza)
		endif
	end


	subroutine dftxyAProx(x,y,Ampl,Phase,N,u,v,Are,Aim)!

cc This less accurate version is optimised for execution time;
cc it is about 3 times as fast as the original.
cc IMPORTANT: Must be initiated with 	call SinArray

c   Simple two-dimensional discrete Fourier transform. Same as dftxy
c exept the input data are amplitude Ampl(x,y) and phase Phase(x,y)
c of the complex data to be transformed. Computed is only one point
c of complex spectrum (Are,Aim) in frequency domain at (u,v).
c The spatial frequencies are assumed to be scaled numbers of cycles
c (fu,fv) per data span (Dx,Dy): u = 2*pi*fu/Dx and v = 2*pi*fv/Dy.
c Output is :   Are + i*Aim = Sum{Ampl*exp[Phase - (u*x + v*y)]}

c Procedura oblicza jeden punkt zespolonego dwuwymiarowego widma 
c Fouriera (Are,Aim) funkcji zespolonej reprezentowanej przez tablice 
c x, y, Ampl i Phase, o wymiarach N. Czestosci kolowe u i v sa
c wyrazone w liczbie okresow (cykli, fu i fv) na przedziale danych
c (Dx i Dy) podzielonej przez te przedzialy i pomnozonej przez 2pi:
c (u=2*pi*fu/Dx i v=2*pi*fv/Dy):
c       Sum{Ampl*exp[Phase - (u*x + v*y)]} 	ewent. *(1/N)
c Procedure mozna uzyc do liczenia jednowymiarowej transformaty tak:
c	call dftxyAP(x,y,Ampl,Phase,N,u,0.d0,Are,Aim)
c gdzie tablica y(N) moze zawierac dowolne wartosci, np. same zera.
c Transformate odwrotna dostaniemy zmieniajac znak argumentow u i v:
c	call dftxyAP(x,y,AmplS,PhaseS,N,-u,-v,Are,Aim)

	implicit real*8 (a-h,o-z)
	parameter (Ns=1223,Ns4=Ns*4,dNs4=Ns4,dNs45=dNs4+.5,
     *dj=Ns/1.5707963267949d0)
	real*4 x(*),y(*),Ampl(*),Phase(*),sinus(-01:9*Ns),sinu
	data pi_2/1.5707963267949d0/ 	!,sinu/0./
c	if(sinu.eq.0.) go to 100
99	Are = 0
	Aim = 0
		if(N.gt.0) then
	do 1 i = 1,N
	faza = u*x(i) + v*y(i) - Phase(i) 
c	j=mod(nint(faza*Ns/1.5707963267949),Ns4)	! slows
	j=dmod(faza*dj,dNs4) + dNs45
	Are = Are +Ampl(i)*sinus(j+Ns)
1	Aim = Aim -Ampl(i)*sinus(j)
		else
	do 2 i = 1,-N
	faza =u*x(i) + v*y(i)
	j=dmod(faza*dj,dNs4) + dNs45!
	Are = Are + Ampl(i)*sinus(j+Ns)	! cos component
2	Aim = Aim - Ampl(i)*sinus(j)		! sin component
		endif
	return

	Entry SinArray
c Prepare table of sines (once, i.e. before the first call only)
100	do 101 i=0,Ns
	sinu=dsin(i*pi_2/Ns)
	sinus(i)=sinu
	sinus(2*Ns-i)=sinu
	sinus(2*Ns+i)=-sinu
	sinus(2*Ns4+i)=sinu
101	sinus(Ns4-i)=-sinu
	do 102 i=1,Ns4
102	sinus(Ns4+i)=sinus(i)
cc	print'(5x,10f7.4)',(sinus(i),i=0,2*Ns4+Ns,Ns/10)
	end



	subroutine SubShape(NDiameters)
c Wish to determine how the subreflector is seen by the feeds
c (whether it is circular or not)
	implicit real*8 (a-h,o-z)
	parameter(NDiamMax=101)
	real*4 feedview(6,0:NDiamMax)
	data pi/3.1415926535898d0/
	equivalence (par(0,0),D),(par(1,0),r1),(par(2,0),f)
     *,(par(3,0),r2s),(par(4,0),hf2)
	common /params/ par(0:4,0:15)
	rhohyp(theta)=(c-a)*(c+a)/(a+c*cos(theta))

	a=par(0,11)
	c=par(1,11)
	amagn=par(4,10)
	N=max(min(NDiameters/2,NDiamMax/2),1)
c 2*N+1 is number of sampled horizontal diameters of aperture
	dD=D/(2*N) ! sampling step

	do 10 i=0,2*N
	rdc=r1+i*dD
	y=dsqrt(dmax1(0.d0,(D/2)**2-(rdc-r1-D/2)**2)) ! half dish diameter at rdc
	rS=dsign(dsqrt(rdc**2+y**2),rdc) ! signed distance to dish edge
	thetaS=2*datan(rS/(2*f))	 ! direction to focus 1
	rhoS = 2*f/(1+cos(thetaS))
	rhohS=rhohyp(thetaS)		 ! focus to subr. edge distance
	yh=y*rhohS/rhoS			 ! half of diameter on hyperboloid
c	tH=dasin(y/rhoS) !=dasin(yh/rhohS) - 1/2 of subtended angle of this diameter
	feedview(6,i+1)=2*yh

	rhS=rhohS*dsin(thetaS)		 ! subr. edge to axis distance
	rhSc=dsqrt(dmax1(rhS**2-yh**2,0d0)) ! same for center of diameter at this level
c	tplane=datan(y/rdc)	! side ray path plane to central plane angle
c	print*,rhSc-rhS*cos(tplane)	! = 0 exactly
	hc=2*c-rhohS*dcos(thetaS)	! height of diameter center above focus 2
	feedview(5,i+1)=rhSc
c	feedview(5,i+1)=hc
	dtfeed=datan(yh/dsqrt(rhSc**2+hc**2)) ! lateral angular extention is twice this
	fifeed=dsign(datan(rhSc/hc),rdc) ! angle axis to center of the diam.
	feedview(1,i+1)=fifeed*180/pi
	feedview(2,i+1)=dtfeed*180/pi *2 ! full width
	fiS=datan(rhS/hc)		! angle: axis to side of the diameter (subr. edge)

	fifeedp=2*datan(rS/(2*f*amagn)) ! angle: axis to side of the diameter
c	fifeedp1=datan(1/(2*c*(a+c*cos(thetaS))/
c     * ((c**2-a**2)*sin(thetaS))-1/tan(thetaS)))
c	print*,(fifeedp-fifeedp1)*180/pi! test OK

c	feedview(4,i+1)=(fifeedp-fifeed)*180000/pi
	feedview(4,i+1)=(fifeedp-fiS)*180/pi
10	continue

	Rf=(feedview(1,2*N+1)-feedview(1,1))/2
	cRf=dcos(Rf*pi/180)
	fic=(feedview(1,2*N+1)+feedview(1,1))/2

	do 11 j=1,2*N+1
c11      feedview(3,j)=feedview(2,j)-2*Rf/N*dsqrt(N*N-(j-N-1.d0)**2)
	fi_fic=(feedview(1,j)-fic)*pi/180
	Dia=2*dacos(dmin1(cRf/dcos(fi_fic),1d0))
11	feedview(3,j)=(Dia*180/pi-feedview(2,j))	!*1d6
	fi1p2=(fic-Rf)*pi/360
	fi2p2=(fic+Rf)*pi/360
	Aratio=(cos(fi2p2)/cos(fi1p2))**2


	if(N.lt.0) print'(a,9f10.3)','  Feed view (fic V H) [deg] Ar:'
     *, fic,2*Rf,feedview(2,N+1), Aratio
c     *,feedview(1,1),cos(fi1p2*2)
	fic=feedview(1,N+1)
	print*,'   Subreflector shape  (angles seen from feed in deg,'
     *//'     last 3 columns in m)',"   i    Phi  "//
     *"  Phi-Phi_m Full_width Off_circ. PhiS'-PhiS  Diam.  r(c) increm."
	feedview(5,0)=feedview(5,1)
	write(*,'(i5,3f10.5,2f11.7,3f7.3)') (j,feedview(1,j),-fic+
     *feedview(1,j),(feedview(i,j),i=2,4),feedview(6,j),feedview(5,j)
     *-feedview(5,N+1),feedview(5,j)-feedview(5,j-1), j=1,2*N+1)
	end


	subroutine RayTrace(rp,fih ,xp,yp, cax,cay,caz, cx,cy,cz,
     *                    Rdish,Xapert,Yapert,pathl)
	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
c     *,ddu,ddv,dNsh ! par 18
     *,cox,coy,coz ! 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
     *, pathRange,ComaL,ComaH,FstLob,TotL0 

	tiltsubX=subtiltX*deg2rad	! required also in the 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)
		if(aN.gt.1d-8) then
	cfxr0 = x*a**2/aN
	cfyr0 = y*a**2/aN
	cfzr0 = -(z+c)*b2/aN
		else		! case of flat subreflector
	cfxr0 = 0.
	cfyr0 = 0.
	cfzr0 = -1d0
		endif
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 rotation pivot fixed in the prime focus!
	zpivY=0	!zpiv
	rpivotY = dsqrt(y**2+(z-zpivY)**2)		! since ypiv=0
	if(rpivotY.ne.0d0) dzetaY = datan2(z-zpivY,y)
	yhp  =   rpivotY*dcos(dzetaY+tiltsubY)
cc	zhp0 =   rpivotY*dsin(dzetaY+tiltsubY)
	zyChng =  zpivY + rpivotY*dsin(dzetaY+tiltsubY) - z
c ...and xz plane
cc	rpivotX = dsqrt((x-xpiv)**2+zhp0**2)
	rpivotX = dsqrt((x-xpiv)**2+(z-zpiv)**2)
cc	if(rpivotX.ne.0d0) dzetaX = datan2(zhp0,x-xpiv)
	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	if(i*j.eq.1) print*,xh-x,yh-y,zh-z
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
c	print'(5f12.7)',xh,xoff,pathl,cx  ! 0. 0.2 9.1480485  -0.0218626
	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

c	print*,' Tutaj2!',f*(zp + f) <- happened to be about -2.D-14
	rpar   = 2*dsqrt(dmax1(0d0,f*(zp + f)))

		if(abs(chz).gt..30d0) then ! to account for chz ~= 0
	dist = (zp - zh)/chz
	if(dist.lt.0) print*,dist, '	- distance	!!!!!!'	
		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 distance to ap. plane perpendicular to central ray
c	dist=-(xp*cox+yp*coy+(zp-zApert)*coz)/(cax*cox+cay*coy+caz*coz)
c or intersection with aperture plane perp. to z-axis at (0,0,zApert)
	dist   = (zApert - zp)/caz
	pathl  = pathl + dist
	Xapert = xp + dist*cax - r1 - D/2	! <=== final
	Yapert = yp + dist*cay			! <===    coordinates
	end


	subroutine Amplitude(Power0,PowerT,Ndat)
	implicit real*8 (a-h,o-z)
	real*4 xd(15000),yd(15000),datAmp(15000),datPha(15000)
     *,xs(0:20300),ys(0:20300)
     *,Sre(0:20300),Sim(0:20300),Sre0(0:20300),Sim0(0:20300)
	data pi/3.1415926535898d0/

	common /arrays/ xd,yd,datAmp,datPha,xs,ys,Sre,Sim,Sre0,Sim0
	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
c     *,ddu,ddv,dNsh ! par 18
     *,cox,coy,coz ! 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
     *, pathRange,ComaL,ComaH,FstLob,TotL0 

	zhyp(theta)=a*(c+a*dcos(theta))/(a+c*dcos(theta))-c
c = -(c-a)*(c+a)*cos(theta)/(a+c*cos(theta))
	rhyp(theta)=(c-a)*(c+a)*sin(theta)/(a+c*cos(theta))
	rhohyp(theta)=(c-a)*(c+a)/(a+c*cos(theta))
	FAreaL(x)=(asin(x) + x*sqrt(1-x**2))*2/pi

	idiffer=idint(dabs(differ))
c Parabolic on pedestal and Gaussian illumination functions constatnts
	taperP = 1-10**(-dabs(TaperdB)/20.)  !used in: 1-taperP*(2*r/D)**2
	taperG = TaperdB*.115129255d0   !*ln10/20; in: exp(taperG*(2*r/D)**2)
c fHPBW - spatial frequency at half power of symmetrical beam
c HPBW of symmetrical (nonoffset) beam according to JWM Baars [2003,
c ALMA Memo 456] is tau*(0.12*tau-0.343)+1.243, where tau=1-taperP
c In our notation, half of this is [in rad]:
	fHPBW=(taperP*(0.12*taperP+0.103)+1.02)/2
	fHPBWd=2*dasin(fHPBW*wavel/D)*180/pi

	call Geometry

c	r2 = r1+D		! outermost edge of dish separation from paraboloid axis
	R0 = r1+D/2		! aperture center distance from main axis
	phi0=2*datan(R0/(2*Feff)) ! feed nominal angle for aperture centre
	phi0d=phi0*180/pi

	dR=D/(2*idint(datNr)+1)	! array step for initial aperture amplitude distr.
	dfX=dtan((phi0+tiltfX/180*pi)/2)*2*Feff-R0
	dfY=dtan(tiltfY/180*pi/2)*2*Feff
c if we want the feed pattern axis 
c parallel to optical axis we need to enlarge dfX and dfY (i.e. tiltfX/Y)
c by central ray inclination to -z-axis in xz, and yz plane

c	dfX= R0*(1-dsqrt(1-tiltfX/100d0))	! for Taper=1-.75*(Ri*2d0/D)**2
c	dphi = 2*datan((R0+dfX)/(2*Feff)) - phi0	
c	if(tiltfX.ne.0.) print'(a,4f10.3)',' Extra feed tilts [deg m]:'
c     *, tiltfX,tiltfY, dfX,dfY

	FST0=1/(1+(R0/(2*Feff))**2)	! free-space taper at aperture centre
	Power0=0
	PowerT=0
	PowerL=0
	Ll=0
	Ls=0
	pathMin=+1d6
	pathMax=-1d6

	call RayTrace(R0,0d0, xp,dfY0, cox,coy,coz,
     * cphifx,cphify,cphifz, Rdish0,Xapert0,Yapert0,pathl)

c Loops computing the amplitude distribution and phase on aperture
c plane according to traced rays from initially regular array

	Ndat=0	! number of sampled points counter
	NRdat=datNr
	Xblock=(r2s+r1s)/2
	Rblock=(r2s-r1s)/2

	do 100 i=-NRdat,NRdat 		! ray-tracing of all the rays
	Xi=i*dR 				! original x coordinates (in m)
	Nxi=dsqrt((D/2)**2-Xi**2)/dR	! 2*Nxi+1 is the number of y-points at Xi
c	Nxi=dsqrt((idint(datNr)+.5)**2-i**2)

	do 100 j=-Nxi,Nxi
	Yi=j*dR 				! original y coordinates (in m)
c these (Xi,Yi) are coordinates on physical circular aperture

	rp=dsqrt((R0+Xi)**2+Yi**2)	! aperture point distance from z axis
c	if(rp.lt.r2s) go to 100!	! skip region blocked by subreflector

		if(rp.eq.0d0) then
	fih=0
		else
	fih=datan2(Yi,R0+Xi) ! angle between central and current meridian
		endif

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

	rpblock=dsqrt((R0+Xapert-Xblock)**2+Yapert**2)
	if(rpblock.lt.Rblock) go to 100!	! skip feed region (below subreflector)
c	rpblock=dsqrt((R0+Xapert+Xapert0)**2+(Yapert+Yapert0)**2)! opposite to feed placement
c	if(rpblock.lt.Rblock) go to 100!	! skip region blocked by subreflector

	phase=dmod(pathl,wavel)/wavel*2*pi	! phase of this ray
	if(pathl.lt.pathMin) pathMin=pathl
	if(pathl.gt.pathMax) pathMax=pathl

c Compute amplitude distribution on the aparture
	Ndat=Ndat+1
	xd(Ndat)=Xapert
	yd(Ndat)=Yapert
c	rp=sqrt((R0+Xapert)**2+Yapert**2)  ! distance from z-axis
c do not use the above since parallel shift does not change waigting
c phihalf=datan(rp/(2*Feff))	! half phi
c Ampl=dcos(phihalf)**2		! following equivalent avoids trigonometry
	Ampl=FST0/(1+(rp/(2*Feff))**2)	! =cos^2(phi/2) free-space taper
	Rt2=(Xi-dfX)**2+(Yi-dfY)**2 ! at feed Ampl corresp. to, Xi not Xappert;
c this means that the feed pattern is directed along central ray (one that
c in symmetrical Cassegrain falls at centre of subreflector and without 
c offsets also to centre of aperture); see note at dfX & dfY (top)
	Ri2=Xi**2+Yi**2
		if(taperdB.ge.0.d0) then
	Taper=1-taperP*Rt2*(2d0/D)**2 ! feed illumination function (ampl.)
	Taper=dmax1(Taper,0d0)
	Taper2=1-taperP*Ri2*4d0/D**2
		else
	Taper=dexp(taperG*Rt2*(2d0/D)**2)
	Taper2=dexp(taperG*Ri2*(2d0/D)**2)
		endif

	PowerT=PowerT+(Ampl*Taper)**2	! summing total power
		if(Rdish.gt.D/2) then	
c Additional condition: .or.r2s**2.gt.xp**2+yp**2 does not work properly
c because of too few rays sampled near the subreflector; in fact it represents
c loss due to shadowing by the subreflector of sky rays coming on its opposite
c side. A separate correction will be made at the end of this routine.
	PowerL=PowerL+(Ampl*Taper)**2 ! power spilled over dish edge
	Ampl=0
	Ll=Ll+1
		endif
		if(r2s**2.gt.xp**2+yp**2) Ls=Ls+1	
	datAmp(Ndat)=Ampl*Taper		! ray amplitude
	datPha(Ndat)=phase		! 	& phase on aperture
c The same for symmetrical illumination (taper function)
	xs(Ndat)=Xi		! temporarily (for use in loop to label 200)
	ys(Ndat)=Yi
c The following Ampl2 gives reference to symmetrical paraboloid
	Ampl2=Taper2/(1+Ri2/(2*Feff)**2) ! =Taper2*cos^2(phi/2)=taper*FS_taper
c	if(Ampl.le.0d0) Power0L=Power0L+Ampl2**2
	Power0=Power0+Ampl2**2	! main lobe total power
	if(idiffer.eq.1) datAmp(Ndat) = Ampl2
	if(idiffer.ge.2) Sim0(Ndat)  = Ampl2
c	if(Ndat.lt.3) print*,' Taper*A,taperP*A',Taper*Ampl,Ampl2
c     *,tiltfX
100	continue
	SpillL=PowerL/PowerT*100
	pathRange=pathMax-PathMin

	if(r1.gt.-D/2+0.001) return	! correction only for symmetric telescopes

c ===========Trace the central ray=============
c*	call RayTrace(R0,0d0, xp,dfY0, cox,coy,coz,
c*     * cphifx,cphify,cphifz, Rdish0,Xapert,Yapert,pathl)
c cox/y/z are direction cosines of central ray, dish to aperture plane
c cphifx/y/z are feed-subreflector_center direction cosines

c	phifx=dasin(cphifx)*180/pi
c	phify=dasin(cphify)*180/pi
cc	dfX0=xp-R0
c	dfY0=yp

c Fraction of area sticking out on one side of two equal circles of diameter
c D and centers s apart: FArea=[asin(x) + x*sqrt(1-x**2)]*2/pi, where x=s/D
c
cc	x=dsqrt(dfX0**2+dfY0**2)/D
cc	FAL=FAreaL(x)
cc	FAL2=FAreaL(Rdish0/D)
cc	Ndat2=pi*D**2/(4*dR**2)
cc	print*,'  Ndat,Ll,Ndat1*FAL:',Ndat,Ll,Ndat2,Ndat2*FAL,Ndat*FAL2
cc     * ,Ls,datLs
c The above approximation for outer part of dish does not work because
c of assymmetry due to the free space taper.

c Fractional loss of rays due to shadowing of the subreflector
	datLs=FAreaL(dmin1(Rdish0/(2*r2s),1d0))*pi*(r2s/dR)**2
	rm=r2s+Rdish0/2				! middle of the area for weight
	Ampl=FST0/(1+(rm/(2*Feff))**2)	! =cos^2(phi/2) free-space taper
		if(taperdB.ge.0.d0) then
	Taper=1-taperP*(rm*2d0/D)**2 !
		else
	Taper=dexp(taperG*(rm*2d0/D)**2)
		endif
	PowerLs=datLs*(Ampl*Taper)**2	! effectively lost power
	SpillL=(PowerL+PowerLs)/PowerT*100

	end








	subroutine RTPattern(phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	implicit real*8 (a-h,o-z)
	real*4 xd(15000),yd(15000),datAmp(15000),datPha(15000)
     *,xs(0:20300),ys(0:20300)
     *,Sre(0:20300),Sim(0:20300),Sre0(0:20300),Sim0(0:20300)
	data pi/3.1415926535898d0/

	common /arrays/ xd,yd,datAmp,datPha,xs,ys,Sre,Sim,Sre0,Sim0
	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
     *, pathRange,ComaL,ComaH,FstLob,TotL0

	if(hf2.ge.1d3) hf2 = f+r2s*((r1+D)/(2*f)-2*f/(r1+D))	! flat subreflector!!!
	if(mod(idint(parNfu),2).ne.0) parNfu=parNfu+dsign(1d0,parNfu) ! they MUST
	if(mod(idint(parNfv),2).ne.0) parNfv=parNfv+dsign(1d0,parNfu) !   be even
 
	idiff=idint(differ)
c If differ < 0 reference for difference is aberration free pattern
c (i.e. for phaseless rays alone) otherwise reference is phaseless offset-free 
c taper (i.e. symmetrical wrt to central ray) and possibly for non-offset dish
c (i.e. with free-space taper computed for symmetrical antenna - see next call)
		if(dmod(datNr,1d0).eq.0d0) then
	call Amplitude(Power0,PowerT,Ndat)
		else
	call Amplitudes(Power0,PowerT,Ndat)
		endif

c ===========Trace the central ray=============
	call RayTrace(r1+D/2,0d0, xp,dfY0, cox,coy,coz,
     * cphifx,cphify,cphifz, Rdish,Xapert,Yapert,pathl)

c cox/y/z are direction cosines of central ray, dish to aperture plane
c cphifx/y/z are feed-subreflector_center direction cosines

	phifx=dasin(cphifx)*180/pi
c	print*,cphifx,phifx
	phify=dasin(cphify)*180/pi
	dfX0=xp-r1-D/2
c	dfY0=yp
cc	dphiX = 2*datan(xp/(2*Feff))*180/pi - phi0d
cc	dphiY = 2*datan((dfY0)/(2*Feff))*180/pi
	dphiX=phifx	
	dphiY=phify

c	print*,cox*D/wavel,coy*D/wavel
c===============================================

c DFT spectrum of data calculated on a rectangular grid
	Nfu=dabs(parNfu)
	Nfv=dabs(parNfv)
	Ns=(Nfu+1)*(Nfv+1) 	! total number of channels
	umax=2*pi*freqmx/(D/2)	! maximum u frequency
	vmax=2*pi*freqmy/(D/2)	! same for v

	if(idiff.lt.2) go to 201	!.or.iaber.eq.1
	do 200 i=0,Nfu/2		! u-axis points
	u=umax*(i-Nfu/2.d0)/max(Nfu,1) 	! i-th point frequency
	do 200 j=0,Nfv		! v-column frequencies for u-row
	v=vmax*(j-Nfv/2.d0)/max(Nfv,1)
	idat=i+(Nfu+1)*j
	idats=NS-1-idat
		if(idat.ge.Ns) print*,i,j,idat,Ns, ' (idat>=Ns,in RTPattern)'
c compute reference pattern for equivalent 'normal' Cassegrain
call dftxyAP(xd,yd,Sim0,Sre,-Ndat,u,v,Ar,Ai)
		if(parNfu.gt.0) then
	call dftxyAP(xs(1),ys(1),Sim0(1),Sre,-Ndat,u,v,Ar,Ai) !(xs,ys)<-(Xi,Yi)
		else
	call dftxyAProx(xs(1),ys(1),Sim0(1),Sre,-Ndat,u,v,Ar,Ai)
		endif
	Sre0(idat)=(Ar**2+Ai**2)*PowerT/Power0
	Sre0(idats)=Sre0(idat)		! DFT of reals is symmetrical
200	continue

c Shift the pattern by the direction of central ray (i.e ~main lobe)
201	idone=0
	du=+cox/wavel*2*pi	! *D
	dv=+coy/wavel*2*pi	! *D

c	do 210 idat=0,Ns-1
205	do 210 i=0,Nfu
	u=umax*(i-Nfu/2.d0)/max(Nfu,1)
	do 210 j=0,Nfv
	v=vmax*(j-Nfv/2.d0)/max(Nfv,1)
	idat=i+(Nfu+1)*j
	xs(idat)=u
	ys(idat)=v
c The pattern of prime interest
		if(parNfv.gt.0d0) then
	call dftxyAP(xd,yd,datAmp,datPha,Ndat,u+du,v+dv,Ar,Ai)
		else
	call dftxyAProx(xd,yd,datAmp,datPha,Ndat,u+du,v+dv,Ar,Ai)
		endif
	Sre(idat)=Ar**2+Ai**2
	if(idiff.eq.2.and.idat.ne.Ns/2) go to 210
c The same without accounting for phase due to ray pathlengths
		if(parNfu.gt.0d0) then
	call dftxyAP(xd,yd,datAmp,datPha,-Ndat,u,v,Ar,Ai)
		else
	call dftxyAProx(xd,yd,datAmp,datPha,-Ndat,u,v,Ar,Ai)
		endif
	Sim(idat)=Ar**2+Ai**2
210	continue


c	print*,' Sre0,Sre,Sim, i=10:',Sre0(10),Sre(10),Sim(10)
209	Smax=0
	Nmax=0
	do 15 i=0,Ns-1
	if(Sre(i).ge.Smax) then
	 Smax=Sre(i)
	 Nmax=i
	endif
15	continue


	  if(idone.ne.1) then
	dNsh=Nmax-Ns/2
	u0=xs(Nmax)	! for correction to coarse tune to maximum
	v0=ys(Nmax)
		if(mod(idint(dabs(dNsh)),Nfu+1).lt.Nfu/4) then
	SreM=Sre(Nmax)
			if(Nfu.gt.10) then	! arbitrarily chosen limit
c Parabolic fit to 3 points near pattern maximum:
	call FindMax(dble(xs(Nmax-1)),dble(Sre(Nmax-1)),dble(xs(Nmax))
     *, SreM,dble(xs(Nmax+1)),dble(Sre(Nmax+1)),u0)!,Sm
			endif

	ddu=u0-xs(Ns/2)	! correction to tune to maximum
	du=du+ddu
cc	print*,2,Nmax-(Nfu+1)

			if(Nfv.ge.Nfu/2) then	! arbitrarily chosen limit
c in v neighbouring values are at Nmax +/- (Nfu+1)
	ivm=max(0,Nmax-(Nfu+1))
	ivp=min(Ns-1,Nmax+(Nfu+1))
	Sre_1=Sre(ivm)
	SreP1=Sre(ivp)
	ys_1=ys(ivm)
	ysP1=ys(ivp)

	call FindMax(ys_1,Sre_1,ysP1,SreP1,dble(ys(Nmax)),SreM,v0)!,Sm
			else

c generate new orthogonal section through the beam (for u=du)
c	print*,3,Nmax-(Nfu+1)
	SreM=0
	do 220 j=0,Nfu		!same number of bins and step as in u
	v=umax*(j-Nfu/2.d0)/max(Nfu,1)
	idat=j+Ns
	ys(idat)=v
	call dftxyAP(xd,yd,datAmp,datPha,Ndat,du,v+dv,Ar,Ai)
	Sre(idat)=Ar**2+Ai**2
	if(Sre(idat).le.SreM) go to 220
	 SreM=Sre(idat)
	 Nmaxv=idat
220	continue
	call FindMax(dble(ys(Nmaxv-1)),dble(Sre(Nmaxv-1)),dble(ys(Nmaxv))
     *, SreM,dble(ys(Nmaxv+1)),dble(Sre(Nmaxv+1)),v0)!,Sm
			endif
		endif
	ddv=v0-ys(Ns/2)
	dv=dv+ddv
	idone=1
	go to 205
	  endif

c Determine width of beam in u
	half=Sre(Nmax)/2
		if(Nmax.lt.Ns-2) then
	do 211 i=Nmax+1,min(Ns-1,Nmax+Nfu/2-1),+1
	if(Sre(i).le.half) go to 212
211	continue
212	u2=xs(i)-(xs(i)-xs(i-1))*(half-Sre(i))/(Sre(i-1)-Sre(i))
		else
	u2=xs(Ns-1)
		endif

c find 1st sidelobe
		if(i+4.gt.Ns) then
	FstLob=Sre(Ns-1)
c	u3=xs(Ns-1)
		else
	do 2121 j=i,min(Ns-1,Nmax+Nfu/2-1),+1
	if(Sre(j).gt.half) go to 2122
2121	half=Sre(j)
2122	do 2123 i=j,min(Ns-1,Nmax+Nfu/2-1),+1
	if(Sre(i).lt.half) go to 2124
2123	half=Sre(i)
2124	FstLob=Sre(i-1)
c	u3=xs(i-1)
		 if(Sre(i).lt.Sre(i-1).and.Sre(i-2).le.Sre(i-1)) then
c Parabola through 3 points near side-lobe maximum:
	call FindMaxE(dble(Sre(i-2)),dble(Sre(i-1)),dble(Sre(i)),FstLob)
c	call FindDer(dble(xs(i-2)),dble(Sre(i-2)),dble(xs(i-1)),
c     *dble(Sre(i-1)),dble(xs(i)),dble(Sre(i)),u3,FstLob,der1,der2)
		 endif
		endif
	FstLob=FstLob/Smax*100	! 1st side-lobe in %

	half=Sre(Nmax)/2
		if(Nmax.gt.1) then
	do 213 i=Nmax-1,max(1,Nmax-Nfu/2+1),-1
	if(Sre(i).le.half) go to 214
213	continue
214	u1=xs(i)+(xs(i+1)-xs(i))*(half-Sre(i))/(Sre(i+1)-Sre(i))
		else
	u1=xs(0)
		endif

c Determine width of beam in v
		if(Nmaxv.le.Ns) then
	do 215 i=Nmax,Ns-1-(Nfu+1),Nfu+1
	if(Sre(i).le.half) go to 216
215	continue
216	i_1=i-(Nfu+1)
	v2=ys(i)-(ys(i)-ys(i_1))*(half-Sre(i))/(Sre(i_1)-Sre(i))
c	if(Nmaxv.eq.0) Nmaxv=Nmax
c	print*,i,Ns,Sre(i-1),Sre(i)!Nmax,Nmaxv,Nfv,Nfu

	do 217 i=Nmax,(Nfu+1),-(Nfu+1)
	if(Sre(i).le.half) go to 218
217	continue
218	i_1=i+(Nfu+1)
	v1=ys(i)+(ys(i_1)-ys(i))*(half-Sre(i))/(Sre(i_1)-Sre(i))
		else
	do 1211 i=Nmaxv+1,Nmaxv+Nfu/2-1,+1
	if(Sre(i).le.half) go to 1212
1211	continue
1212	v2=ys(i)-(ys(i)-ys(i-1))*(half-Sre(i))/(Sre(i-1)-Sre(i))

	do 1213 i=Nmaxv-1,Nmaxv-Nfu/2+1,-1
	if(Sre(i).le.half) go to 1214
1213	continue
1214	v1=ys(i)+(ys(i+1)-ys(i))*(half-Sre(i))/(Sre(i+1)-Sre(i))
		endif


c angle = asin(u*wavel/(2*pi))
	HPBWu=(dasin((u2+du)*wavel/(2*pi))-dasin((u1+du)*wavel/(2*pi)))
     **180/pi
	HPBWv=(dasin((v2+dv)*wavel/(2*pi))-dasin((v1+dv)*wavel/(2*pi)))
     **180/pi

c Main lobe direction
	beamu=dasin(du*wavel/(2*pi))*180/pi
	beamv=dasin(dv*wavel/(2*pi))*180/pi
	Squint=dsqrt(beamu**2+beamv**2)

c find spectral amplitude at frequency (0,0) of not decentered beam
	SmaxI=0
	do 116 i=1,Ndat
116	if(datAmp(i).gt.0.) SmaxI=SmaxI+Sim0(i)
	Gillum=(1-Sim(Ns/2)/SmaxI**2)*100	! loss due to taper shift

c	print*,'Smax,Sim(Ns/2),SmaxI',Smax,Sim(Ns/2),SmaxI**2,Gillum
c     *,Power0,PowerT
cSmax,Sim(Ns/2),SmaxI    450815.687500000         633462.

c Scale spacial frequencies, u and v, back to cycle numbers, so that 
c now they conveniently equal to (D/lambda)*sin(angle)
	uvsc=D/(2*pi)
	Smax0=0
c	Nmax0=0
	du=du*uvsc
	dv=dv*uvsc
	u1=u1*uvsc+du
	u2=u2*uvsc+du
c	u3=u3*uvsc+du
c	dNsh=u3

	do 18 i=0,Ns-1
	xs(i)=xs(i)*uvsc	+du	! added May 9
	ys(i)=ys(i)*uvsc	+dv
	if(Sre0(i).ge.Smax0) then
	 Smax0=Sre0(i)
c	 Nmax0=i
	endif
18	continue
	AberrL=(1-Smax/Sim(Ns/2))*100	! Loss due to aberration
	AberSp=100-(100-AberrL)*(100-SpillL)*1d-2
	TotL = 100-(100-AberrL)*(100-SpillL)*(100-Gillum)*1d-4
	TotL0= (1-PowerT/Power0*(Smax/Sim(Ns/2)))*100 !
c ~= AberrL+SpillL+Gillum

c	print*,Smax,Smax0

	do 19 i=0,Ns-1
	if(idiff.ge.0) Sim0(i)=Sre(i)-Sre0(i)*Smax/Smax0
19	if(idiff.lt.0) Sim0(i)=Sre(i)-Sim(i)

	if(iabs(idiff).lt.2) go to 17
	SmaxD=-1d9
	SminD=+1d9
	do 16 i=0,Ns-1
	if(Sim0(i).ge.SmaxD) SmaxD=Sim0(i)
16	if(Sim0(i).lt.SminD) SminD=Sim0(i)
17	ComaL=SminD/Smax*100
	ComaH=SmaxD/Smax*100
	end




	subroutine FindMax(x1,y1,x2,y2,x3,y3,xm)!,ym
c Parabolic fit to 3 points near pattern maximum uses this solution:
	implicit real*8 (a-h,o-z)

	xm = ((y1-y2)*x3**2 + (y3-y1)*x2**2 + (y2-y3)*x1**2)
     * *.5/((y1-y2)*x3 + (y3-y1)*x2 + (y2-y3)*x1)
c	ym = (x2-xm)*(x3-xm)*y1/((x1-x2)*(x1-x3))
c     *    +(x1-xm)*(x3-xm)*y2/((x2-x1)*(x2-x3))
c     *    +(x1-xm)*(x2-xm)*y3/((x1-x3)*(x2-x3))
c Parabola ax^2+bx+c coefficients:
c      a=y1/((x1-x2)*(x1-x3))+y2/((x1-x2)*(x3-x2))-y3/((x2-x3)*(x3-x1))
c	b= (x2+x3)*y1/((x1-x2)*(x3-x1)) + (x1+x3)*y2/((x1-x2)*(x2-x3)) 
c     *   + (x1 + x2)*y3/((x1 - x3)*(x3 - x2))
c	c=-x2*x3*y1/((x1 - x2)*(-x1 + x3)) +
c     * x1*x3*y2/((x1 - x2)*(-x2 + x3)) - x1*x2*y3/((x2 - x3)*(-x1 + x3))
	return
c*	end

	entry FindMaxE(y1,y2,y3,ym)
c	entry FindMaxE(x1,y1,x2,y2,x3,y3,xm,ym)
c Parabola maximum from 3 equidistant points:
c	xm = x2 + .5*(x2-x1)*(y1 - y3)/(y1 - 2*y2 + y3)
	ym=(2*y1*(4*y2+y3)-y1**2-(y3-4*y2)**2)/(8*(y1+y3)-16*y2)
	return

	entry FindDer(x1,y1,x2,y2,x3,y3,xm,ym,der1,der2)
c Derivative of parabola through 3 points at the 3rd point:
c Parabola ax^2+bx+c coefficients:
	xm = ((y1-y2)*x3**2 + (y3-y1)*x2**2 + (y2-y3)*x1**2)
     * *.5/((y1-y2)*x3 + (y3-y1)*x2 + (y2-y3)*x1)
	ym = (x2-xm)*(x3-xm)*y1/((x1-x2)*(x1-x3))
     *    +(x1-xm)*(x3-xm)*y2/((x2-x1)*(x2-x3)) 
     *    +(x1-xm)*(x2-xm)*y3/((x1-x3)*(x2-x3))
      a = y1/((x1-x2)*(x1-x3))+y2/((x1-x2)*(x3-x2))-y3/((x2-x3)*(x3-x1))
	b = (x2+x3)*y1/((x1-x2)*(x3-x1)) + (x1+x3)*y2/((x1-x2)*(x2-x3)) 
     *   + (x1 + x2)*y3/((x1 - x3)*(x3 - x2))
c	c=x2*x3*y1/((x1 - x2)*(x1-x3)) +
c     * x1*x3*y2/((x1 - x2)*(x3-x2)) + x1*x2*y3/((x2-x3)*(x1-x3))
c	ym1 = (a*xm+b)*xm+c
c	print*,ym,ym1,ym-ym1 	! the same numbers to ~1d-12
	der2=2*a
	der1=der2*x3+b
	end

	subroutine FitParab(x,y,N,a,b,c,e,ea,eb,ec)
c fits function y=a+b*x+c*x^2 to the data (x,y) by the lsqr method
c N - number of data in the series
c a, b, c - parameters of fitted parabola
c e, ea, eb, ec - sigma of the fit and errors of the parameters
	implicit real*8 (a-h,o-z)
	real*8 x(*),y(*)
c	print'(i4/(8f8.4))',N,(x(i),y(i),i=1,N)
	if(N.lt.3) return
	s1=0
	s2=0
	s3=0
	s4=0
	sy=0
	sy2=0
	xy=0
	x2y=0

	do 1 i=1,N
	xi=x(i)
	yi=y(i)
	xi2=xi*xi
	s1=s1+xi
	s2=s2+xi2
	s3=s3+xi2*xi
	s4=s4+xi2*xi2
	sy=sy+yi
	sy2=sy2+yi*yi
	xy=xy+xi*yi
1	x2y=x2y+xi2*yi

	s24=s2*s4-s3*s3
	s23=s2*s3-s1*s4
	s13=s1*s3-s2*s2
	s12=s1*s2-N*s3
	delta=N*s24+s1*s23+s2*s13
	a=(s24*sy+s23*xy+s13*x2y)/delta
	b=(s23*sy+(N*s4-s2*s2)*xy+s12*x2y)/delta
	c=(s13*sy+s12*xy+(N*s2-s1*s1)*x2y)/delta
	s=0
	if(N.gt.3) s=dsqrt(dmax1(0.d0,sy2-a*sy-b*xy-c*x2y)/(N-3))
	ea=s*sqrt((N*s4-s2*s2)/delta)
	eb=s*sqrt(s24/delta)
	ec=s*sqrt((N*s2-s1*s1)/delta)
c if only linear trend evident
		if(ec.gt.3.*abs(c)) then
	sl=xy-s1*sy/N
	sm=s2-s1*s1/N
	b=sl/sm
	a=(sy-b*s1)/N
	c=0
	s=dsqrt(dmax1(0.d0,sy2-sy*sy/N-b*sl)/(N-2))
	ea=s*dsqrt(s2/N/sm)
	eb=s/dsqrt(sm)
	ec=0
		endif
	e=s
	end



	Subroutine Optimize(No,Nm,varIncr)
c Minimizes parameter numbered Nm by varying paremeter numbered No
	implicit real*8 (a-h,o-z)
	parameter (imax=20)
	dimension x(imax),y(imax)
	common /params/ par(0:79)
cc No = (Nod/10)*5+mod(Nod,10) -> Nod = (No/5)*10+mod(No,5)
	Nod=mod(No,5)+(No/5)*10	! numbers for display
	Nmd=mod(Nm,5)+(Nm/5)*10
	dx=.1d0
	if(varIncr.gt.0) dx=varIncr
	x(2)=par(No)
	y(2)=par(Nm)
	x(1)=x(2)-dx
	par(No)=x(1)
	call RTPattern(phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)

	i=1
		if(par(Nm).lt.y(2)) then
	y(1)=y(2)
	x(1)=x(2)
	i=2
		endif
	x(i)=par(No)
	y(i)=par(Nm)
	dxn=x(2)-x(1)
	der=(y(2)-y(1))/dxn
	print'(2x,a,i2.2,a,i3.3,a,a,i2.2)',  'Iter    Par',Nod,' Par'
     *,Nmd,"_min     diff       y'  ",' incrPar',Nod
	print'(2x,i3,5f10.4)',1,x(1),y(1),     0.,    0.,0.
     * ,2,x(2),y(2),y(2)-y(1),der,dxn
		if(abs(der).lt.1d-9) then
	print*,'  Insensitive at this starting parameter_value!'
	return
		endif

	iend=0
	do 1 i=3,imax-1
	x(i)=x(i-1)+dxn
	par(No)=x(i)
	call RTPattern(phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	y(i)=par(Nm)
	dern=(y(i)-y(i-1))/dxn
	print'(2x,i3,5f10.4,f10.6)',i,x(i),y(i),y(i)-y(i-1),dern,dxn
	if(iend.ne.0) go to 2
	  if(y(i).ge.y(i-1)) then 
	dxn=-dxn/2
	iend=1
	go to 1
	  endif
	if(der.ne.0.and.abs((dern-der)/der).lt.5d-2) dxn=2*dxn
1	der=dern
	i=i-1

c value of derivative at this point
		if(iend.lt.0) then	! just to skip it altogether
c	call FindDer(x(nn(i-2)),y(nn(i-2)),x(nn(i-1)),y(nn(i-1))
c     *,x(i),y(i),xm,ym,d(i),der2)
	call FindDer(x(i-2),y(i-2),x(i-1),y(i-1),x(i),y(i),xm,ym,d,der)
c	call FindMax(x(1),y(1),x(2),y(2),x(3),y(3),xm)!,ym
		endif

2		if(i.lt.imax-1) then  !double waight of 2 (of 4) middle values
	i=i+1
	x(i)=x(i-2)
	y(i)=y(i-2)
	i=i+1
	x(i)=x(i-2)
	y(i)=y(i-2)
		endif
	i0=max(i-5,1)
	call FitParab(x(i0),y(i0),i-i0+1,c,b,a,e,ec,eb,ea)
	if(abs(a).le.1d-7.or.y(i-2).gt.y(i)) then
	 par(No)=x(i)
	 print'(a,i3,a)'," Couldn't find minimum in",imax," iterations"
 	 return
	else
c y = ax^2+bx+c, y'=0=2ax+b -> xmin = -b/(2a)
	 par(No)=-b/(2*a)
	 call RTPattern(phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
c	print'(4f15.6)',par(No),par(Nm),a,b
	endif
	end

	subroutine HelpComm
	print*,' To set new value for a parameter type-in its number'//
     *' followed by space or comma and the value. E.g., to set 0.7504'//
     * 'm for feed leftward offset (in y coordinate)  type "21 .7504"'
     *//' (without quotes), where 21 signifies row (#) 2 and column #1.'
	print*,' "520,5.,.75 3 -.2,+.1" set values for 5'//
     *' consecutive parameters (20 to 24).'
	print*,' "-21,140 [step]" attempt to minimize par. No 140 '//
     *'wrt No 21 [search increment;     default is 0.1 m or deg].'
     *//' Par. No alone: show 6 significant digits of value.'
      print*,' "o" optimize (for maximum gain) offsets in'//
     *' feed''s Z and its pattern squints.'
      print*,' "p #i" pointing and gain change due to par. No #i '//
     * '(e.g. 20 for feed X-offset).'
      print*,' "m [11]" secondary mirror angular shape as seen by feed'
     *//' [11 lateral sections].'
      print*,' "$ [ijk]" save all current parameters into the .par '
     *//'file (read on $ at start);    if any of i, j or k equals 1, ' 
     *//'save amplitudes (i=1), power pattern (j=1) and/   /or residua '
     *//'(k=1) into the .dat, .pat and/or .res file, respectively.' 
      print*,' "a" show classical approximations for '//	!optical properties 
     *'feed offset gain losses and beam rates.'
      print*,' "c" clear all the 10 offsets (par. 1i and 2i,'//
     *' where i assumes values 0 to 4).'
      print*,' "g[n]" write .ps and .txt output files'
     *//' [*0 power, *1 voltage, 0* full, 1* 1%].'
      print*,' "ENTER" evaluate (on change) and display derived'//
     *' parameters for current setup.'
      print*,' "q" quit this session. 		 "?" or "h" print this help.'
c     *//''
		end


	Subroutine ReadPar(iGraph,phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
c This program unit displays most important parameters and
c allows user to interact with it.

	implicit real*8 (a-h,o-z)
	character*40 Descr(0:15),help(0:4,0:15)*49,rec*80,byte*1
     *,fname*80
	real*8 parHere(15:19),parSave(15:19)
	data parHere/45d0,-14d0,-14d0,1.2d0,1.2d0/,iRept/1/
	common /params/ par(0:4,0:15)
	common /output/ ie,in,ip,fname 
	data Descr/
     *'Cassegrain: D,X_min,f,Xs_max,h2   [m]',
     *'Subreflector: x,y,z,  tx, ty  [m,deg]',		!1
     *'Feed offsets: x,y,z, Dtx,Dty  [m,deg]',
     *'Ray-tracing: N/R.NxNy,Nu,Nv,maxU,maxV',
     *'Wavel,Taper,zAprt,pltD,pivot [m,dB,-]',		!4
     *'    Reserved for other param''s/filenames etc; now: ',
c ,xpiv,zpiv,ddu,ddv,dNsh - trzeba: oIncr,fnames,inputFile,NrRT,AberRef
     *'Main dish: f/D,f/Dp, D2, G,coordG [m]',		!6
     *' angles: t_V,t_H, t_0, t_1, t_2 [deg]',
     *'Subreflector: Dc,Dm,Ds2,Gs,Xs_min [m]',
     *' feed angles: t_VH,t_c,_0,_1,_2 [deg]',
     *'Cassegrain: f1, f2, fs, F,   M  [m,-]',		!10
     *'Hyperbola: a,c,b,e,  FSAmpRatio [m,-]',
     *'Beam: t_u,t_v,     HPBW_u,_v,_0 [deg]',
     *'  Squint, phi_X,Y, apert_X,Y  [deg,m]',
     *'Loss: Aberr,Spill,A&S,IllDec,Totl [%]',		!14
     *"  pathR,ComaLow,-Hi,1stLob,Totl'[m,%]"/		!15
	data help/	!(0:4,0:7)
     * ' main dish diameter as projected onto aperture, m'
     *,' smallest dish rim to paraboloid axis distance, m'
     *,' paraboloid focal length, m'
     *,' largest secondary rim to axis distance, m'
     *,' height of secondary focus above dish vertex, m'
     *,' lateral offset of subreflector along x-axis, m'
     *,' lateral offset of subreflector along y-axis, m'
     *,' offset of subreflector along z-axis, m'
     *,' subreflector tilt toward +ve x direction, deg'
     *,' subreflector tilt toward +ve y direction, deg'
     *,' lateral offset of feed along x-axis, m'
     *,' lateral offset of feed along y-axis, m'
     *,' offset of feed along z-axis, m'
     *,' feed pattern tilt toward +ve x direction, deg'
     *,' feed pattern tilt toward +ve y direction, deg'
     *,' No of rays/dish radius + No Feeds X/100+Y/10000'
     *,' No of discrete spatial frequencies in u'
     *,' No of discrete spatial frequencies in v'
     *,' +/- range of patterns in u-spacial frequency'
     *,' +/- range of patterns in v-spacial frequency'
     *,' observing wavelength, m'
     *,' signed feed illumination taper at dish edge, dB '
     *,' z-coordinate of assumed aperture plane, m'
     *,' reference pattern index; <0 for phaseless rays'
     *,' rotation center: -1/0/1 = user/focus/vertex'
     *,' x-coordinate of secondary rotation center, m'
     *,' z-coordinate of secondary rotation center, m'
     *,' main lobe shift in u wrt central ray'
     *,' main lobe shift in v wrt central ray'
     *,' main lobe shift in bins wrt central ray'
     *,' dish focal ratio: f/D'
     *,' F# of parent paraboloid: f/D_parent'
     *," elliptical dish longer diameter, m"
     *,' dish largest depth, m'
     *,' deepest point coordinate wrt dish x-edge, m'
     *,' dish subtended angle in vertical plane, deg'
     *,' dish lateral or horizontal subtended angle, deg'
     *,' central ray prime-focal angle wrt "-z"-axis, deg'
     *,' least x of dish focal angle wrt "-z"-axis, deg'
     *,' greatest dish x focal angle wrt "-z"-axis, deg'
     *,' subreflector size at mid of longest diameter, m '	!80, i.e. (0,8)
     *,' lateral secondary size where dish is widest, m'
     *,' vertical subreflector size (longest diameter), m'
     *,' depth of subreflector at its (longit.) center, m'
     *,' subreflector rim to z-axis (least) distance, m'
     *,' secondary opening cone (subtended) angle, deg'
     *,' z-axis to subreflector opening cone axis, deg'
     *,' z-axis to ray that reaches aperture center, deg'
     *,' smallest inclination of rays to z-axis, deg'
     *,' largest inclination of rays to z-axis, deg'
     *,' parent hyperboloid vertex to prime focus, m'	! 100
     *,' parent hyperboloid vertex to secondary focus, m'
     *,' interfocal distance (2c = f1 + f2), m'
     *,' telescope effective focal length (Mf), m'
     *,' telescope magnification (f2/f1)'
     *,' a-constant of canonical equation of hyperbola, m'
     *,' c = (f1+f2)/2 = fs/2,     c^2 = a^2+b^2,       m'
     *,' b-constant of canonical equation of hyperbola, m'
     *,' eccentricity of parent hyperboloid (e=c/a)'
     *,' amplitude ratio (rim/maximum) due to f-s taper'
     *,' vertical downward beam deviation (u-dim''n), deg'	!120
     *,' lateral leftward beam deviation (v-dimen''n), deg'
     *,' vertical half power beam width (ray-trace), deg'
     *,' lateral half power beam width (ray-trace), deg'
     *,' beamwidth calculated from Baars formula, deg'
     *,' absolute value of main beam offset, deg'
     *,' feed pattern inclination to yz plane, deg'
     *,' feed pattern inclination to xz plane, deg'
     *,' aperture x-coordinate of feed pattern maximum, m'
     *,' aperture y-coordinate of feed pattern maximum, m'
     *,' aberration (coma & astigmatism) loss of gain, %'		!140
     *,' loss of gain due to increased spillover, %'
     *," loss of gain due to aberration and spillover, %"
     *,' loss due to decentered illumination, %'
     *,' aberration, spillover & decentering combined, %'
     *,' range of ray path errors: focus to aperture, m'
     *,' lowest difference in patterns due to offsets, %'
     *,' highest difference in patterns due to offsets, %'
     *,' first side lobe power level (+ve u-axis only), %'
     *,' aberration, spillover & decentering combined, %'
     */

	if(iRept.eq.0) go to 1
	call RTPattern(phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)

3	l=-40	!
7	iRept=0
1	print'(a,T41,a,4(6x,1h#,i1)/1x,3(a39,i1,5f8.3),
     *(a39,i1,f9.4,2(f5.0,2x),1x,2f8.3),(a39,i1,3f8.3,f5.0,f8.0))',
     * '        SETTABLE parameters','#     #0',(i,i=1,4),
     * (Descr(i),i,(par(j,i),j=0,4),i=0,4)
c     *, (Descr(i),i,(par(j,i),j=0,3),i=4,4)

		if(l.lt.0) then	! for first entry into this subprog.
			l=-l
			go to 6
		endif
	go to 2
5	ko=k1
4	print'(a,i4,a)','  The number',ko,' is incorrect!  Try again.'
2	print'(a,i3.2,f9.3,a)','  Press ?,ENTER,g,$,c,m,o,a,p #i, or'
     *//' #i value (e.g.',  l,par(mod(l,10),l/10),'): '

c m - secondary mirror shape
c o - optimize z & tx & ty
c c - clear offsets

	read(*,'(a)') rec			!
		if(rec.eq.' ') rec='S'

c ======1-letter commands=======
	read(rec,*) byte	!,byte
c	print*,' byte: "'//byte//'" rec(1:4): "'//rec(1:4)//'"'
		if(byte.eq.'q') stop
		if(byte.eq.'g') then
	  if(iRept.ne.0) call RTPattern
     *                  (phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	iRept=0
			if(rec(2:10).eq.' ') then
	iGraph=0
			else
	read(rec(2:10),*) iGraph
			endif
	return
		endif
		if(byte.eq.'S') then
	if(iRept.gt.0) call RTPattern
     *                  (phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	go to 6
		endif
		if(byte.eq.'a') then
	call Formulae
	go to 2
		endif
		if(byte.eq.'?'.or.byte.eq.'h') then
	call HelpComm
      print*,'Remark: X_min (No=01) is X-coordinate of the '//
     *'point on dish edge closest to the  Z-axis; to work with '
     *//'common, non-offset Cassegrain antenna set it to -D/2.'
	go to 2
		endif
		if(byte.eq.'c') then	! clear all offsets
	do 15 i=0,9
15	par(i,1) = 0d0
	iRept=iRept+1
	go to 1
		endif
		if(byte.eq.'$') then
	  if(rec(2:5).eq.' ') then	! save params to file
	print*,'       Saving all parameters to "working$.par".'
	if(ip.eq.0) open(19,'working$.par')
	if(ip.gt.0) open(19,fname(1:ip)//'working$.par')
	write(19,*) par
	close(19)
	  else				! write data/results to a file
c	print*,'       Saving data according to command: '//rec(1:7)
	read(rec(2:5),*) iopt
	Ns=(par(1,3)+1)*(par(2,3)+1)	! = (Nfu+1)*(Nfv+1)
	call outputArr(iopt,Ndat,Ns)
	  endif
	go to 2
		endif
		if(byte.eq.'m') then
c	print*,' "'//rec(3:10)//'"'
	  if(rec(3:50).eq.' '.or.rec(2:2).ne.' ') then	!.rec(1:1).eq.'m'.and
	print*,' For subr. shape use "m Nd", where Nd [5] is number'
     *//' of diameters (3 to 101).'
	Ndiam=5
	  else
	read(rec,'(a,i3)') byte,Ndiam
	  endif
	call  SubShape(Ndiam)
	go to 6
		endif

		if(byte.eq.'p') then ! pointing error due to a parameter change
	read(rec,*) byte,k
	 dpar=.01
	par(mod(k,10),k/10)=par(mod(k,10),k/10)+dpar
	call RTPattern(phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	 dtu=par(0,12)
	 dtv=par(1,12)
	 deff=par(2,14)
	par(mod(k,10),k/10)=par(mod(k,10),k/10)-dpar
	call RTPattern(phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	 dtu=(dtu-par(0,12))/dpar
	 dtv=(dtv-par(1,12))/dpar
	 deff=(deff-par(2,14))/dpar
	print'(a,f7.1,3h'' &,f7.1,a,f6.2,a)',"  Pointing errors per m"
     * //" or deg in u (elev) & v; dEff:", dtu*60,dtv*60,"';",
     *deff,'%'
	iRept=0
	go to 2
		endif

		if(byte.eq.'o') then ! optimization procedure
	do 100 i=15,19
	parSave(i)=par(i,0)		! save some parameters
100	par(i,0)=parHere(i)
c	par(15)=45		! points on aperture radius
c	par(16)=14		! number of frequencies (plus u=0) in u (Nu)
c	par(17)=14		! number of frequencies in v (Nv)
c	par(18)=-1.2	! spectrum width in u (Nu)
c	par(19)=-1.2	! spectrum width in v (Nv)
	if(iRept.gt.0) call RTPattern		!App
     *                  (phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	call Optimize(12,70,oIncr)	! feed z-offset vs. Aberr.
	if(dabs(par(0,2)).gt.1d-3) call Optimize(13,74,oIncr) ! Dtx vs. Totl.
	if(dabs(par(1,2)).gt.1d-3) call Optimize(14,74,oIncr) ! Dty vs. Totl.
c	par(2,14)=100-(100-par(0,14))*(100-par(1,14))/100
	print'(2(a,f6.3,3f7.3,f6.3))',
     *'  Minimized Aber,Sp,A&S,Tot[%], and Sq[deg]: ',
     *par(0,14),par(1,14),par(2,14),par(4,14),par(0,13)
	iRept=iRept+1
	do 101 i=15,19
101	par(i,0)=parSave(i)
	go to 1
		endif
c ==============================

	read(rec,*) k1
	  if(k1.lt.0) then	!------Option to minimize ------
	iarg=0
	ifirst=0
	do 10 i=80,3,-1
	if(rec(i:i).ne.' ') go to 11
10	if(rec(81-i:81-i).ne.' '.and.ifirst.le.0) ifirst=81-i
11	i=i-1
	if(rec(i:i).eq.' '.or.rec(i:i).eq.',') iarg=iarg+1
	if(i.gt.ifirst) go to 11
	if(iarg.ge.2) read(rec,*) k1,km,oIncr
	if(iarg.eq.1) read(rec,*) k1,km
	if(iarg.lt.1) go to 5
	ko=-k1
		if(ko.ge.43.or.mod(ko,10).gt.4) go to 4
		if(km.lt.60.or.km.gt.154.or.mod(km,10).gt.4) then
	print'(a,i3,a)','  The Nm(2)',km,' is incorrect!  Try again.'
	go to 2
		endif
	do 102 i=15,19
	parSave(i)=par(i,0)		! save some parameters
102	par(i,0)=parHere(i)
c	print*,' ParSave:', parsave
	if(iRept.gt.0) call RTPattern	!App
     *                  (phifx,phify,du,dv,u1,u2,Nmax,Smax,Ndat)
	call Optimize((ko/10)*5+mod(ko,10),(km/10)*5+mod(km,10),oIncr)
	l=k
	do 103 i=15,19
103	par(i,0)=parSave(i)
	iRept=iRept+1
	go to 1
	  endif			!-------------------------

		if(k1.lt.155) then  !******* Normal input (0-155) ********
	do 20 i=2,78
	if(rec(i:i).eq.' '.or.rec(i:i).eq.',') go to 21
20 	continue
		print*,'   No spaces?!: "'//rec//'"'
21		if(rec(i+1:80).ne.' ') go to 22
	print'(a,i4.2,a,1pG13.6,a)',' Parameter No',k1,':'
     * ,par(mod(k1,10),k1/10), help(mod(k1,10),k1/10)
	go to 2
22	read(rec,*) k,value
	if(k/10.lt.6) par(mod(k,10),k/10)=value !***********************
		else
c character data ???


c option to enter a few consecutive param's after normal one
	  kn=k1/100	! number of consecutive param's
	  k=k1-100*kn	! max par. No = k/10*5+mod(k,10)+kn must be <30
	  kn=min(k/10*5+mod(k,10)+kn,30+0)-(k/10*5+mod(k,10))
c	print*,kn,k,(par(mod(i,10),i/10),i=k,k+kn-1)
	read(rec,*) k1,(par(mod(i,10),i/10),i=k,k+kn-1)
		endif
c	if(k.ne.44) 
	iRept=iRept+1
	l=iabs(k)
 	go to 1

6	print'(7x,a,T43,f7.3,4f8.3,
     * (a38,i2,2f8.4,3f8.3),(a38,i2,5f8.3),
     *2(a38,i2,5f8.4),(a38,i2,3f8.4,f8.3,f8.4)
     *4(a38,i2,5f8.4),(a38,i2,4f8.4,f7.3))',
     * ' COMPUTED param''s  & x,z_piv,du,v,N',(par(j,5),
     *j=0,4), (Descr(i),i,(par(j,i),j=0,4),i=6,15)
c     *, (Descr(i),i,(par(j,i),j=0,3),i=15,15)
c	if(byte.eq.'q') return
	go to 2

	end

	subroutine outputArr(ioption,Ndat,Ns)
	character*80 fname
	common /output/ ie,in,ip,fname 
	common /arrays/ xy(15000,2),dat(15000,2),xyS(1:20301,2)
     *,S(1:20301,2),S0(1:20301,2)
c ioption = ijk,  i=1 - ampl. distr. xy,dat -> *.dat
c			j=1 - power spectrum	  -> *.pat
c			k=1 - pottern differ.	  -> *.res
		if(mod(ioption,10).ne.0) then	! .res
	open(8,fname(1:ie-1)//'.res')
	write(8,'(2f10.4,1f12.1)') ((xyS(i,j),j=1,2),S0(i,2),i=1,Ns)
	close(8)
	print*,' New content written to file '//fname(1:ie)
     & //'res  (to be overwritten on next call!!!)'
		endif
		if(mod(ioption,100)/10.ne.0) then	! .pat
	open(8,fname(1:ie-1)//'.pat')
	write(8,'(2f10.4,1f12.1)') ((xyS(i,j),j=1,2),S(i,1),i=1,Ns)
	close(8)
	print*,' New content written to file '//fname(1:ie)
     & //'pat  (to be overwritten on next call!!!)'
		endif
		if(ioption/100.ne.0) then	! .dat
	open(8,fname(1:ie-1)//'.dat')
	write(8,'(2f10.4,2f9.5)') 
     & ((xy(i,j),j=1,2),(dat(i,j),j=1,2),i=1,Ndat)
	close(8)
	print*,' New content written to file '//fname(1:ie)//'dat '
     & //' (to be overwritten on next call!!!)'
		endif
	end

	subroutine Geometry
	implicit real*8 (a-h,o-z)
	data pi/3.1415926535898d0/

	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
c     *,ddu,ddv,dNsh ! par 18
     *,cox,coy,coz ! 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
     *, pathRange,ComaL,ComaH,FstLob,TotL0 

	zhyp(theta)=a*(c+a*dcos(theta))/(a+c*dcos(theta))-c
c = -(c-a)*(c+a)*cos(theta)/(a+c*cos(theta))
	rhyp(theta)=(c-a)*(c+a)*sin(theta)/(a+c*cos(theta))
	rhohyp(theta)=(c-a)*(c+a)/(a+c*cos(theta))
	FAreaL(x)=(asin(x) + x*sqrt(1-x**2))*2/pi

	r2 = r1+D		! outermost edge of dish separation from paraboloid axis
	R0 = r1+D/2		! aperture center distance from main axis
	f_D = f/D		! focal ratio
	fbyD = f/(r2*2)	! focal ratio of parent paraboloid


 	z1 = r1**2/(4*f) - f
 	z2 = r2**2/(4*f) - f
	D2  = dsqrt((z2-z1)**2+D*D)	!     dish length
c E is inclination of plane of the dish elliptical rim
	sinE = (z2-z1)/D2
	cosE = dsqrt(1-sinE**2)
	tanE = (z2-z1)/D

	theta0=2*datan(R0/(2*f))	! center of aperture angle
	theta1=2*datan(r1/(2*f))	! subtanding angles in vertical
	theta2=2*datan(r2/(2*f))	!	plane
	t_0=theta0*180/pi
	t_1=theta1*180/pi
	t_2=theta2*180/pi
	Depth=(f*tanE-r1)*sinE+r1**2*cosE/(4*f)
	xDepth = f*(1 + 1/(cosE*(1+cosE)))*sinE - r1*(cosE+r1*sinE/(4*f))
c which is equal to D2/2 if E=0 (symmetrical telescope)

	zc = (z1+z2)/2		! center defined by middle of lateral diameter D
	rc = 2*dsqrt((zc+f)*f)	! projected onto paraboloid
c this is equivalent to taking zc = [(D/2)^2+((r1+D/2)^2]/(4f) - f
cc	rc=R0				! rc,zc that would lead to narrower
cc	zc = R0**2/(4*f) - f	! 	central diameter/size
	theta3 = 2*datan(dsqrt(zc/f+1))
c = 2*datan(rc/(2f))		! horizonthal dish tips angle
	rhoc=dsqrt(zc*zc+rc*rc)		! distance focus-dish side extreme point
	thetaH=2*dasin(D/(2*rhoc))	! dish subtended angle (horizontally)
	thetaV=theta2-theta1 		! dish subtended angle (vertically)

	t_H=thetaH*180/pi
	t_V=thetaV*180/pi
c	t_0=(theta1+theta2)/2*180/pi
c	print'(a,T30,4f8.3)','  r1 z1 theta1, theta2:', r1,f+z1,t_1,t_2
c	print'(a,T30,10f8.3)',' PARABOL  D D2 f xDepth Depth:',
c     * D, D2, f, xDepth, Depth
c	print'(a,T30,10f8.3)','  thetaV thetaH thetaC:' ,t_V,t_H,t_0

c Compute Cassegrain parameters (hyperboloidal subreflector):
	c = (f - hf2)/2	! hyperbola c-parameter (half foci separation)
c Distance of outer tip of subreflector from primary focus:
	rhoh2 =r2s*(f/r2+.25*r2/f)	! = r2s/dsin(theta2)

c Flat mirror par0-4: 56.000   4.000  60.000   4.000  54.000
c Since rhoh2 = (c**2-a**2)/(a+c*cos(theta2)):
	delta = rhoh2**2 - 4*(rhoh2*dcos(theta2) - c)*c
	a = (dsqrt(delta) - rhoh2)/2	! hyperbola a-parameter 
	b2=(c**2-a**2)
c thus the parent hyperboloid vertex is
	f1 = c-a	! away from primary focus and
	f2 = c+a	! away from secondary focus;
	fs=f1+f2	! this is interfocal distance (= 2*c)
	e=c/dmax1(1d-10,a)	! eccentricity
	b=dsqrt(b2)
c Telescope magnification and effective focal distance are:
	amagn = f2/f1
	Feff = f*amagn

	Aratio = (1+dmax1(r2s,r1)**2/(2*Feff)**2)/(1+(r2/(2*Feff))**2)
c	    = cos^2(phi/2)/cos^2(phi'/2) free-space taper; amplitude
	fi1=atan(r1/(2*Feff))*360/pi
	fi2=atan(r2/(2*Feff))*360/pi
	fiC=(fi1+fi2)/2
	fiV=fi2-fi1

c coordinates of subreflector pivot (of tilt, relative to paraboloid 
c frame; in x-z plane only)
c	theta0=2*datan(R0/(2*f))	! aperture center angle
c	rhohc=rhohyp(theta0)
		if(idint(pivot).eq.1) then	! subrefl. vertex
	xpiv=rhyp(theta0) 
	zpiv=zhyp(theta0)
		endif
		if(idint(pivot).eq.0) then	! prime focus as center of rotation
	xpiv=0
	zpiv=0
		endif
c The other (innermost) end of subreflector rim is 
	rhoh1 = (c**2-a**2)/(a+c*dcos(theta1))	! away from the focus
	if(abs(rhoh1-rhoh2).lt.1d-8) rhoh1=rhoh2	! for symmetrical antenna
c Longer diameter of the subreflector:
	Dh2 = dsqrt(rhoh1**2+rhoh2**2 - 2*rhoh1*rhoh2*dcos(thetaV))
c Coordinates of this diameter center
	r1s=rhoh1*dsin(theta1)
	rhc=(r1s+r2s)/2		! r2s=rhoh2*dsin(theta2)
	zhc=-(rhoh1*dcos(theta1)+rhoh2*dcos(theta2))/2
c Distance from focus to the above center:
c	rhohc=dsqrt(rhc**2+zhc**2)	! the same is obtained with:
c	rhohc = dsqrt((rhoh1-rhoh2)**2/4 
c     * +rhoh1*rhoh2*dcos((theta2-theta1)/2)**2)
c The perpendicular diameter of the subreflector rim at the same height
c ends sqrt(rchsqr) away from the parent hyperboloid axis.
c It however does not correspond to the largest horizontal diameter
c of the main dish; the one corresponding to it is this long:
	Dhm = 2*rhohyp(theta3)*dsin(thetaH/2)
c which derives from: Dhm**2=2*rho**2-2*rho*rho*cos(thetaH)
		if(a.gt.1d-6) then
	rchsqr = ((zhc+c)**2-a**2)*(c**2/a**2-1)
c Length of mentioned lateral diameter at zhc is
	Dhc = 2*dsqrt(rchsqr-rhc**2)
		else	! flat subreflector
	Dhc=Dhm
		endif
c Find subreflector depth, g, at physical (as opposed to optical) center;
c eqs. to solve are: (z+c)^2/a^2-(x^2+y^2)/b^2=1; z=zhc-g cosz;  r=rhc+g cosr
	angle=datan2(-rhoh2*cos(theta2)-zhc,r2s-rhc)
	aa=(c*cos(angle))**2-a**2		!b2
	bb=2*(-b2*cos(angle)*(zhc+c)-a**2*rhc*sin(angle)) !-2b^2(zhc+c)
	cc=b2*zhc*(zhc+2*c)+b2**2-(a*rhc)**2 
	g=(-bb-dsqrt(dmax1(0d0,bb*bb-4*aa*cc)))/(2*aa)
c	g2=(-bb+dsqrt(bb*bb-4*aa*cc))/(2*aa)
c	print*,' g, g2:',g,g2,angle*180/pi,zhc+c-sqrt(zhc*c+a**2),rhc,zhc

	end

	subroutine Formulae
	implicit real*8 (a-h,o-z)
	data pi/3.1415926535898d0/
	equivalence (par(0,0),D),(par(2,0),f),(par(3,0),r2s)	!,(par(1,0),r1)
c     *,(par(4,0),hf2),(par(0,1),xsub),(par(1,1),ysub),(par(2,1),zsub)
c     *,(par(8),subtiltX),(par(9),subtiltY)
     *,(par(0,2),xoff),(par(1,2),yoff),(par(2,2),zoff)
c     *,(par(13),tiltfx),(par(14),tiltfy),(par(15),datNr)
c     *,(par(16),parNfu),(par(17),parNfv),(par(18),freqmx),(par(19),freqmy)
     *,(par(0,4),wavel),(par(1,4),TaperdB), (par(0,8),Dhc)
     *,(par(3,10),Feff),(par(4,10),amagn),(par(1,11),c)
     *,(par(0,5),xpiv),(par(1,5),zpiv),(par(0,10),f1)
	common /params/ par(0:4,0:15)
c Beam deviation function is defined as {r=R/(D/2), t=1-taperP}: 
c Integ[0,1]{(1-(1-t)r^2)r^3/[1+(r/(4f/D))^2]}dr/
c  /Integ[0,1]{(1-(1-t)r^2)r^3}dr which using F=(4f/D)^2 reduces to:
c	BDFac(t,FD) = 3*FD*(1 + t + 2*FD*(1 - t) + 
c     & 2*FD*(FD*(t - 1) - 1)*dLog(1/FD + 1))/(1 + 2*t)
	BDFact(t,fD)=fD*(2*(t*fD+1)*(1-fD*dLog(1/fD+1))-t)/(1-2*t/3.) !!! t=1-tau
	tau=10**(-dabs(taperdB)/20)

		if(taperdB.lt.0) then
	Alfa=taperdB*0.1151292546497d0	!.115...=log(10)/20
c	effIlL=1+2*(1-exp(Alfa))**2/(Alfa*(1-exp(2*Alfa)))
	effIlL=1+2*(1-exp(Alfa))/(Alfa*(1+exp(Alfa)))
	exp_1=exp(-Alfa)-1
	AstCof = 8/Alfa**2-(4/exp_1 + (6 - 4/Alfa))/exp_1
	ComaCof=-24/Alfa**3 - (4 - 12/Alfa*(1 - 2/Alfa))/exp_1
c	print*,alfa,astcof,comacof
		else
	TauP1=1+tau
	effIlL=1-.75d0*TauP1**2/(TauP1+tau**2)
c	AstCof= (11 - 4/(Tau-2)**2 + 2/(Tau-2))/9.	! here Tau=1-tau
c	ComaCof= 1.6 + 1.2/(Tau-2)
	AstCof= (11 - (4/TauP1 + 2)/TauP1)/9.
	ComaCof= 1.6d0 - 1.2d0/TauP1
		endif

		if(taperdB.lt.0) then 	! Gaussian illumination function
	BDF=bdfG(f/D,-taperdB)				! beam deviation factor
	BDFs=bdfG(Feff/D,-taperdB)			! BDF of Cassegrain
		else				! parabolic illumination function
	BDF=BDFact(1-tau,(4*f/D)**2)			! beam deviation factor
	BDFs=BDFact(1-tau,(4*Feff/D)**2)		! BDF of Cassegrain 
c	BDF=BDFac(tau,(4*f/D)**2)
c	BDFs=BDFac(tau,(4*Feff/D)**2) 
		endif
	rad2deg=180/pi
c According to JWM Baars 2003 (ALMA Memo 456) [some signs changed!]
	TranslPrim = -BDF/f*rad2deg !*xoff			! deg/m	-
	RotPrim = +(1+BDF)	    !*angle			! deg/deg	+
	TranslSec = -(BDF-BDFs/amagn)/f*rad2deg !*xoff	! deg/m	+
	RotSecV = +2*c/f*(BDF+BDFs)/(amagn+1) !*angle	! deg/deg	-
	RotSecF = +2*c/f*BDFs/amagn !*angle			! deg/deg	-
c-	RotSec_z = -(zpiv-(zpiv+2*c)/amagn)/f		! deg/deg (from NRAO MEMO246)
	RotSec_z = ((2*c+zpiv)*BDFs/amagn-zpiv*BDF)/f	! deg/deg (KMB)
	TranslFeed = -BDFs/(amagn*f)*rad2deg !*xoff	! deg/m	+

	SideLobe_dB=-24.7+7.9*tau	! first sidelobe in dB [-24 to -17 dB]
c accuracy: a few tenths of decibel for 0<tau<0.8
	PowerLevel=10**(SideLobe_dB/10)*100			! %
	bBaars=1.243+(-.343+.12*tau)*tau		! 0.1<tau<1
	HPBW	= bBaars*wavel/D
	HPBW_Baars = HPBW*rad2deg			! in degrees
	feedOff_HPBM = -HPBW/BDF*f			! xoff/HPBW		+
	feedOff_HPBMc = -HPBW/BDFs*Feff		! xoff/HPBW Cassegr. +

c Theoretical aberration losses in aperture efficiency
c R. Padman: "Optical fundamentals for array feeds," Multi-Feed
c Systems for Radio Telescopes, (D. T. Emerson and J. M. Payne, Eds),
c ASP Conference Series, vol. 7, 1995.
	xoff=par(0,2)
	yoff=par(1,2)
	Dhc=par(0,8)
	alpha = dsqrt(xoff**2+yoff**2)/Feff 	! beam squint on the sky
	piD=(pi*D**3*alpha/(16*f*wavel))**2
	deta_curv=piD*(alpha/Dhc)**2/3
	deta_astig=3*deta_curv/amagn**2 *AstCof
	deta_coma=.5*piD/(6*amagn**2*f)**2 *ComaCof
	gainl=(deta_astig + deta_coma)*100
c total loss after corr. for curv. %
c P-S. Kildal: "The effects of subreflector diffraction on the aperture
c efficiency of a conventional Cassegrain antenna—An analytical approach,"
c IEEE Trans. Antennas Propagat., vol. AP-31, no. 6, pp. 903-909, Nov. 1983.
c  A0=tau=1-taperP  amplitude taper (12 dB) [Gaussian illumination assumed] .25
		if(1.-tau.lt..00001d0) then
	deta_difr=+2/pi*tau*dsqrt(wavel/Dhc)*100
		else
	deta_difr=-2*dlog(tau)*tau/(pi-pi*tau)*dsqrt(wavel/Dhc)*100
		endif
c Radius of curvature of the Petzval surface (JW Lamb, MMA Memo 246)
c	R_Petz = f*2/D*r2s	! not good in general
	R_Petz = f*f1/(f+f1)
c	par(4,5)=f*f1/(f-f1)

	write(6,'(T8,a)')' Some characteristics valid for symmetric'
     *//' reflector antenna'
	write(6,'(1x,a/)')' (classical optics approx.; '//
     *'parabolic on pedestal or Gaussian illumin.)'
	write(6,'(a,T51,f9.4,a)')
     *  ' Half power beam width (Baars 2003) ..............',HPBW_Baars,
     &'          [deg]'
     *, ' First side lobe power level due to taper (Baars) ',PowerLevel,
     &'            [%]'
     *, ' Gain loss due to tapered illumination ...........',
     &effIll*100,'            [%]'
     *, ' Gain loss due to diffraction at subrefl. (Lamb) .',deta_difr,
     &'            [%]'
     *, ' Astigmatic gain loss due to lateral feed offset .',
     &deta_astig*100,'            [%]'
     *, ' Coma gain loss due to lateral feed offset .......',
     &deta_coma*100,'            [%]'
     *, ' Aberration (astigmatism & coma) gain loss .......', gainl,
     &'            [%]'
cNot useful thing!     *, ' Gain loss due to curvature ......................',
c     &deta_curv*100,'            [%]'
     *, ' Petzval surface radius (inverse field curvature) .', R_Petz,
     &'            [m]'
c     *, ' Beam deviation factor, prime focus ..............', BDF,
c     &' dimensionless'
c     *, ' Beam deviation factor, secondary focus ..........', BDFs,
c     &' dimensionless'
     *, ' Beam dev. due to primary translation (Baars) .....',
     & TranslPrim,'*Xoff/m   [deg]'
     *, ' Beam dev. due to primary rotation ...............', RotPrim,
     &'*tilt/deg [deg]'
     *, ' Beam dev. due to secondary translation ..........', TranslSec,
     &'*Xoff/m   [deg]'
     *, ' Beam dev. due to secondary rotation about vertex', RotSecV,
     &'*tilt/deg [deg]'
     *, ' Beam dev. due to secondary rotation about focus .', RotSecF,
     &'*tilt/deg [deg]'
     *, ' Beam dev. due to secondary rotation about "z_piv"', RotSec_z,
     &'*tilt/deg [deg]'
     *, ' Beam dev. due to feed translation in sec. focus .', TranslFeed
     &,'*Xoff/m   [deg]'
     *, ' Feed offset per beam width, prime focus (Baars) .',
     &feedOff_HPBM,'*NoOfBeams  [m]'
     *, ' Feed offset per beam width, secondary focus .....',
     &feedOff_HPBMc,'*NoOfBeams  [m]'
     *, ' '

	end

	real*8 function bdfG(fD,dB)
c Computes the beam deviation factor for Gaussian illumination
c function: 
c Integr_0^1 f(r) r^3/[1+(r/4fD)^2] dr/Integr_0^1 f(r) r^3 dr.
c fD = f/D, i.e. focal distance to diameter ratio
c dB = edge illumination taper in deciBells

	implicit real*8(a-h,o-z)
	C1=0.1151292546497d0*dB			!.115...=log(10)/20
	C2=1.842068074395d0*dB*fD**2
	BDFg=(0.2302585092994d0*dB*fD**2*(-8. + 8.*dExp(C1)
     & + dB*dExp(C1 + C2)*fD**2 *14.73654459516d0*
     & (Ei(-C2) - Ei(-C1 - C2))))/(-1. + 10**(dB/20.) - C1)
	end

	real*8 function Ei(Xn)
c Approximates exponential function for Xn <= 0 with accuracy
c better than 1E-5 (2.3E-7 for 0 >= Xn >= -1)
	implicit real*8(a-h,o-z)
	x=-Xn
		if(x.le.1d0) then
	Ei = +0.00107857d0*(0.67635510460d0 - x)*
     &  (31.7221095234d0 - 9.54127532937d0*x + x**2)*
     &  (24.9432440739d0 + 1.168574739854d0*x + x**2) + dLog(x)
		else
	Ei = -(0.250621d0 + 2.334733d0*x + x**2)/
     &    ((1.681534d0 + 3.330657d0*x + x**2)*dExp(x)*x)
		endif
	end



	subroutine Amplitudes(Power0,PowerT,Ndat)
	implicit real*8 (a-h,o-z)
	real*4 xd(15000),yd(15000),datAmp(15000),datPha(15000)
     *,xs(0:20300),ys(0:20300)
     *,Sre(0:20300),Sim(0:20300),Sre0(0:20300),Sim0(0:20300)
	parameter (MaxSide=31, Maxfeed=MaxSide*MaxSide)
	common /feeds/ feedx(Maxfeed),feedy(Maxfeed),feedz(Maxfeed),
     *feedAmp(Maxfeed),feedPh(Maxfeed)

	data pi/3.1415926535898d0/

	common /arrays/ xd,yd,datAmp,datPha,xs,ys,Sre,Sim,Sre0,Sim0
	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
c     *,ddu,ddv,dNsh ! par 18
     *,cox,coy,coz ! 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
     *, pathRange,ComaL,ComaH,FstLob,TotL0 

	zhyp(theta)=a*(c+a*dcos(theta))/(a+c*dcos(theta))-c
c = -(c-a)*(c+a)*cos(theta)/(a+c*cos(theta))
	rhyp(theta)=(c-a)*(c+a)*sin(theta)/(a+c*cos(theta))
	rhohyp(theta)=(c-a)*(c+a)/(a+c*cos(theta))
	FAreaL(x)=(asin(x) + x*sqrt(1-x**2))*2/pi

	Nsidex=min(MaxSide,idint(dmod(.1d0+dabs(datNr)*1d2,1d2)))
	Nsidey=min(MaxSide,idint(dmod(.1d0+dabs(datNr)*1d4,1d2)))
	if(Nsidey.eq.0) Nsidey=Nsidex
	call FeedOff(Nsidex,Nsidey,Wavel)
	datNr=dsign(idint(dabs(datNr))+Nsidex/100.+Nsidey/1d4,datNr)
	Nfeeds=Nsidex*Nsidey

	xof=xoff
	yof=yoff
	zof=zoff

	idiffer=idint(dabs(differ))
c Parabolic on pedestal and Gaussian illumination functions constatnts
	taperP = 1-10**(-dabs(TaperdB)/20.)  !used in: 1-taperP*(2*r/D)**2
	taperG = TaperdB*.115129255d0   !*ln10/20; in: exp(taperG*(2*r/D)**2)
c fHPBW - spatial frequency at half power of symmetrical beam
c HPBW of symmetrical (nonoffset) beam according to JWM Baars [2003,
c ALMA Memo 456] is tau*(0.12*tau-0.343)+1.243, where tau=1-taperP
c In our notation, half of this is [in rad]:
	fHPBW=(taperP*(0.12*taperP+0.103)+1.02)/2
	fHPBWd=2*dasin(fHPBW*wavel/D)*180/pi

	call Geometry

	R0 = r1+D/2		! aperture center distance from main axis


	phi0=2*datan(R0/(2*Feff)) ! feed nominal angle for aperture centre
	phi0d=phi0*180/pi


	dR=D/(2*idint(datNr)+1)	! array step for initial aperture amplitude distr.
	dfX=dtan((phi0+tiltfX/180*pi)/2)*2*Feff-R0
	dfY=dtan(tiltfY/180*pi/2)*2*Feff

c	dfX= R0*(1-dsqrt(1-tiltfX/100d0))	! for Taper=1-.75*(Ri*2d0/D)**2
c	dphi = 2*datan((R0+dfX)/(2*Feff)) - phi0	
c	if(tiltfX.ne.0.) print'(a,4f10.3)',' Extra feed tilts [deg m]:'
c     *, tiltfX,tiltfY, dfX,dfY

	FST0=1/(1+(R0/(2*Feff))**2)	! free-space taper at aperture center
	Power0=0
	PowerT=0
	PowerL=0
	Ll=0
	Ls=0
	pathMin=+1d6
	pathMax=-1d6

	call RayTrace(R0,0d0, xp,dfY0, cox,coy,coz,
     * cphifx,cphify,cphifz, Rdish0,Xapert,Yapert,pathl)

c Loops computing the amplitude distribution and phase on aperture
c plane according to traced rays from initially regular array

	Ndat=0	! number of sampled points counter
	NRdat=datNr
	Xblock=(r2s+r1s)/2
	Rblock=(r2s-r1s)/2

	do 100 i=-NRdat,NRdat 		! ray-tracing of all the rays
	Xi=i*dR 				! original x coordinates (in m)
	Nxi=dsqrt((D/2)**2-Xi**2)/dR	! 2*Nxi+1 is the number of y-points at Xi
c	Nxi=dsqrt((idint(datNr)+.5)**2-i**2)

	do 100 j=-Nxi,Nxi
	Yi=j*dR 				! original y coordinates (in m)
c these (Xi,Yi) are coordinates on physical circular aperture

	rp=dsqrt((R0+Xi)**2+Yi**2)	! aperture point distance from z axis
c	if(rp.lt.r2s) go to 100!	! skip region blocked by subreflector

		if(rp.eq.0d0) then
	fih=0
		else
	fih=datan2(Yi,R0+Xi) ! angle between central and current meridian
		endif

	Xap=0
	Yap=0
	cAmp=0
	sAmp=0
	sumAmp=0
	do 90 NrFeed = 1,Nfeeds

	xoff=xof+feedx(Nrfeed)
	yoff=yof+feedy(Nrfeed)
	zoff=zof+feedz(Nrfeed)

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

	rpblock=dsqrt((R0+Xapert-Xblock)**2+Yapert**2)
	if(rpblock.lt.Rblock) go to 100!	! skip feed region (below subreflector)

	phase=dmod(pathl,wavel)/wavel*2*pi + feedPh(Nrfeed)	! phase of this ray
	if(pathl.lt.pathMin) pathMin=pathl
	if(pathl.gt.pathMax) pathMax=pathl

c Compute amplitude distribution on the aparture
c	rp=sqrt((R0+Xapert)**2+Yapert**2)  ! distance from z-axis
c do not use the above since parallel shift does not change weigting
c phihalf=datan(rp/(2*Feff))	! half phi
c Ampl=dcos(phihalf)**2		! following equivalent avoids trigonometry
	Ampl=FST0/(1+(rp/(2*Feff))**2)	! =cos^2(phi/2) free-space taper
	Rt2=(Xi-dfX)**2+(Yi-dfY)**2 ! at feed Ampl corresp. to, Xi not Xappert
	Ri2=Xi**2+Yi**2
		if(taperdB.ge.0.d0) then
	Taper=1-taperP*Rt2*(2d0/D)**2 ! feed illumination function (ampl.)
	Taper=dmax1(Taper,0d0)
	Taper2=1-taperP*Ri2*4d0/D**2
		else
	Taper=dexp(taperG*Rt2*(2d0/D)**2)
	Taper2=dexp(taperG*Ri2*(2d0/D)**2)
		endif
	AmplTaper=Ampl*Taper*feedAmp(Nrfeed)
	cAmp=cAmp+AmplTaper*dcos(phase)
	sAmp=sAmp+AmplTaper*dsin(phase)
	sumAmp=sumAmp+feedAmp(Nrfeed)
	Xap=Xap+Xapert*feedAmp(Nrfeed)
90	Yap=Yap+Yapert*feedAmp(Nrfeed)

	AmplTaper=dsqrt(sAmp*sAmp + cAmp*cAmp)/Nfeeds	!!!!!!!!(2*Nfeeds+1)
	phase=datan2(sAmp,cAmp)
	PowerT=PowerT+(AmplTaper)**2	! summing total power
		if(Rdish.gt.D/2) then	
c Additional condition: .or.r2s**2.gt.xp**2+yp**2 does not work properly
c because of too few rays sampled near the subreflector; in fact it represents
c loss due to shadowing by the subreflector of sky rays coming on its opposite
c side. A separate correction will be made at the end of this routine.
	PowerL=PowerL+(AmplTaper)**2 ! power spilled over dish edge
	Ampl=0
	Ll=Ll+1
		endif
	Ndat=Ndat+1
	xd(Ndat)=Xap/sumAmp
	yd(Ndat)=Yap/sumAmp
		if(r2s**2.gt.xp**2+yp**2) Ls=Ls+1	
	datAmp(Ndat)=AmplTaper		! ray amplitude
	datPha(Ndat)=phase		! 	& phase on aperture
c The same for symmetrical illumination (taper function)
	xs(Ndat)=Xi		! temporarily (for use in loop to label 200)
	ys(Ndat)=Yi
c The following Ampl2 gives reference to symmetrical paraboloid
	Ampl2=Taper2/(1+Ri2/(2*Feff)**2) ! =Taper2*cos^2(phi/2)=taper*FS_taper

	if(Nfeeds.gt.1) Ampl2=1.	! reference to uniform illumin.

c	if(Ampl.le.0d0) Power0L=Power0L+Ampl2**2
	Power0=Power0+Ampl2**2	! main lobe total power
	if(idiffer.eq.1) datAmp(Ndat) = Ampl2
	if(idiffer.ge.2) Sim0(Ndat)  = Ampl2
c	if(Ndat.lt.3) print*,' Taper*A,taperP*A',Taper*Ampl,Ampl2
c     *,tiltfX
100	continue
	SpillL=PowerL/PowerT*100
	pathRange=pathMax-PathMin

	xoff=xof
	yoff=yof
	zoff=zof

	if(r1.gt.-D/2+0.001) return	! correction only for symmetric telescopes

c ===========Trace the central ray=============
c*	call RayTrace(R0,0d0, xp,dfY0, cox,coy,coz,
c*     * cphifx,cphify,cphifz, Rdish0,Xapert,Yapert,pathl)
c cox/y/z are direction cosines of central ray, dish to aperture plane
c cphifx/y/z are feed-subreflector_center direction cosines

c	phifx=dasin(cphifx)*180/pi
c	phify=dasin(cphify)*180/pi
cc	dfX0=xp-R0
c	dfY0=yp

c Fraction of area sticking out on one side of two equal circles of diameter
c D and centers s apart: FArea=[asin(x) + x*sqrt(1-x**2)]*2/pi, where x=s/D
c
cc	x=dsqrt(dfX0**2+dfY0**2)/D
cc	FAL=FAreaL(x)
cc	FAL2=FAreaL(Rdish0/D)
cc	Ndat2=pi*D**2/(4*dR**2)
cc	print*,'  Ndat,Ll,Ndat1*FAL:',Ndat,Ll,Ndat2,Ndat2*FAL,Ndat*FAL2
cc     * ,Ls,datLs
c The above approximation for outer part of dish does not work because
c of assymmetry due to the free space taper.

c Fractional loss of rays due to shadowing of the subreflector
	datLs=FAreaL(dmin1(Rdish0/(2*r2s),1d0))*pi*(r2s/dR)**2
	rm=r2s+Rdish0/2				! middle of the area for weight
	Ampl=FST0/(1+(rm/(2*Feff))**2)	! =cos^2(phi/2) free-space taper
		if(taperdB.ge.0.d0) then
	Taper=1-taperP*(rm*2d0/D)**2 !
		else
	Taper=dexp(taperG*(rm*2d0/D)**2)
		endif
	PowerLs=datLs*(Ampl*Taper)**2	! effectively lost power
	SpillL=(PowerL+PowerLs)/PowerT*100

	end

	subroutine FeedOff(Nsidex,Nsidey,Wavel)
	implicit real*8 (a-h,o-z)
	character*300 fname*80,rec
	real RT32Nulls(10)
	parameter (MaxSide=31, Maxfeed=MaxSide*MaxSide)
	real*8 feedpar(Maxfeed,5)
	common /feeds/ feedx(Maxfeed),feedy(Maxfeed),feedz(Maxfeed),
     *feedAmp(Maxfeed),feedPh(Maxfeed)
	equivalence (feedx(1),feedpar(1,1)),(feedy(1),feedpar(1,2)),
     *(feedz(1),feedpar(1,3)),(feedAmp(1),feedpar(1,4)),
     *(feedPh(1),feedpar(1,5))
	common /output/ ie,in,ip,fname 
c	PathErr(dr,R,Phi)=-dr*R/(amagn*f+R*R/(4*amagn*f))*dcos(Phi)
	data RT32Nulls/1.216,2.265,3.185,4.290,5.122,6.4,7.1,8.5,
     *8.95,10.1/,ifirst/0/
14			if(ifirst.eq.0) then
13	ifirst=ifirst+1
	print*,'  Offsets, amp''s and phases of feed array required.'
     *,' Type filename or press ENTER (to generate example): '
	read(*,'(a)') rec
	if(rec.eq.' ') then
c		Nside=max(2,min(MaxSide,Nside))
	Nside=Nsidex
		width=Nside*Wavel*.50
		if(Nside.ge.2) then
	dw=width/(Nside-1)
	shift=width/2
		else
	shift=0
		endif
	do 1 i=1,Nside
	dx=(i-1)*dw-shift
	do 1 j=1,Nside
	dy=(j-1)*dw-shift
	feedx(Nside*(i-1)+j) = dx
1	feedy(Nside*(i-1)+j) = dy

	do 2 i=1,Nside*Nside
	d=dsqrt(feedx(i)**2+feedy(i)**2)/Wavel
c	if(dmod(d,1d0).gt..5d0) then
c	 feedPh(i)=3.141592653589d0
c	else
c	 feedPh(i)=0
c	endif
	u=0.5873*3.141592653589d0/180*d*32	!0.5873*pi/180*d[m]*D/Wavel
		if(u.gt.10.5) then
	iPh=u			!dmod(u+1,2d0)
		else
	do 3 iPh=0,9
	if(RT32nulls(1+iPh).gt.u) go to 4
3	continue
	iPh=iPh-1
		endif
4	feedPh(i)=3.141592653589d0*mod(iPh,2)
c
	 feedPh(i)=3.141592653589d0*mod(idint(2*d),2)
c	feedz(i)=u
2	continue
c Nulls in pattern of uniform distr. for RT32 in u-units [u=(D/Wavel)sinTheta]:
c 1.216  2.265  3.185  4.290  5.122  6.4  7.1  8.5 8.95 ~10.1
cBeam dev. due to feed translation in sec. focus:  -0.5873*Xoff/m   [deg]
c	sinTheta=0.5873*pi/180*d[m] = u/(D/Wavel) -> u=0.5873*pi/180*d[m]*D/Wavel

c	print'(5(1x,F6.3,f6.3))',(feedz(i),feedPh(i),i=1,Nside*Nside)
	do 6 i=1,Nside*Nside
c6	feedz(i)=0
6	feedAmp(i)=1.
		else ! read feeds data from the file typed
c	print*,' ip =',ip, '; '//fname(1:20)
	if(ip.gt.0) open(7,fname(1:ip)//rec(1:12))
	if(ip.eq.0) open(7,rec(1:12))
10	read(7,'(a)',end=11) rec
	go to 12
11	print*,' Cannot find "Array size" in this file!!!'
	close(7)
	go to 13
12	if(rec(1:10).ne.'Array size') go to 10
	read(rec(11:80),*) Nsidex,Nsidey
	read(7,*) rec	! skip the header
	do 15 i=1,Nsidex
	do 15 k=1,5
	read(7,'(a)') rec
15	read(rec(6:300),*) (feedpar(Nsidey*(i-1)+j,k),j=1,Nsidey)
	close(7)
		endif
			else ! on second and later calls to this routine
20	print*,' Warning: This is feed array mode!'//
     *' Type A to change the array.'
	print*,' To modify a parameter of a feed type x,y,z,a or p '
     *//'followed by N(row=x),N(col=y) and the value '//
     *'(e.g. "a 2 3,.512") ["ENTER" to quit]: '
21	read(*,'(a)') rec

	  if(rec(1:1).eq.'A') then
		ifirst=0
		go to 14
	  endif
	if(rec(1:1).eq.' ') return
	read(rec(2:50),*) Nx,Ny,value
	if(rec(1:1).eq.'x') feedx((Nx-1)*Nsidey+Ny)=value
	if(rec(1:1).eq.'y') feedy((Nx-1)*Nsidey+Ny)=value
	if(rec(1:1).eq.'z') feedz((Nx-1)*Nsidey+Ny)=value
	if(rec(1:1).eq.'a') feedAmp((Nx-1)*Nsidey+Ny)=value
	if(rec(1:1).eq.'p') feedPh((Nx-1)*Nsidey+Ny)=value
	print'(a,2i3,f6.3,a)',' Accepted: '//rec(1:1),Nx,Ny,value,
     *'.  Next entry or "ENTER" to quit: ' 
	go to 21
			endif
	end
