#!usr/bin/python

import os
import itertools as i
import numpy as np
import matplotlib.pyplot as plt
import pylab
from numpy import *

# Declaration of constants

K = 3.14159265359
N= 5
p_za=[]
pv_za_temp=[]
pv_za_temp=np.array(pv_za_temp)
#-----------------------
#Allotment of particles
#-----------------------

p_initial = np.linspace(0,2,num=5)

pv_za_temp = (np.array(p_initial))

#Displacement of particles using Zeldovich Approximation

def t_range(start, end, step):
  while start <= end:
    yield start
    start += step

for t in t_range(0,1,0.1):
  print t
  p_za=[]
  pv_za=[]
  
  # Opening file in file_t format
  fname = 'file_' + str(t) + '.dat'
  fo = open(fname,'w')
  
  # p_za.append(p_initial - t*K*np.sin(K*p_initial))
  print 'K=',K
  print 'pv_za_temp =',pv_za_temp
  print '- t*K*np.sin(K*p_initial) = ',- t*K*np.sin(K*p_initial)
  print '-K*np.sin(K*pv_za_temp) = ',-K*np.sin(K*pv_za_temp)
  pv_za.append(-K*np.sin(K*pv_za_temp))
  pv_za_temp = []
  pv_za_temp.append(np.array(pv_za))
      
 
#Imposing Periodic Boundary condition on p_za
  """if p_za < 0:
    p_za = (2-abs(np.array(p_za)) % 2)
      
  if p_za > 2:
    p_za = np.array(p_za) % 2"""

# Periodic Boundary Condition on pv_za

  if pv_za < 0:
    pv_za = (2-abs(np.array(pv_za)) % 2)
      
  if pv_za > 2:
    pv_za = np.array(pv_za) % 2
      
    
      
  
#Calculation of Density


  """temp_za = np.array(p_za/0.02,dtype=int)
  print temp_za

  for j in range(0,100):
    for i in range(0,N):
      counter= (temp_za==j).sum()        
    fo.write('{0:f}  {1:d} \n'.format(j/25.0, counter)) 
  fo.close()""" 

# Preparing Files for Phase Plot
  for j in range(0,N):

 
    fo.write('{0:f}  {1:f} \n'.format(*i.chain(np.array(p_initial),np.array(pv_za)))) 
  fo.close()
 

# Generating Images 
            
for i in t_range(0,1,0.1):
  fname = 'file_' + str(i) + '.dat'
  f=open(fname,"r")
  x,y = np.loadtxt(fname, unpack=True)
  plt.plot(x,y,linewidth=2.0)
  plt.ylim(0,25)
  plt.ylabel('Rho')
  plt.xlabel('X')
  pylab.savefig(fname + '.png')
  plt.cla()
  f.close()