#!/opt/local/stow/python-2.7.1/bin/python import pylab as pb import numpy as np from math import * import scipy, numpy, os import scipy.special from PIL import Image import pyfits w = 600 h = 600 plateScale = 0.1 # arcsec/pixel fov = w*plateScale # arcsec imageUniform = scipy.zeros((h, w)) imageParabolic = scipy.zeros((h, w)) print "Image scale factor = %.3f arcsec/pixel" % (plateScale) wavelengthmm = 1.0 fnumber = 8.0 diameter = 12.0 scalex = fov * (pi*1000*diameter*pi/(3600*180.*wavelengthmm)) scaley = h * scalex / w print 'scalex = ', scalex print 'v = %f to %f' % ( (.5/w-.5)*scalex, ((w-.5)/w - 0.5)*scalex) # make dark areas better visible color_func = sqrt e = 0.75/12. # subreflector diameter/ primary diameter tau = 0.25 # 12dB taper for y in range(h): if (y==h/2): uniformProfile = [] parabolicProfile = [] xprofile = [] for x in range(w): # make xx and yy go from (-0.5 to +0.5) * scalex xx = ((x + .5) / w - .5) * scalex yy = ((y + .5) / h - .5) * scaley v = hypot(xx, yy) xx = ((x + .5) / w - .5) * w yy = ((y + .5) / h - .5) * h vg = hypot(xx, yy) v1 = v2 = 0.5 if v != 0.: v1 = (2*scipy.special.j1(v) / v \ - 2*(e**2)*scipy.special.j1(v*e)/(v*e)) ** 2 v1 /= (1-e**2)**2 v2 = (2*tau*scipy.special.j1(v) / v \ - 2*(e**2)*scipy.special.j1(v*e)/(v*e) \ + 4*(1-tau)*scipy.special.jn(2,v)/(v**2) ) ** 2 v2 /= (1-e**2)**2 imageUniform[y, x] = color_func(v1) imageParabolic[y,x] = color_func(v2) if (y==h/2): uniformProfile.append(v1) parabolicProfile.append(v2) xprofile.append(x) max_val = imageUniform.max() uniformProfile /= np.max(uniformProfile) parabolicProfile /= np.max(parabolicProfile) for ratio in [0.1, 0.25, 0.5]: uniformProfileMax = np.max(uniformProfile) parabolicProfileMax = np.max(parabolicProfile) for i in xprofile: if (uniformProfile[i]>ratio*uniformProfileMax and uniformProfile[i-1]ratio*uniformProfileMax): x1fwhm = i if (parabolicProfile[i]>ratio*parabolicProfileMax and parabolicProfile[i-1]ratio*parabolicProfileMax): p1fwhm = i uniformFWHM = ((x1fwhm-x0fwhm)*fov)/h parabolicFWHM = ((p1fwhm-p0fwhm)*fov)/h print "%.2f: Uniform = %.2f pixels = %.2f arcsec" % (ratio,x1fwhm-x0fwhm, uniformFWHM) print "%.2f: Parabolic = %.2f pixels = %.2f arcsec" % (ratio, p1fwhm-p0fwhm, parabolicFWHM) if (ratio == 0.5): fwhmratio = parabolicFWHM/uniformFWHM print "%.2f: Ratio = %.3f" % (ratio, parabolicFWHM/uniformFWHM) pb.clf() desc = pb.subplot(111) x = xprofile - np.ones(len(xprofile))*w/2 x *= plateScale pb.plot(x,uniformProfile,'-b',x,parabolicProfile,'-r') pb.title('Theoretical beam profile of ALMA 12m antenna with 0.75m blockage') x0 = np.min(x) + 0.05*(np.max(x)-np.min(x)) pb.text(x0,1.13,'blue = Uniform illumination, FWHM = 17.6 arcsec') pb.text(x0,1.07,'red = Parabolic illumination (-12dB edge taper), FWHM=19.8 arcsec') pb.text(x0,1.01,'ratio = %.3f' % (fwhmratio)) pb.ylim([0,1.19]) #pb.xlabel('Pixels (= 0.1 arcsecs)') pb.xlabel('Arcsec') pb.ylabel('Intensity') pb.draw() pb.savefig('profiles.png') image_file = Image.new('L', (w, h)) uniformData = numpy.zeros((w,h), dtype=numpy.float64) parabolicData = numpy.zeros((w,h), dtype=numpy.float64) for y in range(h): for x in range(w): c = int(2**16 * imageUniform[y, x] / max_val) image_file.putpixel((x, y), c) uniformData[x,y] = c c = int(2**16 * imageParabolic[y, x] / max_val) parabolicData[x,y] = c hdu = pyfits.PrimaryHDU(data=uniformData, header=pyfits.Header()) os.system('rm -f uniform.fits') hdu.writeto('uniform.fits') hdu = pyfits.PrimaryHDU(data=parabolicData, header=pyfits.Header()) os.system('rm -f tapered.fits') hdu.writeto('tapered.fits') # imview(raster='airydisk.fits', contour={'file':'airydisk.fits','levels':[0.5]}, zoom = 2) # imview(raster='airydisk.fits', contour={'file':'airydisk.fits','levels':[0.5]}, zoom = 1)