<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Monte Carlo Simulations (Section 8.3 in Text)
Monte Carlo simulations are used to calculate probabilities using random numbers when the probabilities are difficult (or impossible) to calculate by hand or have to be done many times.
Basic idea
Step1:  Q. How would you estimate the time it takes for a photon to travel from the core to the surface?
Step2:  This is a 3D random walk (if you take ASTR3730 you'll see a derivation of this).
http
Step3:  Let's compute tau by using a Monte Carlo experiment.
Area of enclosed circe in a square vs that area of the square is
Step4:  Q. Does "positions" contain the final position of each particle, or the entire trajectory of each particle?
Step5:  Q. How can we get the average position of all particles?
Step6:  Q. And the average traveled distance for all particles?
Step7:  $\sqrt{\frac{4N_s}{\tau}}$ is the theoretical expectation value for the separation for a large number of tests run, in 1D. (different values for higher dimensions).
Step8:  Q. What should this histogram look like? Should it be centrally peaked? If so, at what value? How wide?
Step9:  VECTORIZATION of Implementation
NOTE the difference in arguments!
Inclusive right limit
random.randint(1, 2)
Exclusive right limit
np.random.randint(1, 3)
Step10:  Files
Basic file operation works via the open() function
Step11:  Use so called context manager
<ASSISTANT_TASK:>
Python Code:
from IPython.display import Image
Image(url='http://upload.wikimedia.org/wikipedia/commons/thumb/b/b4/The_Sun_by_the_Atmospheric_Imaging_Assembly_of_NASA%27s_Solar_Dynamics_Observatory_-_20100819.jpg/251px-The_Sun_by_the_Atmospheric_Imaging_Assembly_of_NASA%27s_Solar_Dynamics_Observatory_-_20100819.jpg')
Explanation: Monte Carlo Simulations (Section 8.3 in Text)
Monte Carlo simulations are used to calculate probabilities using random numbers when the probabilities are difficult (or impossible) to calculate by hand or have to be done many times.
Basic idea:  Perform N experiments (numerically), where the outcome is random (or at least unknown), with some event occuring M times, or with a probability M/N each time.  If the experiment is performed many times, M/N yields an estimate of the probability.
There are many examples in astrophysics:
Scattering of a photon by atoms or molecules
Detection of a weak signal (say a galaxy) in the 
   presence of noise (say sky noise and electronics noise)
Observation of a particular power law in galaxy number  
   counts given a model power law in the presence of noise
Markov Chain Monte Carlo for determining models
   that best describe data
Q. Have you encountered any other examples?
End of explanation
Image(url='http://upload.wikimedia.org/wikipedia/commons/thumb/d/d4/Sun_poster.svg/500px-Sun_poster.svg.png')
Explanation: Q. How would you estimate the time it takes for a photon to travel from the core to the surface?
End of explanation
%matplotlib inline
import numpy as np
import matplotlib.pylab as pl
Explanation: This is a 3D random walk (if you take ASTR3730 you'll see a derivation of this).
http://en.wikipedia.org/wiki/Random_walk has a good intro for some of the maths involved.
Random Walk in One Dimension
For simplicity, we will assume the same distance between scatterings 
(and set that distance, or the "mean free path", equal to 1).
Q. What do you think the typical distance a particle will have traveled  (in net distance from its starting point) after a large number of steps ($N_s$)?
In fact, for large $N_s$ the average distance traveled will be $\sqrt{N_s}$.
Q. Put another way, on average, what is the final position of a particle  undergoing a 1D random walk with respect to its starting point?
End of explanation
pl.plot?
N = 10000
radius = 250
tau = 2*np.pi
x = np.random.randint(-radius, radius, N)
y = np.random.randint(-radius, radius, N)
distances = np.sqrt(x**2 + y**2)
phi = np.arange(0, tau, 0.01)
pl.figure(figsize=(5,5))
pl.scatter(x[np.argwhere(distances <= radius)],
           y[np.argwhere(distances <= radius)],
           color='k', s=1);
pl.scatter(x[np.argwhere(distances > radius)],
           y[np.argwhere(distances > radius)],
           color='r', s=1);
pl.plot(radius * np.cos(phi), radius * np.sin(phi), color='g', linewidth=2);
print("tau is", 8.0 * np.sum(distances < radius) / float(N))
tauArray = np.zeros(N)
for num in np.arange(N):
    xArray = np.random.randint(-radius, radius, num + 1)
    yArray = np.random.randint(-radius, radius, num + 1)
    tauArray[num] = 8.0 * np.sum(np.sqrt(xArray**2 + yArray**2) < radius) / float(num + 1)
pl.plot(np.arange(N), tauArray);
pl.axhline(y=tau, lw=2, color='r')
pl.ylim(tau - 0.5, tau + 0.5);
import random
# Number of particles
Np = 100
# Number of steps (per particle)
Ns = 50000
# All particles start at x = 0
positions = np.zeros(Np)  
distances = np.zeros(Np)
# A (randomly drawn) 1 will move the particle to the left
# and a 2 will move it to the right
Left = 1; Right = 2  
# Step Ns times for each particle "p"
for p in range(Np):
    for step in range(Ns):
        
        # Integer random number generator
        direction = random.randint(1, 2) 
        # returns a random integer x such that 1 <= x <= 2
        # (effectively a coin-flip here)
        
        if direction == Left:
            positions[p] -= 1  # Move left
        elif direction == Right:
            positions[p] += 1  # Move right
Explanation: Let's compute tau by using a Monte Carlo experiment.
Area of enclosed circe in a square vs that area of the square is:
$$ \frac{\frac{1}{2}\tau r^2}{\left(2r\right)^2} = \frac{\tau}{8}$$
This means, if we could somehow find the actual ratio of these areas, we could multiply it with 8 to get the value of $\tau$!
End of explanation
print("Positions")
print(positions)
Explanation: Q. Does "positions" contain the final position of each particle, or the entire trajectory of each particle?
End of explanation
print("Average Position", positions.mean())
Explanation: Q. How can we get the average position of all particles?
End of explanation
print("Avg Separation", abs(positions).mean())
print("Expectation    %g" % np.sqrt(4*Ns/tau))
Explanation: Q. And the average traveled distance for all particles?
End of explanation
# Standard deviation
positions.std()
Explanation: $\sqrt{\frac{4N_s}{\tau}}$ is the theoretical expectation value for the separation for a large number of tests run, in 1D. (different values for higher dimensions).
End of explanation
n, bins, patches = pl.hist(positions, 20, facecolor='g')
pl.xlabel('Final Position')
pl.ylabel('Frequency')
pl.title('1-D Random Walk Distance')
pl.grid(True)
Explanation: Q. What should this histogram look like? Should it be centrally peaked? If so, at what value? How wide?
End of explanation
Np = 100    # Number of particles
Ns = 50000  # Number of steps (per particle)
# Draw the move random number steps all at once:
moves = np.random.randint(1, 3, size=Np*Ns)
# FOR np.random.randint THIS RUNS FROM 1 TO 2, INTEGERS ONLY
# Q. What's happening here?
moves = 2 * moves - 3
# Create a 2-D array of moves so that moves[i, j]
# is the "i"th step of particle j:
moves.shape = (Ns, Np)
# Create an array of initial starting positions for each particle
positions = np.zeros(Np)
for step in range(Ns):
    # Select the moves values for the current step:
    positions += moves[step, :]
    # Updates positions for all particles in this step
    # This is vectorized:  I'm not looping over the particles
# Histogram the results
n, bins, patches = pl.hist(positions, bins=np.arange(-500, 500, 50), facecolor='b')
pl.xlabel('Final Position')
pl.ylabel('Frequency')
pl.title('1-D Random Walk Distance')
pl.grid(True)
print("Average Position", np.mean(positions))
print("Avg Separation  ", np.mean(abs(positions)))
print("Expectation      %g" % np.sqrt(4*Ns/tau))
Explanation: VECTORIZATION of Implementation
NOTE the difference in arguments!
Inclusive right limit
random.randint(1, 2)
Exclusive right limit
np.random.randint(1, 3)
End of explanation
f = open('afile')
f = open('afile', 'w')
type(f)
f.write('testing\n')
f.close()
cat afile
f = open('afile')  # default mode: read
f.read()
f.write('more tests\n')
f.close()
f = open('afile', 'w')
f.read()
f.close()
f = open('afile', 'a')  # append!
cat afile
Explanation: Files
Basic file operation works via the open() function:
End of explanation
with open('afile', 'w') as f:
    f.write('testing\nmore tests\n')
cat afile
with open('afile', 'r') as f:
    data = f.read()
data
with open('afile', 'r') as f:
    data = f.readlines()
data
s = 'mystring '
s.strip('m ')
data = [item.strip() for item in data]
data
row = ' '.join(['1','2','3'])
f.write(row+'\n')
Explanation: Use so called context manager:
End of explanation
<END_TASK>