#___________Problem 1__________________________________________________

import numpy as np

'''
   M=50 random 2-D walks of N=100 steps each.
   The length of step in both x and y is equally likely 1 or 2.
   Steps in x and y direction are independent and random.
'''

# this function does one walk of N steps in 2 dimensions
def walk(N):
    x,y = (0,0) # starting from the origin
    for n in range(N):  # N steps
        x += np.random.choice([-2,-1,+1,+2])
        y += np.random.choice([-2,-1,+1,+2])
    return(x*x+y*y)


# main program
# calls function doing one walk of N random walk steps M times
N, M, i, d_aver, d2_aver = 100, 50, 0, 0., 0.
print("\ni      d")
while(i<M):
    i += 1
    d2 = walk(N)
    d = d2**0.5    # final distance of this walk
    d2_aver += d2  # optional computation
    d_aver += d    # compute average distance
    print(i,' ',d) # print 2 columns
print('\n Average distance <d> =',d_aver/M)

# this part is extra

print(' <d2>**0.5 = ',(d2_aver/M)**0.5)
s = np.array([-2,-1,1,2]).std()
sigma2 = N * 2* s**2    
print('expected std deviation of d = sqrt(sigma2) = ',sigma2**0.5,'\n')


#___________Problem 2__________________________________________________
'''
   Curtain is devided into N>>1 segments and its length L is
   summed numerically in a loop.
   It's not the fastest possible way to do it! 
   For efficiency, we could use numpy arrays, realize that the 
   curtain makes 12 identical half-waves between its ends 
   (only one of which needs to be computed) etc. 
   However, this is one of the simplest and clearest algorithms.
   [Plotting at the end was not required.]
'''

def y(x):
    const = 12*np.pi/W
    return (0.4*np.sin(const*x))

# main program for curtain problem
# x changes from 0 to W=6 m.
# y(x) is the sideways deflection of the curtain  
#   
N, W, L = 100000, 6., 0.
dx = W/N
for i in range(N):
    x = i*dx        # x-coordinate of the left end of i0th segment
# y0, y1 are yalues at the begin and end of i-th segment of curtain
    y0, y1 = y(x), y(x+dx)
    dL = ( dx**2 + (y1-y0)**2 )**0.5  # says Pythagoras
    L = L + dL
print('You need to buy',np.around(L,9),'m of fabric for the curtain.')


# optional plotting of curtain
import matplotlib.pyplot as plt
xx = np.linspace(0,W,1000)
plt.plot(xx,y(xx))
plt.plot(xx,0*xx,color=(.8,.8,.8))
plt.plot([0,0],[-2.5,2.5],color=(.9,.9,.9))
plt.plot([6,6],[-2.5,2.5],color=(.9,.9,.9))
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Curtain')
plt.show()