myprojects-hg-reborn

Fork

#! /usr/bin/env python
import scipy as s
import numpy as n
import pylab as p

import csv



#import absorbe as abso  #compiled fortran code for detecting molecule-to-gas collisions
import absorbe_snap as absnap  #compiled fortran code for detecting molecule-to-gas collisions


#! /usr/bin/env python
import scipy as s
import csv


def read_input(filename):


    param_name_list=[]
    param_value_list=[]



    reader = csv.reader(open(filename, 'r'), delimiter = ',')
    line = 0
    # I changed the two following lines into  comments since I do not have a header I should skip
    
    #IndexHeader = reader.next()
    #line = line+1
    for line in reader:
        #print line
        param_name, param_value = line[0], line[1]

        param_name_list.append(param_name)
        param_value_list.append(float(param_value))



    #print "param_name_list is, ", param_name_list
    #print "param_name_list is, ", param_value_list

    param_arr=s.asarray(param_value_list)

    #print "param_arr is, ", param_arr

    #print "param_name_list[2] is, ", param_name_list[2]

    return param_arr





def evolve_by_time_step(amplitude,t_step,N_part,L,pos_at_t):

#     amplitude=self.amplitude
#     t_step=self.t_step
#     N_part=self.N_part
#     L=self.L

    pos_at_t_plus_delta_t=pos_at_t+s.sqrt(amplitude*t_step)*s.random.normal(0.,1.,(3*N_part))

    pos_at_t_plus_delta_t=s.remainder(pos_at_t_plus_delta_t,L) #always fold into the box!!!

    return pos_at_t_plus_delta_t


def evolution_and_absorption(monomer_pos,number_monomers,amplitude,t_step,\
                             N_part,L, pos_at_t,r1_arr,ini_coll_mat):


#     amplitude=self.amplitude
#     t_step=self.t_step
#     N_part=self.N_part
#     L=self.L
#     monomer_positions=self.monomer_positions
#     number_monomers=self.number_monomers

#     r1_arr=self.r1_arr


    n_iter=int(round(t_fin/t_step))

    #my_mod=int(round(t_fin/t_step))/10

    save_every=n_iter/number_save

    #print "N_iter is, ", n_iter


    #print "r1_arr is, ", r1_arr
    
    for i in xrange(n_iter):

        if ((s.remainder(i,save_every) == 0) & (i>0) ):
        
            print "i is, ", i
        #print "monomer_pos are, ", monomer_pos[:,0]
        #print "r1_arr is, ", r1_arr
        if (i==0):

            results=absnap.absorption_snap_calc(pos_at_t, monomer_pos[:,0], \
                                                monomer_pos[:,1], \
                                                monomer_pos[:,2], r1_arr, i, ini_coll_mat)
        
            results_temp=s.copy(results)

#             p.save("results_temp.dat", s.transpose(results_temp))
#             print "saved results_temp"

        elif (i !=0):

 #           if (i == 1):

#                p.save("results_temp_1.dat", s.transpose(results_temp))
                #print "saved results_temp"

  #          if (i == 15):

 #               p.save("results_temp_15.dat", s.transpose(results_temp))
                #print "saved results_temp"



            results=absnap.absorption_snap_calc(pos_at_t, monomer_pos[:,0], \
                                                monomer_pos[:,1], \
                                                monomer_pos[:,2], r1_arr, i, results_temp)
        
            results_temp=s.copy(results)

   #          if (i == 14):

#                 p.save("results_temp_14_updated.dat", s.transpose(results_temp))

            

        # absnap.absorption_snap_calc(initial_pos, monomer_pos[:,0], \
#                                             monomer_pos[:,1], \
#                                             monomer_pos[:,2], mon_rad, k, test1)

            
            
            

        #print "results is ,", results

        pos_at_t=evolve_by_time_step(amplitude,t_step,N_part,L,pos_at_t)

        #print "pos.max(), pos.min() are now, ", pos_at_t.max(), pos_at_t.min()
        #print "pos.mean(), pos.std() are now, ", pos_at_t.mean(), pos_at_t.std()


        #print "pos_at_t is, ", pos_at_t


    return results

def evolution_and_absorption_intermediate_save(monomer_pos,number_monomers,amplitude,t_step,\
                             N_part,L, pos_at_t,r1_arr,ini_coll_mat, number_save):


    n_iter=int(round(t_fin/t_step))



    save_every=n_iter/number_save

    snap_num=0

    
    for i in xrange(n_iter):

  
        if (i==0):

            results=absnap.absorption_snap_calc(pos_at_t, monomer_pos[:,0], \
                                                monomer_pos[:,1], \
                                                monomer_pos[:,2], r1_arr, i, ini_coll_mat)
        
            results_temp=s.copy(results)




        elif (i !=0):

            

            temp=1 #which means I am now saving temporary files during the evolution



            results=absnap.absorption_snap_calc(pos_at_t, monomer_pos[:,0], \
                                                monomer_pos[:,1], \
                                                monomer_pos[:,2], r1_arr, i, results_temp)
        
            results_temp=s.copy(results)

            if ((s.remainder(i,save_every) == 0) & (i>0) ):

                snap_num=snap_num + 1
        
                print "saving snapshot number , ", snap_num

                #Now save the intermediate results

                abso_temp=absorption_individual_monomers(number_mon,results,time, temp, snap_num )

            
                stat_coll=stat_collisions(results, time, 0)

                prefix="temp_%01d"%snap_num
                
                filename="_collision_statistics.dat"

                filename=prefix+filename
                
                p.save(filename,stat_coll )





                #Now we perform some statistical analysis of the data

                number_molecules_vs_time=count_molecules(results,N_part,stat_coll, t_step)

                filename="_number_molecules_vs_time.dat"

                filename=prefix+filename

                p.save(filename,number_molecules_vs_time)

                p.save("intermediate_mol_position.dat",pos_at_t )
                p.save("intermediate_time.dat",(t_step*i*s.ones(1)) )



        pos_at_t=evolve_by_time_step(amplitude,t_step,N_part,L,pos_at_t)




    return results





def absorption_individual_monomers(number_mon,new_abso,time, temp, snap_num):

    for m in xrange(number_mon):

        abso_mon=select_absorption(new_abso,(m+1))  #careful about the m+1 appearing in several places:
        #it comes from
        # the fortran routine way of labelling the molecules 

        if (temp ==1):

            prefix="temp_%01d"%snap_num
                            
            filename="_abso_%01d"%(m+1)

            filename=prefix+filename

        elif (temp != 1):

            filename="abso_%01d"%(m+1)
            
        filename=filename+".dat"

        p.save(filename,abso_mon )





        stat_coll=stat_collisions(abso_mon, time, 0)

        if (temp == 1):

            prefix="temp_%01d"%snap_num

            filename="_collision_statistics_%01d"%(m+1)

            filename=prefix+filename


        elif (temp !=1):

            filename="collision_statistics_%01d"%(m+1)


        filename=filename+".dat"

        p.save(filename,stat_coll )

        return 0






class initialize_system:
    
    def __init__(self, L,N_part):

        self.L=L
        self.N_part=N_part


    def results_ini(self):
        
        results=-s.ones((2,N_part), dtype=int )

        return results



    def place_molecules(self):

        molecule_pos=s.random.uniform(0.0, L,(N_part,3) )

        molecule_pos=s.reshape(molecule_pos,(N_part,3))   #randomly located particles in a box
        return molecule_pos


    def place_molecules_new(self):

        molecule_pos=s.random.uniform(0.0, L,(N_part*3) )

        #molecule_pos=s.reshape(molecule_pos,(N_part,3))   #randomly located particles in a box
        return molecule_pos



    def single_sphere(self, radius):
        
        if (round(radius/L)!= 0.0):
            print "error in choosing box size"



        sphere_pos=s.ones((1,3))*L/2.

        return sphere_pos

    def two_spheres(self, radius):

        my_ratio=radius/L
        my_ratio=my_ratio.astype("int")

        if (my_ratio.any() != 0.0):
                       
            print "error in choosing box size"



        sphere_pos=s.ones((2,3))*L/2.

        sphere_pos[0,0]=L/2.-radius[0]
        
        sphere_pos[1, :]=sphere_pos[0, :]
        sphere_pos[1, 0]=L/2.+radius[1] #this makes the system more symmetric

        return sphere_pos



         

#I should have a box in 3D and generate a Wiener process for many particles in 3D



class particle_3D:  #class needed to generate a single Brownian trajectory
    def __init__(self, ini_pos):

        self.ini_pos=ini_pos  #a 3x1 array with the initial particle position

    def trajectory_in_box(self, L, t_fin, amplitude, t_step):
        ini_pos=self.ini_pos
        N_steps=int(t_fin/t_step)+1
        rand=s.sqrt(amplitude*t_step)*s.random.normal(0.,1.,(3*N_steps))

        rand=s.reshape(rand, (N_steps,3))
        rand[0,:]=0.

        trajectory=s.zeros((N_steps,3))

        
        #trajectory=s.cumsum(rand,axis=1) +ini_pos

        trajectory[:,0]=s.cumsum(rand[:,0])+ini_pos[0]
        trajectory[:,1]=s.cumsum(rand[:,1])+ini_pos[1]
        trajectory[:,2]=s.cumsum(rand[:,2])+ini_pos[2]

        
        #fold trajectory into the box
        trajectory=s.remainder(trajectory, L)

        return trajectory




#Now a derived class to deal with many trajectories

class many_particles_3D(particle_3D):  #this way all the methods of the
                                       #class to generate a single trajectory are made
                                       #available here

    def __init__(self, N_part, ini_pos):

        self.N_part=N_part
        self.ini_pos=ini_pos

        #self.val= val  #this will contain the molecule trajectories
        
        

    def  generate_N_particles(self,L, t_fin, amplitude, t_step):

        global collect_trajectories
        

        #I iterate the result of using the previous class
        
        ini_pos=self.ini_pos

        N_part=self.N_part
        collect_trajectories=s.zeros(((int(t_fin/t_step)+1), (3*N_part))) #array containing
        #many stochastic trajectories 

        
        for i in xrange(N_part):
            
            my_part=particle_3D(ini_pos[i,:])

            #I now iterate the calls to the class generating a single stochastic trajectory


            my_trajectory=my_part.trajectory_in_box(L, t_fin, amplitude, t_step)

            collect_trajectories[: ,(3*i):(3*i+3)]=my_trajectory[: ,0:3]

        return collect_trajectories


    def  generate_N_particles_no_save(self,L, t_fin, amplitude, t_step):

        #Same as the class above, but the generated stochastic trajectories are NOT
        #stored in a global variable. Hence I call this method only from another method to
        #save memory if I need to do some work with the stochastic trajectories but I do not
        #want to store them in memory forever.
        

        #I iterate the result of using the previous class
        
        ini_pos=self.ini_pos

        N_part=self.N_part
        collect_trajectories=s.zeros(((int(t_fin/t_step)+1), (3*N_part)))
        for i in xrange(N_part):

            print "generating particle trajectory number ", i
            
            my_part=particle_3D(ini_pos[i,:])

            my_trajectory=my_part.trajectory_in_box(L, t_fin, amplitude, t_step)

            collect_trajectories[: ,(3*i):(3*i+3)]=my_trajectory[: ,0:3]

            p.save("trajectory_no_save.dat",collect_trajectories )

        return collect_trajectories

    def mean_on_trajectories(self):

        mean_arr=s.zeros(((int(t_fin/t_step)+1),3))

        sel_arr=s.arange(N_part)
        
        mean_arr[:,0]=s.mean(collect_trajectories[:,(sel_arr*3)], axis=1)

        mean_arr[:,1]=s.mean(collect_trajectories[:,(sel_arr*3+1)], axis=1)

        mean_arr[:,2]=s.mean(collect_trajectories[:,(sel_arr*3+2)], axis=1)

        return mean_arr


    def var_on_trajectories(self):

        var_arr=s.zeros(((int(t_fin/t_step)+1),3))

        sel_arr=s.arange(N_part)

        var_arr[:,0]=s.var(collect_trajectories[:,(sel_arr*3)], axis=1)

        var_arr[:,1]=s.var(collect_trajectories[:,(sel_arr*3+1)], axis=1)

        var_arr[:,2]=s.var(collect_trajectories[:,(sel_arr*3+2)], axis=1)

        return var_arr


    def absorption(self,L, t_fin, amplitude,\
                   monomer_positions,mon_radius, t_step, read, fold):

        x_arr=monomer_pos[:,0]
        y_arr=monomer_pos[:,1]
        z_arr=monomer_pos[:,2]

        if (read==0):
        
            trajectories=self.generate_N_particles_no_save(L, t_fin, amplitude, t_step)
        elif (read==1):
            trajectories=n.loadtxt("random_list")
            if (fold == 1):  #In this case I need to fold 
                trajectories=s.remainder(trajectories,L)

        

        my_results=abso.absorption_calc(trajectories, x_arr, y_arr, z_arr, mon_radius)

        #print "the shape if my_results is, ", s.shape(my_results)

        #print "the min is, ", min(my_results.ravel())

        #my_results[:,0]=my_results[:,0]*t_step
        
        p.save("absorption.dat", my_results)
        p.save("absorption_transposed.dat",s.transpose(my_results))


        return my_results


def stat_collisions(collisions, time, save_time):
    myt_step=time[2]-time[1] #detect the time step directly from the time array
    print "myt_step is, ", myt_step
    collision_times=collisions[0,:]
    sel=s.where(collision_times>0.) #get rid of the collisions taking place at t=0
    #(molecules initially overlapping with monomers) and of the trajectories never leading
    #to collisions (an initially absorbed molecule is detected in the first
    #step (k=1) in the loop, hence I need to get rid of absorption "times"
    #smaller than 1)
#    p.save("time[sel].dat", time[sel] )

    collision_time_sel=collision_times[sel]
    #print "collision_sel is, ", collision_sel
    #print "len(collision_sel) is, ", len(collision_sel)

    #Check  carefully this bit!!!!
    
    #t_sorted=s.sort(time[collision_sel]) #hence read only the meaningful times
    t_sorted=s.sort(collision_time_sel) #hence read only the meaningful times

    t_sorted_uni=n.unique1d(t_sorted)  #in case I have several collisions at the same time

    t_sorted_uni=t_sorted_uni

    if (save_time==1):

        p.save("t_sorted_uni.dat", t_sorted_uni*t_step) #I need to multiply the time by t_step in
    #order to get the time array in the proper units
    #count how many collisions I have each time step
    r_bound=t_sorted.searchsorted(t_sorted_uni,side='right')
    l_bound=t_sorted.searchsorted(t_sorted_uni,side='left')

    coll_per_time_step=r_bound-l_bound

    print "the total number of collisions is, ", coll_per_time_step.sum()
    
    results=s.zeros((len(t_sorted_uni),2))

    results[:,0]=t_sorted_uni*t_step #no need to do anything here, I already have
                              #time
    #again, I had to multiply by t_step above in order to get the time in the right units
    
    results[:,1]=coll_per_time_step 

    return results


#Now I count the number of molecules in the system at each time step

def count_molecules(collisions,N_part,results_on_collisions, t_step):
    collision_times=collisions[0,:]
    #print "collision_times is, ", collision_times

    
    
    sel=s.where(collision_times==0) #get rid of the collisions taking place at t=0

    #from the raw data about the collisions, I get rid of the trajectories
    #starting with immediate absorption

    
    N_ini=N_part-len(sel[0]) #initial number of molecules (without those
    #immediately absorbed)

    print "N_ini is, ", N_ini

    N_save=s.array([N_ini])

    p.save("effective_initial_number_molecules.dat", N_save )
    
    N_part_vs_time=s.ones((len(results_on_collisions)+1))*N_ini #array with the
    #number of molecules vs time (just initialized here)

    coll_per_time_step=results_on_collisions[:,1]
    #array telling how many collisions I have at each time step when
    #one or more collisions are recorded.
    #This array comes from the output of stat_collisions(collisions, time) 

    #print "coll_per_time_step is, ", coll_per_time_step

    count_coll=s.ones(len(results_on_collisions))

    #print "len(results_on_collisions) is, ", len(results_on_collisions)
    
    N_part_vs_time[1:]=N_ini-s.cumsum(coll_per_time_step) #subtract a molecule for
    #every collision event

    #print "len(N_part_vs_time) is, ", len(N_part_vs_time)

    #Now return time and the number of molecules vs time.

    results=s.zeros((len(N_part_vs_time),2))
    results[1:,0]=results_on_collisions[:,0]
    results[:,1]=N_part_vs_time

    return results


def select_absorption(new_abso,sel_monomer): #function to extract info about a specific monomer

    select_monomer_collisions=s.where(new_abso[1,:]==sel_monomer)

    monomer_collisions=s.zeros((2,len(select_monomer_collisions[0])))

    monomer_collisions[0,:]=new_abso[0,select_monomer_collisions]
    monomer_collisions[1,:]=new_abso[1,select_monomer_collisions]

    return monomer_collisions


    
    



# def absorption(monomer_positions,mon_radius, t_step, trajectories ):

#     x_arr=monomer_pos[:,0]
#     y_arr=monomer_pos[:,1]
#     z_arr=monomer_pos[:,2]

#     my_results=abso.absorption_calc(trajectories, x_arr, y_arr, z_arr, mon_radius)

#     return my_results

#######################################################################
#######################################################################
#######################################################################
#######################################################################

# L=3.

# t_fin=20.

# t_step=0.001

# amplitude=2.

# N_part=30000


# number_mon=2  #indicate the number of monomers

# mon_rad=s.ones(number_mon)*0.5



input_arr=read_input("input.csv")

#print "input_arr is, ", input_arr


L=input_arr[0]

#print "L is, ", L

t_ini=input_arr[1]

t_fin=input_arr[2]

t_step=input_arr[3]



N_part=int(round(input_arr[4]))
number_mon=int(round(input_arr[5]))
number_save=int(round(input_arr[6]))  #how many intermediate results I am supposed to save

read_ini_pos=int(round(input_arr[7])) #whether I should read the initial position of the molecules \
# from another simulation or not

read_radia=int(round(input_arr[8])) #whether I should read the list of the radia of the monomers from a file

if (read_radia == 1):
    mon_rad=p.load("radia.csv", delimiter=",")


amplitude=input_arr[9]


#I do not need that any longer since I am now using the code in a different way

#fold=1  #whether I should fold the trajectories inside the box or they have already been saved that way



savetime=1 #whether I should save the collision times expliticitely in a file or now


print "amplitude is, ", s.sqrt(amplitude*t_step)

time=s.linspace(0.,t_fin,(int(t_fin/t_step)+1) )
p.save("time.dat", time)

#I am no longe using the "read" parameter: the code no longer reads the trajectories to start with

#read=1  #it means that the code will read some already-generated trajectories
# rather than calculate them on  the fly



#First I need to initialize the system


initialization=initialize_system(L, N_part) #give box size and number of molecules

initial_pos=initialization.place_molecules() #place the molecules within the box

#Now place a single sphere at the centre of the box

if (number_mon == 1):
    monomer_pos=initialization.single_sphere(mon_rad)

#or two of them, with one in the middle

elif (number_mon == 2):
    monomer_pos=initialization.two_spheres(mon_rad)



print "monomer_pos is ", monomer_pos


# my_particles=many_particles_3D(N_part,initial_pos)



# new_abso=my_particles.absorption(L, t_fin, amplitude,\
#                    monomer_pos,mon_rad, t_step, read, fold )


#I replace the bits above with the new evolution of the system


coll_matrix_ini=initialization.results_ini()

test1=s.copy(coll_matrix_ini)

print "coll_matrix_ini is, ", coll_matrix_ini

print "OK here"

#evolution_contained=new_evolution( L,N_part, monomer_pos, \
 #                                  amplitude, t_step, t_fin,mon_rad, number_mon)
print "OK here 2"

initial_pos=initialization.place_molecules_new() #place the molecules within the box

print "initial_pos.mean() is, ", initial_pos.mean()

print "initial_pos.max() and min are, ", initial_pos.max(), initial_pos.min()


# print "initial_pos is, ", initial_pos

# initial_pos=evolve_by_time_step(amplitude,t_step,N_part,L,initial_pos)

# print "initial_pos is, ", initial_pos

# k=0

# print "test1 is, ", test1

# test1=absnap.absorption_snap_calc(initial_pos, monomer_pos[:,0], \
#                                             monomer_pos[:,1], \
#                                             monomer_pos[:,2], mon_rad, k, test1)


# print "test1 is, ", test1

# print "OK after test1"



#new_abso=evolution_and_absorption(monomer_pos,number_mon,amplitude,t_step,\
#                            N_part,L, initial_pos,mon_rad ,coll_matrix_ini)


new_abso=evolution_and_absorption_intermediate_save(monomer_pos,number_mon,amplitude,t_step,\
                             N_part,L, initial_pos,mon_rad,coll_matrix_ini, number_save)


print "OK after test2"


p.save("new_abso_test.dat", s.transpose(new_abso))

# new_abso=evolution_contained.evolution_and_absorption(initial_pos,coll_matrix_ini)

#####################################################################
#####################################################################
#####################################################################
#####################################################################
#####################################################################

# The following part does not need any changing, since it deals with the absorption
#statistics



#Now carry out an iteration on all the monomers  


for m in xrange(number_mon):

    abso_mon=select_absorption(new_abso,(m+1))  #careful about the m+1 appearing in several places: it comes from
    # the fortran routine way of labelling the molecules 

    filename="abso_%01d"%(m+1)

    filename=filename+".dat"

    p.save(filename,abso_mon )





    stat_coll=stat_collisions(abso_mon, time, 0)

    filename="collision_statistics_%01d"%(m+1)

    filename=filename+".dat"

    p.save(filename,stat_coll )


 





print "new_abso is ok!"

#p.save("absorption_data.dat", new_abso)

stat_coll=stat_collisions(new_abso, time, savetime)

p.save("collision_statistics.dat",stat_coll )





#Now we perform some statistical analysis of the data

number_molecules_vs_time=count_molecules(new_abso,N_part,stat_coll, t_step)

p.save("number_molecules_vs_time.dat",number_molecules_vs_time)


fig = p.figure()
axes = fig.gca()

axes.plot(number_molecules_vs_time[:,0],number_molecules_vs_time[:,1],"ro",linewidth=2.)
# axes.plot(time,(mean_stat[:,1]-L/2.),"k+",\
#           label="<y>")
# axes.plot(time,(mean_stat[:,2]-L/2.),"mx",\
#           label="<z> ",linewidth=2.)


p.xlabel('Time')
p.ylabel('Number of Molecules')
p.savefig("number_molecules_vs_time.pdf")
 
p.clf()    





fig = p.figure()
axes = fig.gca()

axes.plot(stat_coll[:,0],stat_coll[:,1],"ro",linewidth=2.)
# axes.plot(time,(mean_stat[:,1]-L/2.),"k+",\
#           label="<y>")
# axes.plot(time,(mean_stat[:,2]-L/2.),"mx",\
#           label="<z> ",linewidth=2.)


p.xlabel('Time')
p.ylabel('Collisions')
p.savefig("collisions_vs_time.pdf")
 
p.clf()    






#The following are just some tests which I commented out.


# initial_pos=s.ones((N_part,3))*L/2.

# my_particles=many_particles_3D(N_part,initial_pos)

# my_simu_trajectories=my_particles.generate_N_particles(L, t_fin, amplitude, t_step)

# print "the trajectories are, ", my_simu_trajectories

# p.save("trajectories.dat",my_simu_trajectories )


# abso_results=absorption(monomer_pos,mon_rad, t_step, my_simu_trajectories)



# print "abso_results is, ", abso_results

# p.save("absorption_data.dat", s.abso_results)

# mean_stat=my_particles.mean_on_trajectories()

# print "mean_stat is, ", mean_stat


# var_stat=my_particles.var_on_trajectories()

# print "var_stat is, ", var_stat

# print "collect_trajectories is, ", collect_trajectories



# fig = p.figure()
# axes = fig.gca()

# axes.plot(time,(mean_stat[:,0]-L/2.),"ro",\
#           label="<x>",linewidth=2.)
# axes.plot(time,(mean_stat[:,1]-L/2.),"k+",\
#           label="<y>")
# axes.plot(time,(mean_stat[:,2]-L/2.),"mx",\
#           label="<z> ",linewidth=2.)


# p.xlabel('Time')
# p.ylabel('Mean position')
# p.savefig("mean_x_y_z.pdf")
 
# p.clf()    



# fig = p.figure()
# axes = fig.gca()

# axes.plot(time,var_stat[:,0],"ro",\
#           label="<delta x sq>",linewidth=2.)
# axes.plot(time,var_stat[:,1],"k+",\
#           label="<delta y sq>")
# axes.plot(time,var_stat[:,2],"mx",\
#           label="<delta z sq> ",linewidth=2.)


# p.xlabel('Time')
# p.ylabel('variance')
# p.savefig("var_x_y_z.pdf")
 
# p.clf()    


# fig = p.figure()
# axes = fig.gca()

# axes.plot(time,(collect_trajectories[:,0]-L/2.),"ro",\
#           label=" x",linewidth=2.)
# axes.plot(time,(collect_trajectories[:,1]-L/2.),"k+",\
#           label="y")
# axes.plot(time,(collect_trajectories[:,2]-L/2.),"mx",\
#           label="z ",linewidth=2.)


# p.xlabel('Time')
# p.ylabel('variance')
# p.savefig("individual_trajectories.pdf")
 
# p.clf()    






print "So far so good"