# FlockingAnalyzeDivergence.py 
# Run as FlockingAnalyzeDiverge.py filename timestep

# Open two data files of flocks with slightly different initial conditions created by FlockingSimulatorDivergence2.py.
# Now consider the vector of dimension 2*N of the difference in positions between the two flock at each time step, delta_P.
# We plot delta_P as a function of a time and compute a single liapunov exponent for the two flocks based on delta_P.
# (Velocities are ignored because they are too irratic). 

# determine filename and timestep
import sys
filename = sys.argv[1]
filenamecopy = filename + 'copy'
timestep = float(sys.argv[2])

# import modules
from numpy import *
from pylab import *
from FlockingGetPositions import *

# Get the position data for the Original Flock
Data = GetPositions(filename)
flock_positions_x = Data[0]
flock_positions_y = Data[1]
bird_positions_x = Data[2]
bird_positions_y = Data[3]
N = Data[4]
num_timesteps = Data[5]

# Get the position data for the 'Perturbed' Flock
DataCopy = GetPositions(filenamecopy)
flock_positions_x_copy = DataCopy[0]
flock_positions_y_copy = DataCopy[1]
bird_positions_x_copy = DataCopy[2]
bird_positions_y_copy = DataCopy[3]

# Determine the total 'position-distance' between the original and modified flocks
Displacements = []
LogDisplacements = []
for t in xrange(num_timesteps):
    D = 0.0
 
    for n in xrange(N):
        D = D + (bird_positions_x[n][t] - bird_positions_x_copy[n][t])**2
        D = D + (bird_positions_y[n][t] - bird_positions_y_copy[n][t])**2
 
    D = D**0.5
    Displacements.append(D)
    LogDisplacements.append(log(D)) 



# Plot the Displacement Data as a Function of Time
T = timestep*arange(num_timesteps) 

figure(1,(8,6))
title('Net Diplacements as a funciton of time')
xlabel('time')
ylabel('Dislacement')
plot(T,Displacements)
savefig(filename + 'Displacements', dpi=600)


# Plot the LogDisplacement Data as a Function of Time
T = timestep*arange(num_timesteps) 

figure(2,(8,6))
title('Log(Net Diplacements) as a funciton of time')
xlabel('time')
ylabel('Log(Dislacement)')
plot(T,LogDisplacements)
savefig(filename + 'LogDisplacements', dpi=600)


# Calculate the liapunov exponent as a Function of time and plot the time series
Initial_Displacement = Displacements[0]
Liapunov_Exponents = []
for t in xrange(num_timesteps):
    if t == 0: 
       Liapunov_Exponents.append(0)
    if t > 0:  
       Liapunov_Exponents.append( (1/(t*timestep)) * log(Displacements[t]/Initial_Displacement)  )

figure(3,(8,6))
title('Liapunov Exponents as a function of time')
xlabel('time')
ylabel('lambda')

plot(T,Liapunov_Exponents)
savefig(filename + 'LiapunovExponents', dpi=600)

show()