myprojects-hg-reborn

Fork

#! /usr/bin/env python

import scipy as s
import numpy as n
import pylab as p




size_cluster=5
density=0.03

box_size=(size_cluster/density)**(1./3.)  #here the box size does not depend on the total particle number

pos_arr=p.load("unfolded_final")

#now I have a n array where each row has the 3 coordinate of a particle.


pos_arr=s.remainder(pos_arr,box_size)

p.save("initial_state_folded", pos_arr)

#Now I want to relocate each cluster in such a way that all the clusters (NOT just the centre of mass
# of the system) is relocated at the centre of the new computational box I will use

new_box_size=200.

centre=new_box_size/2. #this is , along x,y, and z, the centre of the new computational box I am using

dim=s.shape(pos_arr)

n_clus=dim[0]/size_cluster

print "the number of clusters is, ", n_clus

ini_pos_clusters=s.zeros((n_clus,3)) #this will contain the positions (in 3D!!!) of the centre of masses of all the clusters
print "pos_arr folded is, ",pos_arr

print "the shape of pos_arr is, ", s.shape(pos_arr)

for i in xrange(n_clus):
    print "i is, ", i
    ini_pos_clusters[i,0]=pos_arr[(i*size_cluster):((i+1)*size_cluster),0].mean()
    ini_pos_clusters[i,1]=pos_arr[(i*size_cluster):((i+1)*size_cluster),1].mean()
    ini_pos_clusters[i,2]=pos_arr[(i*size_cluster):((i+1)*size_cluster),2].mean()


print "the initial positions of the clusters are, ", ini_pos_clusters

#Now I compute the shifts for every cluster

shift_vec=centre-ini_pos_clusters

print "the shift_vec is, ", shift_vec

#Now I actually perform the shift

for i in xrange(n_clus):
    print "i is, ", i
    pos_arr[(i*size_cluster):((i+1)*size_cluster),0]=pos_arr[(i*size_cluster):((i+1)*size_cluster),0]\
                                                      +shift_vec[i,0]
    pos_arr[(i*size_cluster):((i+1)*size_cluster),1]=pos_arr[(i*size_cluster):((i+1)*size_cluster),1]\
                                                      +shift_vec[i,1]

    pos_arr[(i*size_cluster):((i+1)*size_cluster),2]=pos_arr[(i*size_cluster):((i+1)*size_cluster),2]\
                                                      +shift_vec[i,2]

#Now a test to make sure everything is OK

p.save("initial_state_folded_and_shifted", pos_arr)


for i in xrange(n_clus):
    print "i is, ", i
    ini_pos_clusters[i,0]=pos_arr[(i*size_cluster):((i+1)*size_cluster),0].mean()
    ini_pos_clusters[i,1]=pos_arr[(i*size_cluster):((i+1)*size_cluster),1].mean()
    ini_pos_clusters[i,2]=pos_arr[(i*size_cluster):((i+1)*size_cluster),2].mean()


print "the shifted positions of the clusters are, ", ini_pos_clusters


print "So far so good"