Recursive function call python [closed]

I am a Physics intern at TCD. My summer project is to write a Python code that searches, in a random packing, for "ratthers" i.e. particles that have a number of contacts with other particles at most equal to the degrees of freedom of the given system.

I wrote a code that performs a search for the rattlers and eliminates them from a system of 128 particles, in 2D.

I need to perform this search until all the rattlers are gone, recursively. (If you remove a particle from a system, the neighbours will suffer a decrease of contact number by a factor of 1, so everytime you remove a rattler, chances are that you create some other)

My code is pretty elementary and messy, as I am not familiar with coding nor python.

I hope my question is clear by studying the code: i need to call the function "Final_step" until all rattlers are gone.

Thanks for your help

B-

# -*- coding: utf-8 -*-
"""
Created on Fri Jun 21 13:39:06 2013

@author: blaise
"""

from decimal import Decimal
from numpy import loadtxt, dtype, abs, mean



def Sum_radii(particle1, particle2): #Sum of radii - to be compared with distance
    summ= Decimal(particle1.radius)+Decimal(particle2.radius)
    return summ


def PBCdist(particle1, particle2): #PBC conditions to compute distance for any two     particles
     dx = abs(Decimal(particle1.x)-Decimal(particle2.x))
     dy = abs(Decimal(particle1.y)-Decimal(particle2.y))
     if Decimal(dx) > 0.5:
         dx = Decimal(1.0)-Decimal(dx)
     else: 
         dx = Decimal(dx)
     if Decimal(dy) > 0.5:
         dy = Decimal(1.0)-Decimal(dy)
     else: 
         dy = Decimal(dy)
     DX = Decimal(dx)**Decimal(2)
     DY = Decimal(dy)**Decimal(2)
     dist = Decimal(DX)+Decimal(DY)
     distance = Decimal(dist)**Decimal(0.5)
     return distance



def Rattler_killer(contact_number): #must optimize - appends rattlers to respective     array and removing them from contact_number
    for z in range(len(contact_number)):
         if contact_number[z] <= 2:
             rattlers.append(contact_number[z])
    for j in range(len(rattlers)):
         contact_number.remove(rattlers[j])
    return contact_number


 def Contacts_rattlers_post_mortem(contact_number): #array of indexes of the rattlers with nonzero contacts
     for k in range(len(contact_number)):
         if contact_number[k] <=2 and contact_number[k] >0:
             funny_rattlers.append(k)
     return funny_rattlers


 def Contact_fixer(my_particles): #quite messy, but it should do - reduces neighbour contact number
     funny_particles = Contacts_rattlers_post_mortem(contact_number)
     for k in range(len(funny_particles)):
         particle1 = my_particles[funny_particles[k]]
         for particle2 in my_particles:
         distance = PBCdist(particle1, particle2)
         sum_of_radii = Sum_radii(particle1, particle2)
         if (distance < sum_of_radii) and (distance>0):
              previous = contact_number[particle2.n]
              after = previous - 1
              contact_number[particle2.n] = after


def Rattler_contact_reducer(funny_rattlers): #rounds the contact number per rattler to zero - for the funny ones
    for i in range(len(funny_rattlers)):
        contact_number[funny_rattlers[i]]= 0



def Final_step(contact_number): #combines all actions
    Contact_fixer(my_particles)
    Rattler_killer(contact_number)




class Particle():
    def __init__(self, (x,y), n, radius):
        self.x = x
        self.y = y
        self.radius = radius
        self.n = n


number  = loadtxt("second.txt", usecols=(0,), unpack=True, dtype = int)
c1,c2,r = loadtxt("second.txt", usecols=(1,2,4), unpack=True, dtype=dtype(Decimal))



number_of_particles = len(number)
my_particles        = []
overlap             = []
contact_number      = []
rattlers            = []
funny_rattlers      = []


for i in range(number_of_particles):
    n = number[i]
    x = c1[i]
    y = c2[i]
    radius = r[i]

    particle = Particle((x,y), n, radius)
    my_particles.append(particle)


for particle1 in my_particles:
    count = 0
    for particle2 in my_particles:
         distance = PBCdist(particle1, particle2)
         sum_of_radii = Sum_radii(particle1, particle2)
         if (distance < sum_of_radii) and (distance>0):
             count +=1
             olap = (Decimal(sum_of_radii) - Decimal(distance))
             overlap.append(olap)
    contact_number.append(count)




Final_step(contact_number)

print mean(contact_number)