#!/rpod2/jcpassy/yt-x86_64/bin python2.6

import sys
import commands as C
import numpy as N
import pylab as P
import yt.raven as R
import yt.lagos as L

filebase = sys.argv[-1]

minSnap = 0
maxSnap = 401
skip = 1

t = N.zeros(maxSnap/skip)

X = N.zeros((maxSnap/skip,3,2))
V = N.zeros((maxSnap/skip,3,2))
r = N.zeros(maxSnap/skip)


X2 = N.zeros((maxSnap/skip,3,2))
V2 = N.zeros((maxSnap/skip,3,2))
r2 = N.zeros(maxSnap/skip)

j = 0

for i in range(minSnap,maxSnap,skip):
    
    # File 1

    filen = 'DM_part/'+ 'DD' + '%04i' % i + '/' + filebase + '%04i' % i
    print "loading %s" % filen
    data = L.EnzoStaticOutput(filen,data_style="enzo_packed_3d") # if no gas

    t[j] = data["InitialTime"] * data["years"]

    region = data.h.region([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
            
    i1,i2 = region['particle_index']

    X[j,0,i1],X[j,0,i2] = region['particle_position_x']
    X[j,1,i1],X[j,1,i2] = region['particle_position_y']
    X[j,2,i1],X[j,2,i2] = region['particle_position_z']

    V[j,0,i1],V[j,0,i2] = region['particle_velocity_x']
    V[j,1,i1],V[j,1,i2] = region['particle_velocity_y']
    V[j,2,i1],V[j,2,i2] = region['particle_velocity_z']
    
    r[j] = pow((X[j,0,i1] - X[j,0,i2]),2) + pow((X[j,1,i1] - X[j,1,i2]),2) + pow((X[j,2,i1] - X[j,2,i2]),2)
    r[j] = pow(r[j],0.5)    
    
    j = j + 1