def testlimA(self):
    j = 0
    sigma = math.sqrt(2*math.pi)*math.sqrt(0.5*self.P.intphi_over[j]/(math.pi**2))
    sigma_new = 2*self.P.intphi_over[j]/math.sqrt(2*math.pi)/self.P.intphi_under[j] # chiterstvo
    print "sigma(1) = ", self.P.sigma[j], sigma, sigma_new
    print "h = ", self.P.hval
    print "phi = ", self.P.intphi_over[j]
    limA = self.P.intphi_over[j]/((2*math.pi)**(1.5)*sigma)
    print "phi*limA = ", limA
    phi_r = self.P.intphi_under[j]/(4*math.pi)
    print "phi/(4pi*r) = ", phi_r
    print "*/* =", limA/phi_r
    print "sum = ", sum(self.P.intphi_under[:])/(4*math.pi)
    print "φ(r)/r full =", 0.691797810029
    print "φ(r)/r - sum =", sum(self.P.intphi_under[:])/(4*math.pi) - self.P.gauss[0,2]

self.centres[:] = self.quadphi_under[:,0]
self.C[:] = self.quadphi_under[:,1]
integ.calcomp2()

self.centres[:] = self.quadphi_over[:,0]
self.weights[:] = self.quadphi_over[:,1]
integ.calcomp3()

[self.intphi_over[:]] = self.withWrapMemo([self.axes,
                                           self.node_coordinates,
                                           self.dim_quad],
                                          'integ.calcomp()',
                                          [self.intphi_over[:]],
                                          self.flagMemo)

[self.gauss[:,3]] = self.withWrapMemo([self.axes,
                                       self.node_coordinates,
                                       self.dim_quad,
                                       self.k_wave],
                                      'integ.calcsing()',
                                      [self.gauss[:,3]],
                                      self.flagMemo)

def withWrapMemo(self,largs, body, larrs, flagMemo):
    @memoize
    def withWrapMemoInternal(largs,body):
        logging.info(body + ' is saved to memo')
        eval(body)
        return larrs
    if flagMemo:
        logging.info(body + ' uses memo')
        return withWrapMemoInternal(largs,body)
    else:
        logging.info(body + ' does not use memo')
        eval(body)
        return larrs

@@ -2057,28 +2057,1018 @@
 (body)%0A        return larrs%0A
+%0A    def norm(self,x):%0A        return (x%5B0%5D)**2 + (x%5B1%5D)**2 + (x%5B2%5D)**2%0A%0A    def delta_numbers(self):%0A        x0 = zeros((3))%0A%0A        smp = set()%0A        sbp = set()%0A        lm = 0%0A        lb = 0%0A        s = 0%0A        for j in range(0, self.numnodes,1):%0A            x0%5B:%5D = self.node_coordinates%5Bj,:%5D%0A            sb = set()%0A            sm = set()%0A            def trav(i):%0A                if self.norm(x0) %3E self.norm(self.node_coordinates%5Bi,:%5D):%0A                    sb.add(i)%0A                else:%0A                    sm.add(i)%0A%0A            map(trav, range(0,self.numnodes,1))%0A            if j%3E0:%0A                s = s + len(sm - smp)%0A                lm = max(lm,len(sm - smp))%0A                lb = max(lb,len(sb - sbp))%0A            smp = sm%0A            sbp = sb%0A        print lm, lb, 1.0*self.numnodes**2/s%0A        return lm, lb%0A%0A        # self.P.delta_numbers()%0A        # tock = datetime.now()%0A        # self.diff = tock - tick - self.diff%0A        # print %22seconds: %22, self.diff.seconds%0A

add delta_numbers