Случайная прогулка фотона по солнцу

Question

Для проекта я пытаюсь определить, сколько времени потребуется фотону, чтобы покинуть Солнце. Однако у меня возникли проблемы с моим кодом (см. Ниже).

Точнее говоря, я настроил цикл for с оператором if, и если какая-то случайно сгенерированная вероятность меньше вероятности столкновения, это означает, что фотон сталкивается, и он меняет направление.

У меня возникли проблемы с установкой условия, при котором цикл 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')

Любая критика моего кода приветствуется, это для оценки, и я хочу научиться делать это правильно, поэтому, пожалуйста, сообщите мне, если вы заметили какие-либо другие ошибки.

Answer 1

Вот способ размышления о проблеме, который может оказаться для вас полезным. В каждый момент фотон имеет положение (x, y) и направление (dx, dy). Переменные (dx, dy) являются коэффициентами единичного вектора, поэтому sqrt(dx**2 + dy**2) = 1.0. Расстояние, пройденное фотоном за один шаг, составляет направление path_length *.

На каждом шаге вы делаете 4 вещи:

1. рассчитать новую позицию фотона
2. выяснить, покинул ли фотон Солнце, вычислив его расстояние от центральной точки
3. определить с помощью одного случайного числа, сталкивается ли фотон или нет. Если это так, вы случайным образом генерируете новое направление.
4. Добавление текущей позиции фотона в список. Возможно, вы захотите сделать это как функцию расстояния от центра, а не х, у.

В конце вы строите построенный список.

Вы должны также выбрать случайное направление в самом начале.

Я не знаю, как вы прервете цикл, потому что фотон никогда не гарантированно покинет солнце - как в реальном мире. В принципе, программа может работать вечно (или до тех пор, пока не сгорит солнце).

Существует небольшая неточность в том, что фотон может столкнуться в любой момент, а не только в конце одного шага. Но поскольку шаги невелики, то и ошибка, вносимая этим упрощением.

Я укажу, что вам ничего не нужно для этого, кроме, возможно, заключительного сюжета. Стандартная библиотека Python имеет все математические функции, которые вам нужны. Numpy, конечно, отлично подходит для манипулирования массивами данных, но единственный массив, который у вас здесь будет, это тот, который вы строите, шаг за шагом, положение фотона в зависимости от времени.

Как я указывал в одном из моих комментариев, вы моделируете солнце как двумерный объект. Если вы хотите выполнить этот расчет в трех измерениях, вам не нужно менять этот базовый подход.

Answer 2

Я не понимаю мотивацию для формул, которые вы внедрили.Я объясню свою собственную мотивацию здесь, но если ваш инструктор сказал вам делать что-то еще, я думаю, вы должны вместо этого их слушать.

Если бы я собирался это сделать, я бы сгенерировал последовательность движенийфотона, останавливаясь, когда расстояние фотона до центра солнца больше солнечного радиуса.Каждое движение является выборкой из распределения, которое имеет два компонента: один для расстояния и один для направления.Я предполагаю, что они независимы (это может быть поставлено под сомнение при более тщательном моделировании).

Представляется вероятным, что распределение расстояния является экспоненциальным распределением с параметром 1 / (средняя длина свободного пробега).Тогда плотность будет p (d) = (1 / MFP) exp (-d / MFP).Его cdf равен 1 - exp (-d / MFP), а обратный к cdf - журнал -MFP (1 - p), где p = cdf (d).Теперь вы можете сделать выборку из распределения расстояний: пусть p = rand (0, 1), где rand = равномерное случайное число, и вставьте его в обратный cdf, чтобы получить d.Это называется методом обратной выборки cdf;веб-поиск найдет больше информации об этом.

Что касается направления, вы можете указать angle = rand (0, 2 * pi), а затем (x, y) = (cos (angle), sin (угол)).

Теперь вы можете построить ряд позиций.Из исходного местоположения пусть новое местоположение = предыдущий + d * (x, y).Остановитесь, когда расстояние от центра до радиуса больше радиуса.

Выглядит как большая проблема!Удачи и приятного времяпровождения.Дайте мне знать, если у вас есть какие-либо вопросы.