通过太阳随机漫步光子

Question

对于一个项目,我试图确定一个光子离开太阳所需的时间.但是,我的代码有问题(见下文).

更具体地说,我用if语句设置了一个for循环,如果一些随机生成的概率小于碰撞概率,那意味着光子碰撞并改变方向.

我遇到的问题是设置一个条件,如果光子逃逸(当距离>太阳半径时),for循环停止.我已经建立的那个似乎不起作用.

我使用太阳半径的非常小比例的测量,因为如果我不这样做,光子在我的模拟中需要很长时间才能逃脱.

    from numpy.random import random as rng # we want them random numbers
    import numpy as np # for the math functions
    import matplotlib.pyplot as plt # to make pretty pretty class

    mass_proton = 1.67e-27
    mass_electron = 9.11e-31
    Thompson_cross = 6.65e-29
    Sun_density = 150000
    Sun_radius = .005

    Mean_Free = (mass_proton + mass_electron)/(Thompson_cross*Sun_density*np.sqrt(2))

    time_step= 10**-13 # Used this specifically in order for the path length to be < Mean free Path
    path_length = (3e8)*time_step


    Probability = 1-np.exp(-path_length/Mean_Free) # Probability of the photon colliding


    def Random_walk():
        x = 0   # Start at origin (0,0)
        y = 0
        N = 1000
        m=0 # This is a counter I have set up for the number of collisions

        for i in range(1,N+1): 
            prand = rng(N+1)    # Randomly generated probability

            if prand[i] < Probability:  # If my prand is less than the probability
                                # of collision, the photon collides and changes
                                # direction
                x += Mean_Free*np.cos(2*np.pi*prand)    
                y += Mean_Free*np.sin(2*np.pi*prand)    
                m += 1  # Everytime a collision occurs 1 is added to my collision counter


            distance = np.sqrt(x**2 + y**2) # Final distance the photon covers

            if np.all(distance) > Sun_radius:   # if the distance the photon travels
                break                           # is greater than the Radius of the Sun,
                                        # the for loop stops, meaning the 
                                        #photon escaped

        print(m)    

        return x,y,distance


    x,y,d = Random_walk()
    plt.plot(x,y, '--')
    plt.plot(x[-1], y[-1], 'ro')

我欢迎任何批评我的代码,这是一个年级,我想知道如何正确地做到这一点,所以如果你发现任何其他错误,请告诉我.