diff --git a/zdm/beam_generator/DSA110_beamfile.dat b/zdm/beam_generator/DSA110_beamfile.dat
new file mode 100644
index 0000000..3de5b32
--- /dev/null
+++ b/zdm/beam_generator/DSA110_beamfile.dat
@@ -0,0 +1,9 @@
+1.405e9        #MHz obs frequency
+3.4         # FWHM [deg]
+34.37 0     # pointing angle (zenith, azimuth) of primary beam. Centre is 71.6 deg - 37:14:02 = 34.37
+48          # number of antennas
+-200 0      # ant 1
+200 0       # ant 48
+256         # number of beams
+-3.764     # azimuth angle of beam0 [deg]. 2.125 / sin(34.37)
+3.764      # azimuth angle of final beam [deg] Adding 255 / 60. to above
diff --git a/zdm/beam_generator/FAST.dat b/zdm/beam_generator/FAST.dat
new file mode 100644
index 0000000..b29779b
--- /dev/null
+++ b/zdm/beam_generator/FAST.dat
@@ -0,0 +1,29 @@
+2 -1 1 1.04 60.
+# above are: beam method (2: read tabe below)
+# -1: individual half power beam widths
+# 1: beam widths are HPBW
+# 1.04: scale beam widths by this value (1.3 GHz to 1.25 GHz)
+# 60. units of below, here arcminutes
+# data from https://iopscience.iop.org/article/10.1088/1674-4527/20/5/64/pdf
+# bX    bY      HPBW1300    A1250
+-0.04   -0.02   3.01        0.622
+5.76    -0.01   3.04        0.581
+2.86    -4.98   3.03        0.617
+-2.89   -5.01   3.05        0.598
+-5.78   -0.02   3.03        0.639
+-2.92   4.98    3.07        0.604
+2.85    4.97    3.02        0.602
+11.55   -0.02   3.04        0.554
+8.65    -5.01   3.04        0.577
+5.78    -10.02  3.08        0.586
+0.00    -10.04  3.05        0.580
+-5.78   -10.07  3.07        0.576
+-8.67   -5.05   3.04        0.604
+-11.61  -0.02   3.10        0.573
+-8.68   4.99    3.06        0.568
+-5.83   10.02   3.10        0.529
+-0.03   10.00   3.02        0.565
+5.76    10.00   3.05        0.557
+8.63    4.98    3.03        0.569
+
+
diff --git a/zdm/beam_generator/README.txt b/zdm/beam_generator/README.txt
new file mode 100644
index 0000000..7e5b85a
--- /dev/null
+++ b/zdm/beam_generator/README.txt
@@ -0,0 +1,24 @@
+This folder contains files used to generate simulate telescope beamshapes,
+and convert them into Beamdata. These programs are independent of the zdm code,
+but provide useful into to it.
+
+The relevant programs are:
+
+###### sim_fast_mutibeam.py ########
+This program simulates the FAST multibeam system using data
+contained in "FAST.dat", which is taken from
+https://iopscience.iop.org/article/10.1088/1674-4527/20/5/64/pdf
+With little adaptation, it could also simulate Parkes and Arecibo
+multibeam systems.
+
+###### sim_askap_beams.py ########
+This program is used to simulate the envelope formed by
+ASKAP beams.
+
+
+###### sim_DSA_beams.py ######
+Used to simulate the beamhape from tied beams generated by
+an array of coherently added antennas. It is used to simulate
+e.g. DSA beamshapes - although the exact locations of DSA
+antennas are not known. It could easily be adapted to an arbitrary
+tied array (e.g. MeerKAT)
diff --git a/zdm/beam_generator/compare_fast_parkes.py b/zdm/beam_generator/compare_fast_parkes.py
new file mode 100644
index 0000000..1751892
--- /dev/null
+++ b/zdm/beam_generator/compare_fast_parkes.py
@@ -0,0 +1,52 @@
+"""
+Compares FAST and PARKES multibeam profiles
+
+Parkes is scaled to FAST according to simple scaling
+relations.
+
+Note that the Parkes profile is originally taken
+from a full beam simulation.
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+
+def main():
+    
+    
+    
+    FASTb = np.load("FAST_log_bins.npy")
+    FASTh = np.load("FAST_log_hist.npy")
+    
+    Parkesb = np.load("../data/BeamData/parkes_mb_log_bins.npy")
+    Parkesh = np.load("../data/BeamData/parkes_mb_log_hist.npy")
+    
+    lPb = np.log10(Parkesb)
+    lPb = lPb[:-1]
+    lPb += (lPb[1]-lPb[0])/2.
+    Parkesb = 10**lPb
+    
+    lFb = np.log10(FASTb)
+    lFb = lFb[:-1]
+    lFb += (lFb[1]-lFb[0])/2.
+    FASTb = 10**lFb
+    
+    total = np.sum(FASTh)
+    print("Total FAST is ",total)
+    # scaling of Parkes hist
+    Parkesh *= 19/13 * (64/300)**2.
+    
+    plt.figure()
+    
+    plt.plot(Parkesb,Parkesh,label="Scaled Parkes")
+    plt.plot(FASTb,FASTh,label="Gaussian FAST sim")
+    plt.xlabel("B")
+    plt.ylabel("$\\Omega(B)$")
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig("parkes_fast_comparison.png")
+    plt.close()
+    
+main()
+    
diff --git a/zdm/beam_generator/plot_askap_avals.py b/zdm/beam_generator/plot_askap_avals.py
new file mode 100644
index 0000000..78e5610
--- /dev/null
+++ b/zdm/beam_generator/plot_askap_avals.py
@@ -0,0 +1,44 @@
+"""
+Simple script to plot the "avals": values of
+beam sensitivity for ASKAP as a function of
+distance from boresight.
+
+These "avals" are taken from James et al.
+
+"""
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+
+#### Defines theta offsets from boresights [deg] #####
+thetas = np.linspace(0,4,51)
+
+##### Loads in avals from Hotan et al #####
+data = np.loadtxt("avals.dat")
+x = data[:,0]
+y = data[:,1]
+h_avals = np.interp(thetas,x,y)
+
+##### Calculates avals from James et al. #####
+theta_off = 0.8
+sigma_theta = 3.47
+
+small = np.where(thetas < theta_off)
+
+ftheta = np.exp(-0.5 * ((thetas-theta_off)/sigma_theta)**2.)
+ftheta[small] = 1.
+
+##### plots the values ######
+plt.figure()
+
+plt.plot(thetas,ftheta,label="James et al")
+plt.plot(thetas,h_avals,label="Hotan et al")
+plt.xlim(0,4.)
+plt.ylim(0.63,1.06)
+plt.xlabel("Offset from boresight [deg]")
+plt.ylabel("Relative beam amplitude [B/B(0)]")
+plt.legend()
+plt.tight_layout()
+plt.savefig('avals.pdf')
+plt.close()
diff --git a/zdm/beam_generator/sim_DSA_beams.py b/zdm/beam_generator/sim_DSA_beams.py
new file mode 100644
index 0000000..7cf0e45
--- /dev/null
+++ b/zdm/beam_generator/sim_DSA_beams.py
@@ -0,0 +1,394 @@
+"""
+This program simulates a single Gaussian beam, and
+assumes N formed beams from M identical antennas are formed.
+
+Values relevant to DSA are hard-coded.
+"""
+
+import numpy as np
+
+def main():
+    """
+    main program. Apologies that inputs and outputs are hard-coded for now!
+    """
+    
+    Freq,FWHM,Nants,axs,ays,Nbeams,bzs,bas,pz,pa = read_simple_beamfile('DSA110_beamfile.dat')
+    
+    Nbeams = 1
+    bzs = bzs[0:1]
+    bas = bas[0:1]
+    
+    gridz,grida,weights = create_grid(pz,pa,300,6.)
+    
+    envelope = sim_formed_beams(axs,ays,Freq,bzs,bas,gridz,grida)
+    
+    primary = apply_primary_beamshape(pz,pa,FWHM,gridz,grida)
+    
+    combined = primary * envelope
+    
+    combined /= np.max(combined)
+    
+    bins = np.logspace(-4,0,401)
+    h,b=np.histogram(combined.flatten(),bins,weights=weights.flatten())
+    
+    np.save('DSA_log_hist.npy',h)
+    np.save('DSA_log_bins.npy',b)
+    
+    plot=True
+    if plot:
+        
+        from matplotlib import pyplot as plt
+        # note that b will have one more element
+        bcs = b[:-1] + (b[1]-b[0])/2.
+        plt.figure()
+        plt.plot(bcs,h)
+        plt.xlabel('B')
+        plt.ylabel('$\\Omega(B)$')
+        plt.xscale('log')
+        plt.tight_layout()
+        plt.savefig('dsa_omega_b.pdf')
+        plt.close()
+        
+        plt.figure()
+        plt.imshow(envelope,extent=(grida[0,0],grida[0,-1],gridz[0,0],gridz[-1,0]),origin='lower')
+        plt.xlabel('azimuth angle [deg]')
+        plt.ylabel('zenith [deg]')
+        plt.tight_layout()
+        plt.savefig('dsa_tied_beam_pattern.pdf')
+        plt.close()
+    
+        
+        plt.figure()
+        plt.imshow(primary,extent=(grida[0,0],grida[0,-1],gridz[0,0],gridz[-1,0]),origin='lower')
+        plt.xlabel('azimuth angle [deg]')
+        plt.ylabel('zenith [deg]')
+        plt.tight_layout()
+        plt.savefig('dsa_primary_beam_pattern.pdf')
+        plt.close()
+        
+        plt.figure()
+        plt.imshow(combined,extent=(grida[0,0],grida[0,-1],gridz[0,0],gridz[-1,0]),origin='lower')
+        plt.xlabel('azimuth angle [deg]')
+        plt.ylabel('zenith [deg]')
+        plt.tight_layout()
+        plt.savefig('dsa_beam_pattern.pdf')
+        plt.close()
+
+def apply_primary_beamshape(pz,pa,FWHM,gridz,grida):
+    """
+    Applies primary beamshape correction to formed beams
+    
+    pz [float]: zenith value of beam pointing direction [degrees]
+    
+    pa [float]: azimuth value of beam pointing direction [degrees]
+    
+    FWHM [float]: full width, half-max of the Gaussian beam in degrees
+    
+    gridz [array of floats]: zenith angle values of each grid point [degrees]
+    
+    grida [array of floats]: azimuth angle values of each grid point [degrees]
+    
+    Returns: primary beam values on the grid
+    """
+    
+    # calculates distances from grid points to primary
+    # beam centre
+    px,py,pz = get_xyz(pz,pa)
+    
+    # gets unit vectors in grid directions
+    gx,gy,gz = get_xyz(gridz,grida)
+    
+    cosines = px*gx + py*gy + pz*gz
+    angles = np.arccos(cosines)
+    
+    sigma = (np.pi/180. * FWHM/2.)*(2*np.log(2))**-0.5
+    
+    primary = Gauss(angles,sigma)
+    return primary
+
+def Gauss(r,sigma):
+    """
+    Simple Gaussian function, normalised to a peak
+    amplitude of unity
+    
+    Inputs:
+        r [float] is a radial offset from centre
+        sigma [float] is a std deviation in the same units
+    
+    Return:
+        returns value of Gaussian at r points
+    """
+    
+    return np.exp(-0.5 * (r/sigma)**2)
+
+  
+def get_xyz(z,a):
+    """
+    returns x,y,z coordinates of given zenith, azimuth position
+    
+    z [float]: zenith angle [degrees]
+    
+    a [flot]: azimuth angle [degrees]
+    
+    Returns: x,y,z unit vectors in direction of a,z
+    """
+    zr = z * np.pi/180.
+    ar = a * np.pi/180.
+    
+    cz = np.cos(zr)
+    sz = np.sin(zr)
+    ca = np.cos(ar)
+    sa = np.sin(ar)
+    
+    x = sz * sa
+    y = sz * ca
+    z = cz
+    return x,y,z
+
+def sim_formed_beams(axs,ays,nu,bz,ba,gz,ga):
+    """
+    This program simulated the formed beamshape
+    given a grid of antenna positions, a frequency,
+    and a pointing direction (zenith, azimuth)
+    
+    Inputs:
+        axs: numpy array of antenna x positions (East, m)
+        ays: numpy array of antenna y positions (North, m)
+        nu: central frequency (Hz)
+        zenith: pointing zenith angle (deg)
+        azimuth: pointing azimuth, East from North (deg)
+        bz: beam phase centre zenith (deg)
+        ba: beam phase centre azimuth (deg)
+        gz: grid of zenith angles (deg)
+        ga: grid of azimuth angles (deg)
+    """
+    
+    # get number of beams. Allows for there to be a single beam or multiple
+    bz = np.array(bz)
+    ba = np.array(ba)
+    Nbeams = bz.size
+    
+    # gets shape of sky grid, gz.
+    # NOTE: shape of 'gz' and 'ga' arrays does NOT need to be
+    # in simple deltaz, deltaa format - but it helps to label them as such
+    # could just be in x-y format
+    Nz,Na = gz.shape
+    
+    wavelength = 2.99792458e8/nu
+    
+    # apparent baselines in units of lambda
+    #xbaseline = axs * xcosine / wavelength
+    #ybaseline = ays * ycosine / wavelength
+    
+    # first: calculate direction cosines in x and y directions for
+    # each grid point
+    
+    gxcosines,gycosines = get_cosines(gz,ga)
+    
+    # these are the distances incurred for each m in the x and y directions respectively
+    # now multiply by x and y distances in units of phase factors
+    phase_factor = 2. * np.pi / wavelength
+    
+    #calculate delay corresponding to phase centre
+    bxcosines,bycosines = get_cosines(bz,ba)
+    
+    # set up grid for all beams
+    formed_beams = np.zeros([Nbeams,Nz,Na],dtype='complex')
+    
+    # loop over antennas, calculating delay for each
+    for iant,ax in enumerate(axs):
+        dax = ax * phase_factor
+        day = ays[iant] * phase_factor
+        
+        # gives delays to each position for this antenna
+        # units are phase
+        delays = dax * gxcosines + day * gycosines
+        
+        
+        # gives base delay for each beam
+        # units also phase
+        delay0s = dax * bxcosines + day * bycosines
+        
+        # loops over beam
+        for ibeam,delay0 in enumerate(delay0s):
+            # delays to each position for this beam
+            bdelays = delays - delay0
+            formed_beams[ibeam,:,:] += np.exp(bdelays * 1j)
+            
+    
+    mags = np.abs(formed_beams)
+    envelope = np.max(mags,axis=0)
+    
+    return envelope
+    
+    
+ 
+def create_grid(zen,az,Ngrid,extent):
+    """
+    creates a grid of sky positions,
+    along with steradians of each cell
+    
+    The grid is centred on position zen,az
+    and is equally spaced in zenith and azimuth coordinates
+    
+    extent is the total degree extent in each direction about the centre
+    
+    Returns the grid zenith, grid azimuth,
+    and solid angle values of the grid
+    """ 
+    dgrid = extent / (Ngrid-1.)
+    
+    zmin = zen - (Ngrid-1.)/2.*dgrid
+    zmax = zmin + (Ngrid-1.)*dgrid
+    zvals = np.linspace(zmin,zmax,Ngrid)
+    
+    # could scale these by zenith angle perhaps?
+    amin = az - (Ngrid-1.)/2.*dgrid
+    amax = amin + (Ngrid-1.)*dgrid
+    avals = np.linspace(amin,amax,Ngrid)
+    
+    apoints = np.repeat(avals,Ngrid)
+    apoints = apoints.reshape([Ngrid,Ngrid])
+    apoints = apoints.T
+    
+    zpoints = np.repeat(zvals,Ngrid)
+    zpoints = zpoints.reshape([Ngrid,Ngrid])
+    
+    # zenith angle increment
+    dz = dgrid * np.pi/180.
+    # aziumuth angle increment
+    da = dgrid * np.pi/180. * np.sin(zvals*np.pi/180.)
+    wgt = da * dz
+    
+    weights = np.repeat(wgt,Ngrid)
+    weights = weights.reshape([Ngrid,Ngrid])
+    
+    return zpoints,apoints,weights
+    
+def read_simple_beamfile(inputfile):
+    """
+    Reads a file containing information on the telescope
+    """
+    with open(inputfile) as infile:
+        lines = infile.readlines()
+        
+        Freq=float(lines[0].split()[0])
+        FWHM=float(lines[1].split()[0])
+        pz=float(lines[2].split()[0])
+        pa=float(lines[2].split()[1])
+        
+        Nants = int(lines[3].split()[0])
+        ax0 = float(lines[4].split()[0])
+        ay0 = float(lines[4].split()[1])
+        ax1 = float(lines[5].split()[0])
+        ay1 = float(lines[5].split()[1])
+        
+        Nbeams = int(lines[6].split()[0])
+        ba0 = float(lines[7].split()[0])
+        ba1 = float(lines[8].split()[0])
+    
+    # initialises beams
+    bas = np.linspace(ba0,ba1,Nbeams)
+    bzs = np.full([Nbeams],pz)
+    
+    # initialises antenna locations
+    axs = np.linspace(ax0,ax1,Nants)
+    ays = np.linspace(ay0,ay1,Nants)
+    
+    return Freq,FWHM,Nants,axs,ays,Nbeams,bzs,bas,pz,pa
+
+def read_complex_beamfile(inputfile):
+    """
+    Reads a file containing information on the telescope
+    """
+    with open(inputfile) as infile:
+        lines = infile.readlines()
+        for i,line in enumerate(lines):
+            words = line.split()
+            if i==0:
+                # metadata
+                Nants = int(words[0])
+                Nbeams = int(words[1])
+                Freq = float(words[2])
+                
+                axs = np.zeros([Nants])
+                ays = np.zeros([Nants])
+                
+                bzs = np.zeros([Nbeams])
+                bas = np.zeros([Nbeams])
+                
+            elif i==1:
+                # prmary beam data
+                pz = float(words[0]) # primary beam zenith angle [deg]
+                pa = float(words[1]) # primary beam azimuth angle [deg]
+                pw = float(words[2]) # primary beam width [deg]
+                pm = int(words[3]) # method for specifying primary beam
+                
+                # converts to Gaussian sigma
+                psigma = process_beam_width(pw,pm)
+            elif i < Nants +2:
+                # this is antenna data
+                ax = float(words[0]) # primary beam zenith angle [deg]
+                ay = float(words[1])
+                axs[i-2] = ax
+                ays[i-2] = ay
+                
+            else:
+                # this is formed beam data
+                bz = float(words[0])
+                ba = float(words[1])
+                bzs[i-2-Nants] = bz
+                bas[i-2-Nants] = ba
+                
+def process_beam_width(value,method):
+    """
+    converts different measures of beam width
+    to a standard Gaussian sigma
+    
+    Input:
+        value [float]: value of beamwidth
+        method [int]: definition of value
+    
+    Return:
+        sigma [float]: std dev in degrees of Gaussian beam     
+    """
+    
+    
+    if method == 0:
+        # actual sigma of beam
+        sigma = value
+    elif method == 1:
+        # sigma is FWHM (aka HPBW) of beam:
+        sigma = (value/2.)*(2*np.log(2))**-0.5
+    else:
+        raise ValueError("Invalid method for beam width method found")
+    
+    return sigma
+
+    
+def get_cosines(z,a):
+    """
+    Get direction cosines corresponding to the zentih angle z and azimuth angle a
+    """       
+    # grid zenith and azimuth angles in radians
+    zr = z * np.pi/180.
+    ar = a * np.pi/180.
+    
+    
+    
+    # cosines and sines of those angles
+    cz = np.cos(zr)
+    ca = np.cos(ar)
+    sa = np.sin(ar)
+    
+    
+    
+    # direction cosines: azimuth is East from North,
+    # zenith is from zenith, x is East and y is North
+    xcosines = cz * sa
+    ycosines = cz * ca
+    
+    return xcosines,ycosines
+
+
+main()
diff --git a/zdm/beam_generator/sim_askap_beams.py b/zdm/beam_generator/sim_askap_beams.py
new file mode 100644
index 0000000..4682ed1
--- /dev/null
+++ b/zdm/beam_generator/sim_askap_beams.py
@@ -0,0 +1,378 @@
+"""
+This program simulates 36 Gaussian beams on the sky,
+nominally from ASKAP, and produces their beam pattern.
+
+It assumes that all beams are perfect Gaussians, 
+and includes a tapering of beam performance away from
+boresight according to Hotan et al (an alternative is
+provided in "askap_beam_amps", but not used)
+
+It writes out beamfile information for use by the zdm code
+"""
+
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+
+def main(plot=True,pattern="square",pitch=0.9,freq=1400,gsize=800,spacing=0.02):
+    """
+    
+    Pattern:
+        Square: square array
+        Closepack: triangular grid
+    
+    Sep:
+        Beam separation in degrees
+    
+    Frequency:
+        frequency in MHz
+        
+    gsize [int]: number of grid elements on which to calculate beam pattern
+    
+    spacing [float]: spacing in degrees between calculation points. Total
+        area is (gsize * spacing)^2
+        
+    All files produced have pattern_sep_freq notation
+    """
+    
+    
+    # defines grid in x,y - this will be waaaay too small!
+    gsize = gsize
+    xgrid,ygrid = make_grid(gsize,spacing)
+    
+    ###### defines antenna pattern ########
+    beamx, beamy, beamw, beama = gen_askap_beams(pattern,pitch,freq*1e6)
+    
+    Nbeams = beamx.size
+    
+    # loops through beams, filling the grid
+    beams = np.zeros([Nbeams,gsize,gsize])
+    for i in np.arange(Nbeams):
+        beams[i,:,:] = beama[i] * grid_beam(xgrid,ygrid,beamx[i],beamy[i],beamw[i])
+    
+    envelope = np.max(beams,axis=0)
+    
+    if plot:
+        plot_envelope(envelope,xgrid[0,:])
+    
+    
+    # makes a beam histogram
+    bins = np.logspace(-4,0,401)
+    
+    weights = calculate_weights(xgrid,ygrid)
+    
+    # calculates single value
+    contributions = weights * envelope**1.5
+    total = np.sum(contributions) * 100
+    
+    tsys = get_tsys(freq)
+    
+    adjusted = total * (80./tsys)**1.5
+    
+    print(pattern,freq/1e6,pitch,total,adjusted)
+    
+    
+    ########### Makes a histogram of the beamshape ###########
+    tsys0 = get_tsys(1272.5) # normalising frequency
+    envelope *= tsys0/tsys # adjusts so B=1 is same sensitivity at 1272.5 GHz
+    # makes a beam histogram
+    bins = np.logspace(-4,0,401)
+    h,b=np.histogram(envelope.flatten(),bins,weights=weights.flatten())
+    
+    np.save('ASKAP/'+pattern+'_'+str(pitch)+'_'+str(freq)+'.npy',h)
+    np.save('ASKAP/log_bins.npy',b)
+    
+    if False:
+        db = bins[1]/bins[0]
+        bcs = b[:-1]*db**0.5
+        plt.figure()
+        plt.plot(bcs,h)
+        lat50b = np.load('../data/BeamData/lat50_log_bins.npy')
+        lat50h = np.load('../data/BeamData/lat50_log_hist.npy')
+        plt.plot(10**lat50b,lat50h*1000/400*35/40)
+        
+        plt.xscale('log')
+        plt.yscale('log')
+        plt.show()
+    
+    #h,b=np.histogram(envelope.flatten(),bins,weights=weights.flatten())
+    
+    #if plot:
+    #    plot_histogram(h,b)
+
+def get_tsys(freq,plot=False):
+    """
+    reads in interpolated data from Hotan et al
+    
+    Freq: frequency in MHz
+    """   
+    
+    data = np.loadtxt("askap_tsys.dat")
+    x = data[:,0]
+    y = data[:,1]
+    deg = 8
+    coeffs = np.polyfit(x,y,deg)
+    
+    
+    if plot:
+        xvals = np.linspace(700,1800,201)
+        yvals = np.polyval(coeffs,xvals)
+        plt.figure()
+        plt.plot(xvals,yvals,color='grey')
+        plt.plot(x,y,linestyle="",color='red',marker='o')
+        plt.ylim(0,130)
+        plt.savefig("askap_tsys.pdf")
+        plt.close()
+    
+    Tsys = np.polyval(coeffs,freq)
+    return Tsys
+    
+def gen_askap_beams(pattern,sep,freq):
+    """
+    Generates beam centres
+    
+    pattern [string]: square or closepack
+    
+    sep [float, deg]: beam separation
+    
+    freq [float]: frequency in Hz
+    """
+    
+    if pattern == "square":
+        # square grid
+        x=np.linspace(0,5.*sep,6)
+        x=np.repeat(x,6)
+        x = x.reshape([6,6])
+        y=np.copy(x)
+        y=y.T
+        x -= sep*2.5
+        y -= sep*2.5
+    elif pattern == "closepack":
+        # triangular grid
+        x = np.zeros([6,6])
+        y = np.zeros([6,6])
+        for iy in np.arange(6):
+            if iy % 2 == 0:
+                x0 = 0.
+            else:
+                x0 = 0.5 * sep
+            x[iy,:] = np.arange(6) * sep + x0
+            y[iy,:] = iy*sep*3**0.5/2.
+        x -= np.sum(x)/36.
+        y -= np.sum(y)/36.
+    
+    x = x.flatten()
+    y = y.flatten()
+    
+    # implements function of distance
+    dists = (x**2 + y**2)**0.5
+    
+    # from Hotan et al
+    avals = get_avals(dists)
+    
+    
+    # function of frequency only
+    # Gaussian sigma
+    diam = 12.
+    width = get_beamwidth(freq)
+    beamw = np.full([36],width)
+    
+    return x, y, beamw, avals
+        
+def get_beamwidth(freq):
+    """
+    Gets beamwidth with frequency in Hz
+    returns Gaussian sigma in degrees
+    
+    freq [float]: frequency in HHz
+    """
+    D = 12.
+    c_light = 299792458.
+    HPBW=1.09*(c_light/(freq))/D # from RACS
+    sigma=(HPBW/2.)*(2*np.log(2))**-0.5
+    deg_sigma = sigma * 180./np.pi
+    return deg_sigma
+
+def askap_beam_amps(theta):
+    """
+    Equation 7 from lat50 paper of James et al.
+    
+    """
+    
+    theta_off = 0.8
+    sigma_theta = 3.47
+    
+    # gets locations where amplitude is unity
+    shape = theta.shape
+    theta = theta.flatten()
+    small = np.where(theta < theta_off)
+    
+    # sets entire array to Gaussian
+    amps = np.exp(-0.5 * ((theta-theta_off)/sigma_theta)**2)
+    # resets small values to unity
+    amps[small] = 1.
+    
+    amps = amps.reshape(shape)
+    return amps
+
+def calculate_weights(xvals,yvals):
+    """
+    calculates histogram weights, accounting for sine
+    """
+    # converts to steradians
+    xvals *= np.pi/180.
+    yvals *= np.pi/180.
+    
+    # calculates distances from the origin
+    dists = (xvals**2 + yvals**2)**0.5
+    
+    cosfactor = np.cos(dists) #deformation factor due to curvature of sky
+    base = (xvals[0,1]-xvals[0,0])**2 # pixel size
+    weights = cosfactor * base
+    return weights
+
+def plot_histogram(h,b):
+    """
+    makes a simple plot of Omega(b)
+    """
+    
+    from matplotlib import pyplot as plt
+    # note that b will have one more element
+    bcs = b[:-1] + (b[1]-b[0])/2.
+    plt.figure()
+    plt.plot(bcs,h)
+    plt.xlabel('B')
+    plt.ylabel('$\\Omega(B)$')
+    plt.xscale('log')
+    plt.tight_layout()
+    plt.savefig('omega_b.pdf')
+    plt.close()
+    
+
+def plot_envelope(env,coords):
+    """
+    Makes simple plot of beamshape
+    """
+    
+    from matplotlib import pyplot as plt
+    dc = coords[1]-coords[0]
+    themin = coords[0]-dc/2.
+    themax = coords[-1] + dc/2.
+    
+    plt.figure()
+    plt.imshow(env,extent=(themin,themax,themin,themax),origin='lower')
+    plt.xlabel('x [deg]')
+    plt.ylabel('y [deg]')
+    plt.tight_layout()
+    plt.savefig('askap_beam_pattern.pdf')
+    plt.close()
+   
+def make_grid(Npoints,ddeg):
+    """
+    Generates a grid in which to store beam values.
+    
+    Inputs:
+        Npoints [int]: number of points in each dimension.
+            Final array will be Npoints X Npoints
+        ddeg [float]: Degree increment between points
+    """
+    # generates xpoints 
+    extent = (Npoints-1.)*ddeg/2.
+    degrees = np.linspace(extent,-extent,Npoints)
+    ypoints = np.repeat(degrees,Npoints)
+    ypoints = ypoints.reshape([Npoints,Npoints])
+    
+    degrees = np.linspace(-extent,extent,Npoints)
+    xpoints = np.repeat(degrees,Npoints)
+    xpoints = xpoints.reshape([Npoints,Npoints])
+    xpoints = xpoints.T
+    
+    return xpoints,ypoints
+
+def process_beam_width(value,method):
+    """
+    converts different measures of beam width
+    to a standard Gaussian sigma
+    
+    Input:
+        value [float]: value of beamwidth
+        method [int]: definition of value
+    
+    Return:
+        sigma [float]: std dev in degrees of Gaussian beam     
+    """
+    
+    
+    if method == 0:
+        # actual sigma of beam
+        sigma = value
+    elif method == 1:
+        # sigma is FWHM (aka HPBW) of beam:
+        sigma = (value/2.)*(2*np.log(2))**-0.5
+    else:
+        raise ValueError("Invalid method for beam width method found")
+    
+    return sigma
+
+
+def get_avals(thetas):
+    """
+    Values extrated from Hotan et al.
+    """
+    
+    data = np.loadtxt("avals.dat")
+    x = data[:,0]
+    y = data[:,1]
+    avals = np.interp(thetas,x,y)
+    return avals
+    
+def grid_beam(xvals,yvals,beamx,beamy,beamw):
+    """
+    Fills grid with sensitivities of a beam
+    
+    Inputs:
+        xvals [np 2d array]: x coordinates (degrees)
+            of grid centres (not bin edges)
+        yvals [np 2d array]: y coordinates (degrees)
+            of grid centres (not bin edges)
+        beamx: coordinates of beam in x-direction
+        beamy: coordinates of beam in y-direction
+        beamw: beamwidth in degrees
+    
+    Returns:
+        grid (np 2d array), filled with beam values
+    """
+    
+    
+    offsets = ((xvals - beamx)**2 + (yvals - beamy)**2 )**0.5
+    values = Gauss(offsets,beamw)
+    return values
+
+
+def Gauss(r,sigma):
+    """
+    Simple Gaussian function, normalised to a peak
+    amplitude of unity
+    
+    Inputs:
+        r [float] is a radial offset from centre
+        sigma [float] is a std deviation in the same units
+    
+    Return:
+        returns value of Gaussian at r points
+    """
+    
+    return np.exp(-0.5 * (r/sigma)**2)
+
+# original
+# this is a hard-coded list of *every* combination of frequency and pitch angle
+# used by ASKAP to date
+flist = [ 819.,  832.,  864.,  920.,  936., 950.,  970., 980., 990., 1032., 1100.,
+            1112., 1266., 1272., 1297., 1320., 1407., 1418., 1632., 1641., 1656.]
+plist = [0.72, 0.75, 0.84, 0.9,  1.05, 1.1 ]
+array = ["closepack","square"]
+for freq in flist:
+    for p in plist:
+        for config in array:
+            main(freq=freq,pattern = config, pitch=p)
diff --git a/zdm/beam_generator/sim_fast_multibeam.py b/zdm/beam_generator/sim_fast_multibeam.py
new file mode 100644
index 0000000..af9ef5e
--- /dev/null
+++ b/zdm/beam_generator/sim_fast_multibeam.py
@@ -0,0 +1,284 @@
+"""
+This program simulates N Gaussian beams on the sky,
+takes their envelope, then converts them into
+beamdata files for use with surveys in the zdm code
+
+It is best utilised for multibeam beam patterns.
+
+E.g. use 'FAST.dat' for input. In theory,
+this could be easily adapted to other
+multibeam systems
+
+"""
+
+
+import numpy as np
+
+
+def main(plot=True):
+    
+    gsize = 700
+    xgrid,ygrid = make_grid(gsize,0.001)
+    
+    beamx, beamy, beamw, beama = read_beamfile('FAST.dat')
+    
+    Nbeams = beamx.size
+    
+    # loops through beams, filling the grid
+    beams = np.zeros([Nbeams,gsize,gsize])
+    for i in np.arange(Nbeams):
+        beams[i,:,:] = beama[i] * grid_beam(xgrid,ygrid,beamx[i],beamy[i],beamw[i])
+    
+    envelope = np.max(beams,axis=0)
+    
+    if plot:
+        plot_envelope(envelope,xgrid[0,:])
+    
+    
+    # makes a beam histogram
+    bins = np.logspace(-4,0,401)
+    
+    weights = calculate_weights(xgrid,ygrid)
+    h,b=np.histogram(envelope.flatten(),bins,weights=weights.flatten())
+    
+    
+    if plot:
+        plot_histogram(h,b)
+    
+    np.save('FAST_log_hist.npy',h)
+    np.save('FAST_log_bins.npy',b)
+
+def calculate_weights(xvals,yvals):
+    """
+    calculates histogram weights, accounting for sine
+    """
+    # converts to steradians
+    xvals *= np.pi/180.
+    yvals *= np.pi/180.
+    
+    # calculates distances from the origin
+    dists = (xvals**2 + yvals**2)**0.5
+    
+    cosfactor = np.cos(dists) #deformation factor due to curvature of sky
+    base = (xvals[0,1]-xvals[0,0])**2 # pixel size
+    weights = cosfactor * base
+    return weights
+
+def plot_histogram(h,b):
+    """
+    makes a simple plot of Omega(b)
+    """
+    
+    from matplotlib import pyplot as plt
+    # note that b will have one more element
+    bcs = b[:-1] + (b[1]-b[0])/2.
+    plt.figure()
+    plt.plot(bcs,h)
+    plt.xlabel('B')
+    plt.ylabel('$\\Omega(B)$')
+    plt.xscale('log')
+    plt.tight_layout()
+    plt.savefig('omega_b.pdf')
+    plt.close()
+    
+
+def plot_envelope(env,coords):
+    """
+    Makes simple plot of beamshape
+    """
+    
+    from matplotlib import pyplot as plt
+    dc = coords[1]-coords[0]
+    themin = coords[0]-dc/2.
+    themax = coords[-1] + dc/2.
+    
+    plt.figure()
+    plt.imshow(env,extent=(themin,themax,themin,themax),origin='lower')
+    plt.xlabel('x [deg]')
+    plt.ylabel('y [deg]')
+    plt.tight_layout()
+    plt.savefig('beam_pattern.pdf')
+    plt.close()
+   
+def make_grid(Npoints,ddeg):
+    """
+    Generates a grid in which to store beam values.
+    
+    Inputs:
+        Npoints [int]: number of points in each dimension.
+            Final array will be Npoints X Npoints
+        ddeg [float]: Degree increment between points
+    """
+    # generates xpoints 
+    extent = (Npoints-1.)*ddeg/2.
+    degrees = np.linspace(extent,-extent,Npoints)
+    ypoints = np.repeat(degrees,Npoints)
+    ypoints = ypoints.reshape([Npoints,Npoints])
+    
+    degrees = np.linspace(-extent,extent,Npoints)
+    xpoints = np.repeat(degrees,Npoints)
+    xpoints = xpoints.reshape([Npoints,Npoints])
+    xpoints = xpoints.T
+    
+    return xpoints,ypoints
+    
+def read_beamfile(beamfile):
+    """
+    Reads an input beam file.
+    The file must be in the correct format.
+    """
+    beamx = np.empty([0])
+    beamy = np.empty([0])
+    beamw = np.empty([0])
+    beama = np.empty([0])
+    
+    with open(beamfile) as infile:
+        
+        for i,line in enumerate(infile):
+            words = line.split()
+            
+            if len(line) <= 1:
+                break
+            
+            if i==0:
+                method = int(words[0]) # method of describing beams
+                if method != 1 and method != 2:
+                    print("invalid method")
+                width = float(words[1])
+                # If width == -1: read in from file
+                
+                # FWHM or sigma
+                wmethod = int(words[2])
+                
+                wscale = float(words[3])
+                # scale widths by value according to frequency
+                
+                unit = float(words[4]) # units - divide by this value
+                
+            else:
+                if line[0] == '#':
+                    continue
+                if method == 1:
+                    bx,by,bw,ba = gen_circular_beams(words,width)
+                elif method == 2:
+                    bx,by,bw,ba = read_beam_coordinates(words,width)
+                
+                beamx = np.append(beamx,bx)
+                beamy = np.append(beamy,by)
+                beamw = np.append(beamw,bw)
+                beama = np.append(beama,ba)
+    
+    beamx = beamx / unit
+    beamy = beamy / unit
+    beamw = beamw * wscale / unit
+    beamw = process_beam_width(beamw,wmethod)
+    beama = beama / np.max(beama)
+    
+    return beamx, beamy, beamw, beama
+    
+def process_beam_width(value,method):
+    """
+    converts different measures of beam width
+    to a standard Gaussian sigma
+    
+    Input:
+        value [float]: value of beamwidth
+        method [int]: definition of value
+    
+    Return:
+        sigma [float]: std dev in degrees of Gaussian beam     
+    """
+    
+    
+    if method == 0:
+        # actual sigma of beam
+        sigma = value
+    elif method == 1:
+        # sigma is FWHM (aka HPBW) of beam:
+        sigma = (value/2.)*(2*np.log(2))**-0.5
+    else:
+        raise ValueError("Invalid method for beam width method found")
+    
+    return sigma
+    
+    
+    
+def gen_circular_beams(words):
+    """
+    Generates beam centres for a circle of N
+    beams a distance r from the overall centre
+    and offset by some fractional spacing
+    
+    Inputs:
+        words (line describing inputs)
+    """
+    Nbeams = int(words[0])
+    r = float(words[1])
+    offset = float(words[2])
+    angles = (np.arange(Nbeams) + offset)* 2. * np.pi/ Nbeams
+    
+    beamx = r * np.cos(angles)
+    beamy = r * np.sin(angles)
+    return beamx,beamy
+
+def read_beam_coordinates(words,width):
+    """
+    Just reads in the beam centre coordinates
+    """
+    
+    beamx = float(words[0])
+    beamy = float(words[1])
+    
+    if width == -1:
+        if len(words) < 3:
+            raise ValueError("Cannot find beam width")
+        else:
+            bw = float(words[2])
+    else:
+        bw = width
+    
+    if len(words) == 4:
+        amp = float(words[3])
+    else:
+        amp == 1.
+    
+    return beamx,beamy,bw,amp
+
+
+def grid_beam(xvals,yvals,beamx,beamy,beamw):
+    """
+    Fills grid with sensitivities of a beam
+    
+    Inputs:
+        xvals [np 2d array]: x coordinates (degrees)
+            of grid centres (not bin edges)
+        yvals [np 2d array]: y coordinates (degrees)
+            of grid centres (not bin edges)
+        beamx: coordinates of beam in x-direction
+        beamy: coordinates of beam in y-direction
+        beamw: beamwidth in degrees
+    
+    Returns:
+        grid (np 2d array), filled with beam values
+    """
+    
+    offsets = ((xvals - beamx)**2 + (yvals - beamy)**2 )**0.5
+    values = Gauss(offsets,beamw)
+    return values
+
+def Gauss(r,sigma):
+    """
+    Simple Gaussian function, normalised to a peak
+    amplitude of unity
+    
+    Inputs:
+        r [float] is a radial offset from centre
+        sigma [float] is a std deviation in the same units
+    
+    Return:
+        returns value of Gaussian at r points
+    """
+    
+    return np.exp(-0.5 * (r/sigma)**2)
+    
+main()