diff --git a/src/madness/mra/funcimpl.h b/src/madness/mra/funcimpl.h
index 3e4a2acb6b2..07d93ab66bf 100644
--- a/src/madness/mra/funcimpl.h
+++ b/src/madness/mra/funcimpl.h
@@ -1709,12 +1709,6 @@ template<size_t NDIM>
         coeffT parent_to_child_NS(const keyT& child, const keyT& parent,
                                   const coeffT& coeff) const;
 
-        /// Get the scaling function coeffs at level n starting from NS form
-        // N=2^n, M=N/q, q must be power of 2
-        // q=0 return coeffs [N,k] for direct sum
-        // q>0 return coeffs [k,q,M] for fft sum
-        tensorT coeffs_for_jun(Level n, long q=0);
-
         /// Return the values when given the coeffs in scaling function basis
         /// @param[in] key the key of the function node (box)
         /// @param[in] coeff the tensor of scaling function coefficients for function node (box)
diff --git a/src/madness/mra/mra.h b/src/madness/mra/mra.h
index e14be5fb59d..84ca9731170 100644
--- a/src/madness/mra/mra.h
+++ b/src/madness/mra/mra.h
@@ -841,15 +841,6 @@ namespace madness {
         }
 
 
-        /// Get the scaling function coeffs at level n starting from NS form
-        Tensor<T> coeffs_for_jun(Level n, long mode=0) {
-            PROFILE_MEMBER_FUNC(Function);
-            make_nonstandard(true, true);
-            return impl->coeffs_for_jun(n,mode);
-            //return impl->coeffs_for_jun(n);
-        }
-
-
         /// Clears the function as if constructed uninitialized.  Optional fence.
 
         /// Any underlying data will not be freed until the next global fence.
diff --git a/src/madness/mra/mraimpl.h b/src/madness/mra/mraimpl.h
index e18e05b27d2..146314d1020 100644
--- a/src/madness/mra/mraimpl.h
+++ b/src/madness/mra/mraimpl.h
@@ -713,124 +713,6 @@ namespace madness {
         return result;
     }
 
-    /// Get the scaling function coeffs at level n starting from NS form
-    // N=2^n, M=N/q, q must be power of 2
-    // q=0 return coeffs [N,k] for direct sum
-    // q>0 return coeffs [k,q,M] for fft sum
-    template <typename T, std::size_t NDIM>
-    typename FunctionImpl<T,NDIM>::tensorT FunctionImpl<T,NDIM>::coeffs_for_jun(Level n, long q) {
-        MADNESS_ASSERT(is_nonstandard() && NDIM<=3);
-        tensorT r,r0;
-        long N=1<<n;
-        long M = (q ? N/q: N);
-        if (q==0) {
-            q = 1;
-            long dim[2*NDIM];
-            for (std::size_t d=0; d<NDIM; ++d) {
-                dim[d ] = N;
-                dim[d+NDIM] = cdata.k;
-            }
-            tensorT rr(2*NDIM,dim);
-            r0=r=rr;
-            //NNkk->MqMqkk, since fuse is not allowed. Now needs to move back to 2*NDIM, since tensor max dim is 6
-            //for (int d=NDIM-1; d>=0; --d) r.splitdim_inplace_base(d,M,q);
-        } else {
-            long dim[2*NDIM];
-            for (std::size_t d=0; d<NDIM; ++d) {
-                //dim[d+NDIM*2] = M;
-                dim[d+NDIM ] = N;
-                dim[d ] = cdata.k;
-            }
-            tensorT rr(2*NDIM,dim);
-            r0=rr;
-            /*vector<long> map(3*NDIM);
-              for (int d=0; d<NDIM; ++d) {
-              map[d]=d+2*NDIM;
-              map[NDIM+d]=2*d+1;
-              map[2*NDIM+d]=2*d;
-              }
-              r.mapdim_inplace_base(map);
-              //print(rr);
-              //for (int d=1; d<NDIM; ++d) rr.swapdim_inplace_base(2*NDIM+d,NDIM+d); //kkqqMM->kkqMqM
-              //print(rr);
-              //for (int d=0; d<NDIM; ++d) rr.swapdim_inplace_base(NDIM+2*d,NDIM+2*d-1); //kkqMqM->kkMqMq
-              //print(rr);
-              //for (int d=0; d<NDIM; ++d) rr.fusedim_inplace_base(NDIM+d); //->kkNN
-              //seems that this fuse is not allowed :(
-
-              //print(rr);
-              */
-            r=rr.cycledim(NDIM,0,-1); //->NNkk or MqMqkk
-        }
-        print("faking done M q r(fake) r0(real)",M,q,"\n", std::vector<long> (r.dims(),r.dims()+6),std::vector<long> (r0.dims(),r0.dims()+6));
-        ProcessID me = world.rank();
-        Vector<long,NDIM> t(N);
-
-        Vector<long,NDIM> powq, powN, powM;
-        long NDIM1 = NDIM-1;
-        powM[NDIM1]=powq[NDIM1]=powN[NDIM1]=1;
-        for (int d=NDIM1-1; d>=0; --d) {
-            powM[d] = powM[d+1]*M;
-            powq[d] = powq[d+1]*q;
-            powN[d] = powN[d+1]*N;
-        }
-        long powMNDIM = powM[0]*M;
-
-        for (IndexIterator it(t); it; ++it) {
-            keyT key(n, Vector<Translation,NDIM>(*it));
-            if (coeffs.owner(key) == me) {
-                typename dcT::iterator it = coeffs.find(key).get();
-                coeffT qq;
-
-                if (it == coeffs.end()) {
-                    // must get from above
-                    typedef std::pair< keyT,coeffT > pairT;
-                    Future<pairT> result;
-                    sock_it_to_me(key, result.remote_ref(world));
-                    const keyT& parent = result.get().first;
-                    //                        const tensorT& t = result.get().second.full_tensor_copy();
-                    const coeffT& t = result.get().second;
-
-                    qq = (parent_to_child(t, parent, key));
-                } else {
-                    qq = copy(it->second.coeff());
-                }
-                std::vector<Slice> s(NDIM*2);
-                long ll = 0;
-                for (std::size_t d=0; d<NDIM; ++d) {
-                    Translation l = key.translation()[d];
-                    long dum = long(float(l)/q);
-                    ll += (l - dum*q)*powMNDIM*powq[d] + dum*powM[d];
-                    //ll += (l % q)*powM[NDIM]*pow((double)q,NDIM-d-1) + (l/q)*pow((double)M,NDIM-d-1);
-
-                    //print("translation",l);
-                    //s[d       ] = Slice(l,l,0);
-                    //s[d+NDIM  ] = Slice(l%q,l%q,0);
-                    //s[d+NDIM] = Slice(0,k-1,1);
-                }
-                //long dum = ll;
-                for (std::size_t d=0; d<NDIM; ++d) {
-                    Translation l = Translation(float(ll) / powN[d]);
-                    //Translation l = ll / pow((double)N,NDIM-d-1);
-                    s[d ] = Slice(l,l,0);
-                    s[d+NDIM] = Slice(0,k-1,1);
-                    ll = ll - l*powN[d];
-                    //ll = ll % long(pow((double)N,NDIM-d-1));
-                }
-                //print(s, dum, key.translation());
-                coeffT qqq=qq(cdata.s0);
-                r(s) = qqq.full_tensor_copy();
-
-            }
-        }
-
-        world.gop.fence();
-        world.gop.sum(r0);
-        //print(r,r0);
-
-        return r0;
-    }
-
     /// truncate tree at a certain level
     template <typename T, std::size_t NDIM>
     void FunctionImpl<T,NDIM>::erase(const Level& max_level) {