diff --git a/README.md b/README.md
index 76f1bfc..9de6e7f 100644
--- a/README.md
+++ b/README.md
@@ -50,7 +50,7 @@ Curved-sky maps
 * Primary modules: `curvedsky`, `lensing`
 * Internal dependencies: `sharp` + [flat-sky maps]
 * External dependencies: [`libsharp`](http://sourceforge.net/projects/libsharp/), [`cython`](http://cython.org)
-* To build: `make sharp`
+* To build: `make sharp`, `make interpol` 
 
 Mapmaking
 ---------
diff --git a/lensing.py b/lensing.py
index 38fe5ec..046cf9b 100644
--- a/lensing.py
+++ b/lensing.py
@@ -74,21 +74,31 @@ def lens_map_flat(cmb_map, phi_map):
 
 ######## Curved sky lensing ########
 
-def rand_map(shape, wcs, ps_lensinput, lmax=None, maplmax=None, dtype=np.float64, seed=None, oversample=2.0, spin=2, output="l", geodesic=True, verbose=False, delta_theta=None):
+def rand_map(shape, wcs, ps_lensinput, lmax=None, maplmax=None, dtype=np.float64, seed=None, oversample=2.0, spin=2, output="l", geodesic=True, verbose=False, delta_theta=None, phi_seed=None, separate_phi_from_cmb=False):
 	import curvedsky, sharp
 	ctype   = np.result_type(dtype,0j)
 	# Restrict to target number of components
 	oshape  = shape[-3:]
 	if len(oshape) == 2: shape = (1,)+tuple(shape)
+
+	#van Engelen added this:
+
 	# First draw a random lensing field, and use it to compute the undeflected positions
-	if verbose: print("Generating alms")
-	alm, ainfo = curvedsky.rand_alm(ps_lensinput, lmax=lmax, seed=seed, dtype=ctype, return_ainfo=True)
-	phi_alm, cmb_alm = alm[0], alm[1:1+shape[-3]]
+	if not separate_phi_from_cmb:
+
+		#AVE - this was the default option.
+		if verbose: print("Generating alms")
+		alm, ainfo = curvedsky.rand_alm(ps_lensinput, lmax=lmax, seed=seed, dtype=ctype, return_ainfo=True)
+		phi_alm, cmb_alm = alm[0], alm[1:1+shape[-3]]
+		del alm
+        else:
+		if verbose: print("Generating alms, separating phi from cmb")
+		phi_alm, phi_ainfo = curvedsky.rand_alm(ps_lensinput[0, 0, :], lmax=lmax, seed=phi_seed, dtype=ctype, return_ainfo=True)
+		cmb_alm, cmb_ainfo = curvedsky.rand_alm(ps_lensinput[1:, 1:, :], lmax=lmax, seed=seed, dtype=ctype, return_ainfo=True)
 	# Truncate alm if we want a smoother map. In taylens, it was necessary to truncate
 	# to a lower lmax for the map than for phi, to avoid aliasing. The appropriate lmax
 	# for the cmb was the one that fits the resolution. FIXME: Can't slice alm this way.
 	#if maplmax: cmb_alm = cmb_alm[:,:maplmax]
-	del alm
 	if delta_theta is None: bsize = shape[-2]
 	else:
 		bsize = utils.nint(abs(delta_theta/utils.degree/wcs.wcs.cdelt[1]))