#*********************************************************************************
#
# Name: Peyton Capristo
# Pledge: Peyton Capristo...
#
# This program will perform a monte carlo simulation on a photon traveling 
# through the sun and find how long the photon takes to travel a certain distance
#
# Inputs: mfp = mean free path of a photon in meters
#         distances = a list containing different distance for the photon to travel
#
# Outpus: counter = the number of steps the photon took
#         time_elapsed = The average time taken per step
#
#*********************************************************************************


# Make an assumption that the typical mean free path of a photon is 2mm. 
# center of sun to radiation zone is 0.714 solar radii =(approx) 4.967x10^8 m 


# import libraries
import random as random
import math as math
import numpy as np
import matplotlib.pyplot as plt

c = 3.0*10**8
N = 100
distances = [0.01, 0.1, 1, 10]
total_times = []

for distance in distances:
    final = [] 
    for i in range(0, N): 
        # Initialize some variables 
        mfp = 0.002 
        position_tracker = 0
        x = 0
        y = 0
        z = 0
        counter = 0
# Do the random walk
        while position_tracker < distance:
            cos_phi = random.uniform(-1, 1)  # generates a random angle 
            sin_phi = math.sqrt(1 - cos_phi ** 2)
            theta = random.random() * 2 * math.pi   # generates a random angle 
            x = x + sin_phi * math.cos(theta)  * mfp  # multiply by a large number so the program runs shorter than 30,000-100,000 years.
            y = y + sin_phi * math.sin(theta) * mfp  # multiply by a large number so the program runs shorter than 30,000-100,000 years.
            z = z + cos_phi * mfp
            position_tracker = math.sqrt(x**2 + y**2 + z**2)
    
            counter += 1
        final.append(counter)
    avg_steps_taken = sum(final) / len(final)

    time_per_step = mfp / c
    time_elapsed = time_per_step * avg_steps_taken
    print("time_elapsed =", time_elapsed)
    total_times.append(time_elapsed)

# Set ups and displays the graph
plt.plot(distances, total_times) # Prepare a plot of time elapsed vs displacement 
plt.title('Time elapsed vs Displacement') # Give a title for the graph
plt.xlabel('Distance Travled (m)') # Give x axis label for the graph
plt.ylabel('Time Elapsed (s)') # Give y axis label for the graph
plt.grid(True, which='both')  # Show grid lines in both directions

plt.show() # Finally, show the plot