#!/usr/bin/python
# -*- coding: utf-8 -*-

# 2015-06-09 MW
# program to exclude the contacts on the boundary of the DEM assembly
#   for reasons of comparison with the images from DEM2Image
#   the conctacts on the boundary (periodic bc) will not appear in the images...

import numpy
import os, sys, time

# input
folder = "cubic_loose_cycles/"
filenames = [   "triax.-1.eps=0.spheres_image_gaussian=08_noise=03",
                "triax.0.eps=-0.0025.spheres_image_gaussian=08_noise=03",
                "triax.1.eps=-0.005.spheres_image_gaussian=08_noise=03",
                "triax.2.eps=-0.01.spheres_image_gaussian=08_noise=03",
                "triax.3.eps=-0.02.spheres_image_gaussian=08_noise=03",
                "triax.4.eps=-0.03.spheres_image_gaussian=08_noise=03",
                "triax.5.eps=-0.04.spheres_image_gaussian=08_noise=03",
                "triax.6.eps=-0.05.spheres_image_gaussian=08_noise=03",
                "triax.7.eps=-0.06.spheres_image_gaussian=08_noise=03",
                "triax.8.eps=-0.07.spheres_image_gaussian=08_noise=03",
                "triax.9.eps=-0.08.spheres_image_gaussian=08_noise=03",
                "triax.10.eps=-0.09.spheres_image_gaussian=08_noise=03",
                "triax.11.eps=-0.1.spheres_image_gaussian=08_noise=03",
                "triax.12.eps=-50124.6.spheres_image_gaussian=08_noise=03",
                      ]              

NBRUN = len(filenames)


def excluding_contacts(filename):
    filename_positions    = folder + "positions/"+filename+".spheres"
    filename_contacts     = folder + "fabric/"+filename+".contacts"
    num_lines = sum(1 for line in open(filename_positions+".txt"))


    #print "\t ------------------------------------------"
    #print "\t Program to turn exclude contacts on the periodic boundary"
    #print "\t ------------------------------------------\n"


    #print "\t loading input files ...",
    # loading the dem file of particles
    box_size_dem = numpy.genfromtxt(filename_positions+".txt",skip_footer=(num_lines-1),comments='%',usecols=(1,2,3))
    particles = numpy.loadtxt(filename_positions+".txt",skiprows=0,usecols=(0,2,3,4,5))

    # loading the dem file of contacts
    contacts = numpy.loadtxt(filename_contacts+".txt")
    #print "done"

    contacts_new = numpy.copy(contacts)

    # check contact based
    for i in range(0,len(contacts)):
      position_1 = particles[int(contacts[i,0]),1:4]
      position_2 = particles[int(contacts[i,1]),1:4]
      radius_1 = particles[int(contacts[i,0]),4]
      radius_2 = particles[int(contacts[i,1]),4]
      distance = numpy.linalg.norm(position_1 - position_2)
      radii = radius_1+radius_2
      if (distance > radii):
        #print "ALERT! This contact is on the boundary and will be excluded!"
        index = numpy.where((contacts_new[:,0]==contacts[i,0]) & (contacts_new[:,1]==contacts[i,1]))
        contacts_new = numpy.delete(contacts_new,index,0)

    myfile = open(filename_contacts+'-exclude.txt', 'w')
    myfile.write( "# idA idB nx ny nz Fn uN cellDistX cellDistY cellDistZ \n" )
    for i in range(0,len(contacts_new)):
        myfile.write( "%i %i "%( contacts_new[i,0], contacts_new[i,1] ) )
        myfile.write( "%g %g %g "%( contacts_new[i,2], contacts_new[i,3], contacts_new[i,4] ) )
        myfile.write( "%g %g "%( contacts_new[i,5], contacts_new[i,6] ) )
        myfile.write( "%i %i %i "%( contacts_new[i,7], contacts_new[i,8], contacts_new[i,9] ) )
        myfile.write( "\n" )
        
    myfile.close()
    
if __name__ == '__main__':
  for i in range (0, NBRUN):
    print ''
    print '-------------------------------------'
    print "\nRun Number %03i/%03i"%( i+1, NBRUN )
    excluding_contacts(filenames[i])
