Accelerated-Radiative-Diffusion/physics.py at main · hankyleh/Accelerated-Radiative-Diffusion · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
import numpy
import tools
import math

H = numpy.float128(4.135667696e-15) # eV * s
K = numpy.float128(8.617333262e-5)  # eV / K
C = numpy.float128(2.99792458e+10)   # cm / s
SIG_R = numpy.float128(3.53916934e+7) # eV / (cm^2 * s * K^4)
A_R = numpy.float128(4.72215928e-3)  # eV / (cm^3  * K^4)

ev_to_erg = (6.24150934e11)**-1

pi = numpy.pi


def cumulative_sigma(mesh, T, nu):
    # evaluates integral(B) d nu, from 0 to nu
    z = (H*nu)/(K*T)
    sigma =(
        (z <= 2) *
            ((z**3) * ((1/3) - (z/8) + ((z**2)/62.4))) +
        (z > 2)  *
            (6.4939 - (numpy.exp(-z) * ((z**3)+(3*(z**2))+(6*z)+7.28)))
    )
    # print(nu)
    # print("nu")
    # print(T)
    # print("T")
    # print(z)
    # print("z values")
    # print((z <= 2) + (z > 2))
    # print("valid z")
    return sigma

def group_planck(mesh, T):
    # evaluates integral(B) d nu, along each group bound
    # "beta" is the planck function integrated over 4pi,
    #  beta = 4pi * (2h nu^3 / c^2)/(e^h nu/kT - 1)
    groups = mesh.groups[:, numpy.newaxis].reshape((mesh.ng+1, 1))

    if numpy.size(T) != 1:
        T.reshape((1, mesh.nx))
        T[T<0] = K*mesh.eps
        t = numpy.tile(T, (mesh.ng+1, 1))
    else:
        t = numpy.tile(T, (mesh.ng+1, mesh.nx))
    g = numpy.tile(groups, (1, mesh.nx))

    result = (4*pi)*(1/pi) *numpy.diff(cumulative_sigma(mesh, t, g), axis=0) * t[1:]**4 * ((2*pi*(K)**4)/(H**3*C**2))

    if (sum(math.isnan(result[i,j]) for i in range(0, mesh.ng) for j in range(0, mesh.nx)) > 0):
        print(result[:, 0:3])
        print("Planck result sample")
        print(T)
        print("Temperature")
        # print(t[1:]**4)
        # print("T^4")
        print(g)
        print("groups")

        print(cumulative_sigma(mesh, t, g))
        print("Sigma values")
        raise ValueError("Nan Planck encountered")


    return result

def cumulative_dsigma_dT(mesh, T, nu):
    z = (H*nu)/(K*T)
    dzdt = -(H*nu)/(K*(T**2))
    d = (
        (z <= 2) *
            (((3*z**2) * ((1/3) - (z/8) + ((z**2)/62.4)))
             + ((z**3) * ((-1/8) + ((2*z)/62.4)))) +
        (z > 2)  *
            ((numpy.exp(-z) * ((z**3)+(3*(z**2))+(6*z)+7.28)) -
             (numpy.exp(-z) * ((3*z**2)+(6*z)+(6))))
    )
    return d*dzdt

def group_dB_dT(mesh, T):

    groups = mesh.groups[:, numpy.newaxis].reshape((mesh.ng+1, 1))
    if numpy.size(T) != 1:
        T.reshape((1, mesh.nx))
        T[T<0] = K*mesh.eps
        t = numpy.tile(T, (mesh.ng+1, 1))
    else:
        t = numpy.tile(T, (mesh.ng+1, mesh.nx))
    g = numpy.tile(groups, (1, mesh.nx))

    result = (

        (4*pi)*(1/pi)*((2*pi*(K)**4)/(H**3*C**2))*(
            (t[1:]**4)*numpy.diff(cumulative_dsigma_dT(mesh, t, g), axis=0)
            +
            (4*t[1:]**3)*numpy.diff(cumulative_sigma(mesh, t, g), axis=0)
        )
    )


    if (sum(math.isnan(result[i,j]) for i in range(0, mesh.ng) for j in range(0, mesh.nx)) > 0):
        print(result[:, 0:3])
        print("dB/dT result sample")
        print(T)
        print("Temperature")
        print(numpy.diff(cumulative_sigma(mesh, t, g), axis=0))
        print("group sigma")

        print(cumulative_sigma(mesh, t, g))
        print("Sigma values")


        raise ValueError("Nan dB/dT encountered")

    return result