# -*- coding: utf-8 -*-
import math
import numpy
from matplotlib import pyplot as plt
R = 145.5436
Sigma_t = 0.05
Sigma_s = 0.03
nv_Sigma_f = 0.0225
accuracy = 1e-5
# S4
mu = [-0.86113, -0.33998, 0.33998, 0.86113]
w = [0.34785, 0.65214, 0.65214, 0.34785]
# S6
#mu = [-0.93246, -0.66120, -0.23861, 0.23861, 0.66120, 0.93246]
#w = [0.17132, 0.36076, 0.46791, 0.46791, 0.36076, 0.17132]
# S12
#mu = [-0.9815606342, -0.9041172564, -0.7699026742, -0.5873179543, -0.3678314990, -0.1252334085, 0.1252334085, 0.3678314990, 0.5873179543, 0.7699026742, 0.9041172564, 0.9815606342]
#w = [0.0471753364, 0.1069393260, 0.1600783285, 0.2031674267, 0.2334925365, 0.2491470458, 0.2491470458, 0.2334925365, 0.2031674267, 0.1600783285, 0.1069393260, 0.0471753364]
# S16
#mu = [-0.98940, -0.94457, -0.86563, -0.75540, -0.61787, -0.45801, -0.28160, -0.09601, 0.09601, 0.28160, 0.45801, \
# 0.61787, 0.75540, 0.86563, 0.94457, 0.98940]
#w = [0.02715, 0.06225, 0.09515, 0.12462, 0.14959, 0.16915, 0.18260, 0.18945, 0.18945, 0.18260, 0.16915, 0.14959, \
# 0.12462, 0.09515, 0.06225, 0.02715]
K = 10000 # 段数
array_number = K + 1 #数组维数
delta_r = R / K
r = []
r_i = 0
# 定义初始量
for i in range(array_number):
r.append(r_i)
r_i += delta_r
i = 0
A = numpy.zeros([array_number])
while i < K:
i += 1
A[i] = r[i] * r[i]
i = 0
V = numpy.zeros([K])
while i < K:
V[i] = (r[i + 1]**3 - r[i]**3) / 3
i += 1
alpha = numpy.zeros([len(mu) + 1])
for i in range(len(mu)):
alpha[i + 1] = - (mu[i] * w[i]) + alpha[i]
k_avg = [0]
print r
print A
print V
print alpha[4]
S_his = []
def main():
a = numpy.ones([2 * len(mu) + 1,array_number + K])
S = numpy.ones([K])
phi = numpy.ones([2 * len(mu) + 1,array_number + K])
cycle = 0
S_his.append(S[1])
while 1:
# 算源
k = 0
while k < K:
sum = 0
j = 2
while j < (2 * len(mu) + 1):
sum = sum + phi[j - 1][2 * k + 1] * w[j / 2 - 1]
j += 2
S[k] = ((Sigma_s + nv_Sigma_f) / 2) * sum
k += 1
# 算通量
i = 0
k = K
phi[0][2 * K] = 0
while k > 0:
k -= 1
phi[0][2 * k + 1] = (V[k] * S[k] + (A[k + 1] + A[k]) * phi[0][2 * k + 2]) / (2 * A[k] + Sigma_t * V[k])
phi[0][2 * k] = 2 * phi[0][2 * k + 1] - phi[0][2 * k + 2]
for i in range(2,(2 * len(mu) + 1),2):
if (mu[i / 2 - 1] < 0):
k = K
while k > 1:
k -= 1
phi[i - 1][2 * k + 1] = (V[k] * S[k] + abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) * phi[i - 1][2 * k + 2] + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) * phi[i - 2][2 * k + 1] / w[i / 2 - 1]))\
/ (abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) / w[i / 2 - 1]) + Sigma_t * V[k])
phi[i - 1][2 * k] = 2 * phi[i - 1][2 * k + 1] - phi[i - 1][2 * k + 2]
phi[i][2 * k + 1] = 2 * phi[i - 1][2 * k + 1] - phi[i - 2][2 * k + 1]
if (mu[i / 2 - 1] > 0):
k = 0
phi[i - 1][0] = phi[2 * len(mu) - (i - 1)][0]
while k < K:
phi[i - 1][2 * k + 1] = (V[k] * S[k] + abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) * phi[i - 1][2 * k] + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) * phi[i - 2][2 * k + 1] / w[i / 2 - 1]))\
/ (abs(mu[i / 2 - 1]) * (A[k + 1] + A[k]) + ((A[k + 1] - A[k]) * (alpha[i / 2] + alpha[i / 2 - 1]) / w[i / 2 - 1]) + Sigma_t * V[k])
phi[i - 1][2 * k + 2] = 2 * phi[i - 1][2 * k + 1] - phi[i - 1][2 * k]
phi[i][2 * k + 1] = 2 * phi[i - 1][2 * k + 1] - phi[i - 2][2 * k + 1]
k += 1
cycle = cycle + 1
S_his.append(S[1])
print phi[3][K - 1]
k_avg.append(abs(S_his[cycle]/S_his[cycle - 1]))
if abs(S_his[cycle] - S_his[cycle - 1]) < accuracy:
print k_avg[cycle]
break
i = 0
#while i < (2 * len(mu) + 1):
'''
phi_new = numpy.zeros([2 * len(mu) + 1, K])
r_new = numpy.zeros([K])
for i in range(0, 2 * len(mu) + 1):
phi_new[]
for index in range(len(r) - 1):
r_new[index] = ((r[index] + r[index + 1]) / 2)
phi_new[i][index] = phi[i][2 * index + 1]
plt.plot(r_new,phi_new[i])
# i += 2
'''
phi_new = numpy.zeros([K])
r_new = numpy.zeros([K])
for j in range(0, K):
phi_new[j] = phi[5][2 * j + 1]
r_new[j] = (r[j] + r[j + 1]) / 2
plt.plot(r_new, phi_new)
plt.show()
if __name__ == '__main__':
main()
