Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/BOX COUNTING/BoxPerfectSphere.py
86c9f01a7d54master

Last Change
Raw File

BoxPerfectSphere.py
View Options

import numpy as num
import matplotlib.pyplot as plt
import math
import time
import csv
from mpl_toolkits.mplot3d import Axes3D
import random
from random import randint
#The code calculate the number of boxes intersected for chosen scaling factors for a perfect sphere
random.seed(10)
m = 1000; # number of points of the object (sphere in this case)
d = num.ones(m)
lat = num.zeros(m); # latitude in degrees, initialization
lon = num.zeros(m); # longitude in degrees, initialization
for i in range(m):
lat[i] = randint(-90, 90) # generate random latitudes and longitudes angles (theta and phi)
lon[i] = randint(-180, 180)
xc = num.zeros(m); # initialization of the cartesian coordinates
yc = num.zeros(m);
zc = num.zeros(m);
for i in range(m): # conversion from spherical coordinates into cartesian coo.
xc[i] = d[i] * math.sin(math.radians(lon[i])) * math.cos(math.radians(lat[i]))
yc[i] = d[i] * math.sin(math.radians(lon[i])) * math.sin(math.radians(lat[i]))
zc[i] = d[i] * math.cos(math.radians(lon[i]))
XBoxmin=1.0000001*min(min(xc), min(yc), min(zc)); #Min. coordinate of the first box
XBoxmax=1.0000001*max(max(xc), max(yc), max(zc)); #Max. coordinate of the first box
print('number of points:')
print(m)
print('Coordinates min. and max. of the first box:')
print(XBoxmin)
print(XBoxmax)
LBox = XBoxmax - XBoxmin; #Side length of the first box
print('Length of the side of the first box:')
print(LBox)
steps=18; #Number of steps for the computation
#sf=num.zeros(steps); #scaling factors
sf=[5,7,10,15,20,25,30,40,50,60,70,80,100,120,150,200,250,300];
NBox=num.zeros(steps,dtype=int); #nr of total boxes for each scaling factor
for i in range(steps):
NBox[i]=int(sf[i]**3);
print('scaling factors:')
print(sf)
print('Nr. of total boxes:')
print(NBox)
NN=num.zeros(steps,dtype=int);
t = time.time()
for n in range(0,steps): #ciclo scaling factor
lbox=LBox/sf[n]
array=num.zeros(NBox[n],dtype=int)
for i in range(0,m):
ii=int( (xc[i]-XBoxmin)/lbox)
ij=int( (yc[i]-XBoxmin)/lbox)
ik=int( (zc[i]-XBoxmin)/lbox)
ind=int(ii+sf[n]*ij+sf[n]*sf[n]*ik)
array[ind]=array[ind]+1
for s in range(0,NBox[n]):
if array[s]>0:
NN[n]=NN[n]+1
print('Nr of boxes intersected:')
print(NN)
elapsed = time.time() - t
print('time elapsed [s]:')
print(elapsed)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xc, yc, zc, c='r', marker='o')
plt.show()
c4science · Help