third.py

   1 # -*- coding: utf-8 -*
   2
   3 from numpy import *
   4 from params import *
   5
   6 LLY = LY**2
   7 LLZ = LZ**2
   8 LLAMBDA = LAMBDA**2
   9 K = 2 * pi * C / LAMBDA
  10
  11 def setParameters(_ly, _lz, _c, _lambda):
  12     LY = _ly
  13     LLY = _ly**2
  14     LZ = _lz
  15     LLZ = _lz**2
  16     C = _c
  17     LAMBDA = _lambda
  18     LLAMBDA = _lambda**2
  19     K = 2 * pi * C / LAMBDA
  20
  21 def w(m):
  22     return pi * C * sqrt(1/LLY + m**2 / LLZ)
  23
  24 def psi(y, z, t):
  25     return (LZ - z) / LZ * sin(pi * y / LY) * sin(K * t)
  26
  27 def D(m):
  28     val = 2/(m * pi)
  29     val *= LLZ * (4 * LLY - LLAMBDA)
  30     val /= (LLZ * LLAMBDA + m**2 * LLY * LLAMBDA - 4 * LLY * LLZ)
  31     return val
  32
  33 def gamma(m, t):
  34     val = D(m) * sin(K*t)
  35     val -= K / w(m) * (2 / (m * pi) + D(m)) * sin(w(m) * t)
  36     return val
  37
  38 def um(y, z, t, m):
  39     return gamma(m, t) * sin(pi * y / LY) * sin(m * pi * z / LZ)
  40
  41 def Um(y, z, t, _M):
  42     v = 0.0
  43     for m in xrange(1, _M + 1):
  44         v += um(y, z, t, m)
  45     return v
  46
  47 def U(y, z, t):
  48     return Um(y, z, t, M)
  49
  50 def Ex(y, z, t):
  51     return U(y, z, t) + psi(y, z, t)
  52
  53 def Exm(y, z, t, _M):
  54     return Um(y, z, t, _M) + psi(y, z, t)
  55
  56 def remainder_length(eps):
  57     G = sqrt(4 * LY**2 - LAMBDA**2)/LAMBDA * (LZ / LY)
  58     M = LAMBDA * G**2 + 2 * LZ
  59     M += sqrt((LAMBDA * G**2 + 2 * LZ)**2 + 8 * pi * LAMBDA * LZ * eps * G**2)
  60     M /= (pi * LAMBDA * eps)
  61     if M < 20.0:
  62         M = 20
  63     return int(round(M + 0.49))
  64
  65 def analyze_remainder(M, eps, y, z, t):
  66     last = 0
  67     prev_sum = 0
  68     curr_sum = 0
  69     for m in xrange(1, M+1):
  70         curr_sum += um(y, z, t, m)
  71         if (abs(curr_sum - prev_sum) > eps):
  72             last = m
  73             prev_sum = curr_sum
  74     print prev_sum
  75     print curr_sum
  76     return last