myprojects-hg-reborn

Fork

#! /usr/bin/env python

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

import ray_count as rc



class particle_2D:
	#I initialize the position of the particle
    def __init__(self, x0,y0):
        self.x0 = x0
        self.y0 = y0

    def position(self):
        pos=s.array([self.x0, self.y0])

        return pos
        


class circle_2D:
    def __init__(self, x0, y0, R, N):

        self.x0=x0
        self.y0=y0
        self.R=R
        self.N=N


    def gen_theta_arr(self):
        N=self.N
        theta_arr=s.linspace(0.,(2.*s.pi), N) #array with the angles
        return theta_arr


    def  points_on_circle(self):

        N=self.N
        x0=self.x0
        y0=self.y0
        R=self.R
        
        theta_arr=self.gen_theta_arr()

        
        #Now generate the 2D coordinates of every point along the circle in 2D

        x_arr=R*s.cos(theta_arr)
        y_arr=R*s.sin(theta_arr)

        points=s.transpose(s.vstack((x_arr, y_arr)))

        return points

class ray_trace:
    
    def __init__(self, x0, y0, R, N, N_on_ray, points):

        self.x0=x0
        self.y0=y0
        self.R=R
        self.N=N

        self.N_on_ray=N_on_ray
        self.points=points

    

    def gen_theta_arr(self):
        N=self.N
        theta_arr=s.linspace(0.,(2.*s.pi), N) #array with the angles
        return theta_arr


    def rays(self):
        #For each point (i.e. a set of coordinates in 2D)
        #I need to generate a ray made up N_on_ray points pointing towards (x0,y0)


        points=self.points
        x0=self.x0
        y0=self.y0
        N_on_ray=self.N_on_ray
        R=self.R
        N=self.N

        theta_arr=self.gen_theta_arr()
        ray_arr=s.linspace(0.,R, N_on_ray)
        line_arr=s.zeros((N_on_ray, (2*len(points))))

        for i in xrange(len(points)):
            line_arr[:,(2*i)]=ray_arr*s.cos(theta_arr[i])
            line_arr[:,(2*i+1)]=ray_arr*s.sin(theta_arr[i])

        return line_arr
        
        




#Coordinates of the centre of the monomer I want to study

x_0=0.
y_0=0.

r1=1.  #sphere radius



part0=particle_2D(x_0,y_0)
print "part0 position is, ", part0.position()


#Coordinates of the centres of neighboring spheres

x_1=2.
y_1=0.


x_2=-2.
y_2=0.


x_3=0.
y_3=2.


x_4=2.
y_4=0.


R=3.5  #radius of the circle
N=100 #number of pts along the circle (which is also the number or rays!!!)

N_on_ray=50 #  number of pts on each ray

#arrays with the x and y coordinates of the centres of the neighboring spheres

# x_arr=s.array([x_1,x_2,x_3,x_4])
# y_arr=s.array([y_1,y_2,y_3,y_4])


x_arr=s.array([x_1,x_2,x_3])
y_arr=s.array([y_1,y_2,y_3])



my_circle=circle_2D(x_0,y_0,R,N)

points=my_circle.points_on_circle()

print "the points on circle are, ", points 




my_rays=ray_trace(x_0, y_0, R, N, N_on_ray, points)

ray_arr=my_rays.rays()



my_ratio=rc.ray_count(ray_arr, x_arr, y_arr, r1)

print "the shielding factor is, ", my_ratio


#See that I am really covering a circle

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

            
axes.plot(points[:,0],points[:,1],"ro", ray_arr[:,12], ray_arr[:,13], linewidth=2.)
p.xlabel('Rcos(theta)')
p.ylabel('Rsin(theta)')

p.title('Points on circle')
p.grid(True)

p.savefig("points_on_circle.pdf")
             
p.clf()




print "So far so good"