<jupyter_start><jupyter_text>define density profiles: density_delta used hereafter<jupyter_code>def density_alpha(radius):
    p0 = 1000
    p = p0*(radius/radius)
    return p
def density_delta(radius, r_star):
    p0 = 100
    frac_radius = r/r_star
    p = p0*(1-frac_radius**2)
    return p
r = np.arange(1e-4, 6.957e10+100000, 100000)
density_d = density_delta(r, 6.957e10)
plt.plot(r, density_d, label = '$\\rho_{\delta}$')
plt.xlabel('r[cm]')
plt.ylabel('$\\rho$ [$\\frac{g}{cm^3}$]')
plt.show()<jupyter_output><empty_output><jupyter_text>derive pressure profiles using equations of stellar structure<jupyter_code>G = 6.67e-8
k = 1.38e-16
mu = 1.3
mp = 1.67e-24

def pressure_alpha(radius, density, M_star, R_star):
    P_c = (2/3)*np.pi*G*(density)**2*R_star**2
    P_r = P_c*(1-(radius**2)/(R_star**2))
    
    return P_r

def pressure_delta(radius, central_density, R_star): 
    """
    P_c = (176/225)*np.pi*G*R_star*central_density**2
    
    P_r = P_c - ((4*np.pi*G*(central_density**2))*((radius/3)-((8*radius**3)/(45*R_star**2))+((radius**5)/(25*R_star**4))))
    """
    P_c = (4/15)*np.pi*G*(central_density**2)*(R_star**2)
    term1 = (r**2)/6
    term2 = (2*r**4)/(15*R_star**2)
    term3 = (r**6)/(30*R_star**4)
    
    P_r = P_c - (4*np.pi*G*(central_density**2)*(term1-term2+term3))
    
    return P_r
r = np.arange(1e-4, 6.957e10+100000, 100000)
pressure_d = pressure_delta(r, 100, 6.957e10)
plt.plot(r,pressure_d)
plt.xlabel('r[cm]')
plt.ylabel('P [dyn $cm^{-2}$]')
plt.show()<jupyter_output><empty_output><jupyter_text>temperature profile - derived using equation of state<jupyter_code>def temp_alpha(density, pressure):

    T_r = (mu*mp*pressure)/(k*density)
    
    return T_r

temp = temp_alpha(density_d, pressure_d)

plt.plot(r,temp)
plt.xlabel('r[cm]')
plt.ylabel('T[K]')
#plt.xlim(0, 0.99e11)
plt.show()<jupyter_output><empty_output><jupyter_text>Different methods of defining, cross section of interaction, and opacity. rosseland_mean_opacity function from Carroll and Ostile used to find mean free path.<jupyter_code>def rosseland_path(density, temperature):
    #gives cross section of interaction
    X = 0.7381 #solar value
    Z = 0.0134 #solar value
    g_ff = 1.0
    C = 6.25E-2
    K_r = C*(1+X)*(Z**2)*g_ff*(density)*((temperature)**(-3.5))
    return K_r

def free_free_opacity(density, temperature):
    C_ff = 4e22
    gff = 1
    X = 0.7381
    Y = 0.2485
    Z = 0.0134
    K_ff = C_ff*(X+Y)*(1+X)*(density/(temperature**3.5))
    return K_ff

def rosseland_mean_opacity(density, temperature):
    X = 0.7381
    Y = 0.2485
    Z = 0.0134
    
    C_ff = 3.68e18
    gff = 1
    
    K_ff = C_ff*(X+Y)*(1+X)*(density/(temperature**3.5))
    
    C_bf = 4.34e21
    g_bf = 1
    t = 1
    
    K_bf = C_bf*(g_bf/t)*Z*(1+X)*(density/(temperature**3.5))
    
    C_es = 0.02
    K_es = C_es*(1+X)
    
    #we're able to ignore H- opacity
    """
    if temperature <= 6000 and density <= 1e-5:
        C_H_ = 7.9e-34
        K_H_ = C_H_*(Z/0.02)*(density**-0.5)*(temperature**9)
    elif temperature >= 3000 and density >= 1e-10:
        C_H_ = 7.9e-34
        K_H_ = C_H_*(Z/0.02)*(density**-0.5)*(temperature**9)
    else:
        K_H_ = 0
    
    #K = np.mean([K_ff, K_bf, K_es]) 
    #K = ((K_ff**-1)+(K_bf**-1)+(K_es**-1))**-1
    """
    K = K_ff+K_bf+K_es
    
    return K
   
r = np.arange(1e-4, 6.957e10, 100000)
density_interaction = density_delta(r, 6.957e10)
pressure_interaction = pressure_delta(r, 100, 6.957e10)
temp_interaction = temp_alpha(density_interaction, pressure_interaction)
interaction = rosseland_path(density_interaction, temp_interaction)

plt.semilogy(r, interaction)
plt.xlabel('r[m]')
plt.ylabel('$\\sigma$ [$cm^2$]')
plt.show()

def number_density(density):
    n = density/mp
    return n<jupyter_output><empty_output><jupyter_text>fig = plt.figure(figsize = (10,10))ax = plt.axes(projection='3d')a = 0.75b = 0.45c = 0.25u, v = np.meshgrid(np.arange(0,2*np.pi, 0.01), np.arange(0,np.pi,0.005))    ellipsoid_x = (a * np.cos(u) * np.sin(v)) + 0.2ellipsoid_y = (b * np.sin(u) * np.sin(v)) - 0.4ellipsoid_z = (c * np.cos(v)) - 0.1ax.plot_wireframe(ellipsoid_x, ellipsoid_y, ellipsoid_z)ax.set_xlim(-a,a)ax.set_ylim(-a,a)ax.set_zlim(-a,a)plt.show() mean free path can de defined in different ways, (1) $l = \frac{1}{\sigma * n}$, or (2) $l = \frac{1}{\kappa * \rho}$. I chose to use definition(2).<jupyter_code>def mean_free_path(cross_section, n_density):
    l = 1/(cross_section*n_density)
    return l
def path_probability(l, dtheta, dphi, e):
    r_theta = ((l/(1+e))*(1-e**2))/(1-e*np.cos(dtheta))
    r_phi = ((l/(1+e))*(1-e**2))/(1-e*np.cos(dphi))
    l = np.mean([r_theta, r_phi])
    #l = r_theta
    return l
def mean_free_path_beta(opacity, density):
    l = 1/(density*opacity)
    return l

n_density = number_density(density_interaction)
path_length = mean_free_path(interaction, n_density)
plt.plot(r, path_length)
plt.ylabel('l[cm]')
plt.xlabel('r[cm]')
plt.show()

path_length_thom = mean_free_path(thomson_interaction, n_density)
plt.semilogy(r, path_length_thom)
plt.ylabel('l[cm]')
plt.xlabel('r[cm]')
plt.show()

opacity_interaction = free_free_opacity(density_interaction, temp_interaction)
path_length = mean_free_path_beta(opacity_interaction, density_interaction)
plt.plot(r, path_length)
plt.show()

array = np.column_stack((density_interaction, temp_interaction))


mean_path = []
for row in array: 
    rosseland_opacity = rosseland_mean_opacity(row[0], row[1])
    path = mean_free_path_beta(rosseland_opacity, row[0])
    mean_path.append(path)
plt.semilogy(r, mean_path)
plt.show()<jupyter_output><empty_output><jupyter_text>Now we have the work-horses. Below are different methods I tried to have the photon randomly walk out of the sun, and determine how long that process takes.<jupyter_code>#method (1): Purely random walk
N_steps = 1e6
N_layers = 10
N_walkers = 1000
c = 3e10

x = 0
y = 0
z = 0

x_val = [0]
y_val = [0]
z_val = [0]

r_star = 6.957e10

section = r_star/N_layers

epsilon= 1e-4

radius_bisection = np.arange(epsilon, r_star, section)
t_est = []
median_velocity = []

for radius in radius_bisection:
    f = open('rand_{:.0f}_layer.txt'.format(float(radius/section)), 'a+')
    for walker in range(N_walkers):
        x = radius
        y = 0
        z = 0
        t = 0
        r = radius
        #x_walk = [x]
        #y_walk = [y]
        #z_walk = [z]
        l_all = []
        for n in range(N_steps):
            
            density = density_delta(r, r_star)
            pressure = pressure_delta(r, 100, r_star)
            temperature = temp_alpha(density, pressure)
            opacity = rosseland_mean_opacity(density, temperature)

            theta = rand.uniform(0,2*np.pi)
            phi = rand.uniform(0,2*np.pi)
            
            l = mean_free_path_beta(opacity, density)
            #l = path_probability(l_tot, theta, phi, 0.85)
            l_all.append(l)
            
            dx = l * np.sin(phi) * np.cos(theta)
            dy = l * np.sin(phi) * np.sin(theta)
            dz = l * np.cos(phi)
           
            x = x + dx
            y = y + dy
            z = z + dz
            #x_walk.append(x)
            #y_walk.append(y)
            #z_walk.append(z)
            
            r = np.sqrt(x**2 + y**2 + z**2)

            d = r - radius

            dt = l/c
            
            t = t + dt
        """
        fig = plt.figure(figsize = (10,10))
        ax = plt.subplot(111, projection='3d')
        ax.plot3D(x_walk,y_walk,z_walk, '.', color = 'orange')
        ax.scatter3D(x_walk[0],y_walk[0],z_walk[0], color = 'green')
        ax.scatter3D(x_walk[-1],y_walk[-1],z_walk[-1], '.', color = 'red')
        ax.plot3D(x_walk,y_walk,z_walk, color = 'orange')
        plt.savefig('example_photon_walk_rand.pdf')
        plt.show()
        """
        l_all = np.array(l_all)
        np.savetxt(f, np.array([d,t, np.mean(l_all)]))
    f.close()


N_steps = 100
N_layers = 10
N_walkers = 100
c = 3e10

x = 0
y = 0
z = 0

x_val = [0]
y_val = [0]
z_val = [0]

r_star = 6.957e10

section = r_star/N_layers

epsilon= 1e-4

radius_bisection = np.arange(epsilon, r_star, section)
t_est = []
median_velocity = []

for radius in radius_bisection: 
    f = open('dot_prod_{:.0f}_layer.txt'.format(float(radius/section)), 'w+')
    
    for walker in range(N_walkers):
        
        x = radius
        y = 0
        z = 0
        t = 0
        r = radius
        dt = 0
        tot_l = 0
        for n in range(N_steps):
            
            density = density_delta(r, r_star)
            pressure = pressure_delta(r, 100, r_star)
            temperature = temp_alpha(density, pressure)
            opacity = rosseland_mean_opacity(density, temperature)

            theta = rand.uniform(0,2*np.pi)
            phi = rand.uniform(0,2*np.pi)
            
            l = mean_free_path_beta(opacity, density)
            
            tot_l_new = tot_l + l
            
            mean_l = tot_l_new/(n+1)
            
            tot_l = tot_l_new

            dx = l * np.sin(phi) * np.cos(theta)
            dy = l * np.sin(phi) * np.sin(theta)
            dz = l * np.cos(phi)
           
            x = x + dx
            y = y + dy
            z = z + dz
            
            r = np.sqrt(x**2 + y**2 + z**2)

            d = r - radius

            dt_new = l/c
            
            mean_dt = (t + dt)/(n+1)
            
            dt = dt_new
            t = t + dt
        
        data = np.array([mean_l, mean_dt])
        f.write(str(data)+ "\n")
        
    f.close()
#method (2): dot product walk

N_steps = 1e6
N_layers = 10
N_walkers = 1000
c = 3e10

x = 0
y = 0
z = 0

x_val = [0]
y_val = [0]
z_val = [0]

r_star = 6.957e10

section = r_star/N_layers

epsilon= 1e-4

radius_bisection = np.arange(epsilon, r_star, section)

t_est = []
median_velocity = []

for radius in radius_bisection: 
    f = open('dot_product_{:.0f}_layer.txt'.format(float(radius/section)), 'w+')
    for walker in range(N_walkers):
        #w = open('walker_data/walker{:.0f}_layer{:.0f}.txt'.format(walker, float(radius/section)), 'w+')
        x = radius/np.sqrt(3)
        y = radius/np.sqrt(3)
        z = radius/np.sqrt(3)
        norm_vector = np.array([1/np.sqrt(3),1/np.sqrt(3),1/np.sqrt(3)])
        #coordinate = np.array([x,y,z])
        r0 = radius
        t = 0
        r = radius
        l_all = []
        #x_walk = [x]
        #y_walk = [y]
        #z_walk = [z]
        for n in range(N_steps):
            density = density_delta(r, r_star)
            pressure = pressure_delta(r, 100, r_star)
            temperature = temp_alpha(density, pressure)
            opacity = rosseland_mean_opacity(density, temperature)
            l = mean_free_path_beta(opacity, density)
            l_all.append(l)
            theta = rand.uniform(0,2*np.pi)
            phi = rand.uniform(0,2*np.pi)
 
            dx = l * np.sin(phi) * np.cos(theta)
            dy = l * np.sin(phi) * np.sin(theta)
            dz = l * np.cos(phi)
            
            dr = np.array([dx,dy,dz])
            
            r_add = r + (np.dot(dr, norm_vector))
            
            r = r_add 
            
            #x = r/np.sqrt(3)
            #y = r/np.sqrt(3)
            #z = r/np.sqrt(3)
            #x_walk.append(x)
            #y_walk.append(y)
            #z_walk.append(z)
            
            #coordinate = np.array([r/np.sqrt(3), r/np.sqrt(3), r/np.sqrt(3), l, dt])
            #np.savetxt(w, coordinate)
            
            dt = l/c
  
            t = t + dt
        #w.close()
        d = r - radius
        l_mean = np.mean(np.array(l_all))
        """
        fig = plt.figure(figsize = (10,10))
        ax = plt.subplot(111, projection='3d')
        ax.plot3D(x_walk,y_walk,z_walk, '.', color = 'orange')
        ax.scatter3D(x_walk[0],y_walk[0],z_walk[0], color = 'green')
        ax.scatter3D(x_walk[-1],y_walk[-1],z_walk[-1], '.', color = 'red')
        ax.plot3D(x_walk,y_walk,z_walk, color = 'orange')
        plt.savefig('example_photon_walk_dot.pdf')
        plt.show()
        """
        np.savetxt(f, np.array([d,t, l_mean]))
    f.close()
all_layers =  sorted(glob.glob("rand*.txt"))
tau = []
tau_beta = []
avg_t = []
avg_d = []
avg_v = []
avg_l = []
for layer in all_layers:
    data = np.genfromtxt(layer)
    n_layer = layer.split('_')[1]
    reshaped = data.reshape(N_walkers,3)
    d_walkers = reshaped[:,0]
    t_walkers = reshaped[:,1]
    l_walkers = reshaped[:,2]
    avg_d_layer = np.mean(d_walkers)
    avg_t_layer = np.mean(t_walkers)
    avg_l_layer = np.mean(l_walkers)
    
    plt.title('distance walker traveled layer {0}'.format(n_layer))
    plt.hist(d_walkers, bins = 20)
    plt.show()
    plt.title('time walker traveled layer {0}'.format(n_layer))
    plt.hist(t_walkers, bins = 20)
    plt.show()
    plt.title('average mfp walker layer {0}'.format(n_layer))
    plt.hist(l_walkers, bins = 20)
    plt.show()
    
    avg_v_layer = avg_d_layer/avg_t_layer
    avg_d.append(avg_d_layer)
    avg_t.append(avg_t_layer)
    avg_v.append(avg_v_layer)
    avg_l.append(avg_l_layer)
    
    tau_section = section/avg_v_layer
    tau.append(tau_section)
    
    N_steps_needed = (section/avg_l_layer)**2
    tau_section_beta = N_steps_needed*avg_t_layer
    tau_beta.append(tau_section_beta)
plt.title('distance per layer')   
plt.plot(np.arange(0,10,1), avg_d, '.')
plt.xlabel('layer')
plt.ylabel('avg. distance walker travelled [cm]')
plt.savefig('rand_distance_layer.pdf')
plt.show()
plt.title('time spent per layer')
plt.plot(np.arange(0,10,1), avg_t, '.')
plt.xlabel('layer')
plt.ylabel('avg. time walker spent [s]')
plt.savefig('rand_time_layer.pdf')
plt.show()
tau = np.array(tau)
tau_beta = np.array(tau_beta)
escape_time=np.sum(tau)
escape_time_beta = np.sum(tau_beta)
print('photon escape time = {0}[s], from random walk'.format(escape_time))
print('photon escape time = {0}[yr], from random walk'.format(escape_time/(np.pi*1e9)))
print('photon escape time = {0}[s], from random walk'.format(escape_time_beta))
print('photon escape time = {0}[yr], from number steps'.format(escape_time_beta/(np.pi*1e9)))

f = open('results_rand_walk.txt', '+w')
f.write('photon escape time = {0}[s], from random walk'.format(escape_time)+'\n'+'photon escape time = {0}[yr], from random walk'.format(escape_time/(np.pi*1e9))+'\n'+'photon escape time = {0}[s], from random walk'.format(escape_time_beta)+'\n'+'photon escape time = {0}[yr], from number steps'.format(escape_time_beta/(np.pi*1e9)))
f.close()
all_layers = sorted(glob.glob("dot_product*.txt"))
tau = []
tau_beta = []
avg_t = []
avg_d = []
avg_v = []
avg_l = []
for layer in all_layers:
    data = np.genfromtxt(layer)
    n_layer = layer.split('_')[2]
    reshaped = data.reshape(N_walkers,3)
    d_walkers = reshaped[:,0]
    t_walkers = reshaped[:,1]
    l_walkers = reshaped[:,2]
    avg_d_layer = np.mean(d_walkers)
    avg_t_layer = np.mean(t_walkers)
    avg_l_layer = np.mean(l_walkers)
    
    plt.title('distance walker traveled layer {0}'.format(n_layer))
    plt.hist(d_walkers, bins = 20)
    plt.show()
    plt.title('time walker traveled layer {0}'.format(n_layer))
    plt.hist(t_walkers, bins = 20)
    plt.show()
    plt.title('average mfp walker layer {0}'.format(n_layer))
    plt.hist(l_walkers, bins = 20)
    plt.show()
    
    avg_v_layer = avg_d_layer/avg_t_layer
    avg_d.append(avg_d_layer)
    avg_t.append(avg_t_layer)
    avg_v.append(avg_v_layer)
    avg_l.append(avg_l_layer)
    
    tau_section = section/avg_v_layer
    tau.append(tau_section)
    
    N_steps_needed = (section/avg_l_layer)**2
    tau_section_beta = N_steps_needed*avg_t_layer
    tau_beta.append(tau_section_beta)
plt.title('distance per layer')   
plt.plot(np.arange(0,10,1), avg_d, '.')
plt.xlabel('layer')
plt.ylabel('avg. distance walker travelled [cm]')
plt.savefig('dot_distance_layer.pdf')
plt.show()
plt.title('time spent per layer')
plt.plot(np.arange(0,10,1), avg_t, '.')
plt.xlabel('layer')
plt.ylabel('avg. time walker spent [s]')
plt.savefig('dot_time_layer.pdf')
plt.show()
tau = np.array(tau)
tau_beta = np.array(tau_beta)
escape_time=np.sum(tau)
escape_time_beta = np.sum(tau_beta)
print('photon escape time = {0}[s], from random walk'.format(escape_time))
print('photon escape time = {0}[yr], from random walk'.format(escape_time/(np.pi*1e9)))
print('photon escape time = {0}[s], from random walk'.format(escape_time_beta))
print('photon escape time = {0}[yr], from number steps'.format(escape_time_beta/(np.pi*1e9)))

f = open('results_dot_walk.txt', '+w')
f.write('photon escape time = {0}[s], from random walk'.format(escape_time)+'\n'+'photon escape time = {0}[yr], from random walk'.format(escape_time/(np.pi*1e9))+'\n'+'photon escape time = {0}[s], from random walk'.format(escape_time_beta)+'\n'+'photon escape time = {0}[yr], from number steps'.format(escape_time_beta/(np.pi*1e9)))
f.close()<jupyter_output>photon escape time = -694.3552763578667[s], from random walk
photon escape time = -2.2102014898858707e-07[yr], from random walk
photon escape time = 4009912293167.525[s], from random walk
photon escape time = 1276.394725645138[yr], from number steps
StarcoderdataJupyterStructuredCleanDedup