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

def intensitymaker(phases):
    #given an array of phases, returns electrical intensity
    summedphase = 0
    for i in range(500):
        for k in range(500):
            sumone = math.cos(phases[i] - phases[k]) #maybe should not be cosine?
        summedphase = sumone +summedphase
    constants = 1 #needs calculation. set as one
    intensity = constants * summedphase
    return(intensity)
    
def timestep(x_pos,y_pos,z_pos):
    #returns x y and z positions after a time step
    sd = 1 #fake standard deviation for now
    for i in range(len(x_pos)):
        x_pos[i] = x_pos[i] + np.random.normal(0,sd) #need to calculate standard deviation here
        y_pos[i] = y_pos[i] + np.random.normal(0,sd)
        z_pos[i] = z_pos[i] + np.random.normal(0,sd)
    return(x_pos,y_pos,z_pos)

def phasemaker(x_pos,y_pos,z_pos):
    #given positions of particles relative to their starting ppositions returns phases
    #hold on, I do not think I need z_pos? Because phase is dot product of vector K and r, and the z component ok K is 0
    #so K is K = kx + ky, where k is lambda/2pi. K dot r is kx*x +ky*y +0*z
    k = 1 #set as one for now
    phase = []
    for i in range(len(x_pos)):
        phase_one = x_pos[i] + y_pos[i]
        phase.append(phase_one)
    return(phase)

x_position = np.zeros(1000)
y_position = x_position
z_position = x_position

intensity_time = np.array([])
time = np.array([])

for i in range(500):

    time = np.append(time, int(i))
    anything = phasemaker(x_position, y_position, z_position)
    value = intensitymaker(anything)
    intensity_time = np.append(intensity_time, value)
    x_position, y_position, z_position = timestep( x_position, y_position, z_position)

intensity_time[0] = 0    
#x = np.correlate(intensity_time, intensity_time)
#plt.acorr(intensity_time)
plt.plot(time, intensity_time)
plt.show()

import numpy as np
x = np.zeros(100)
print(x)