Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/SPHARM/SPHARMCode.py
e6c5d1a9727emaster

Last Change
Raw File

SPHARMCode.py
View Options

from matplotlib import pyplot as plt
import numpy as num
import math
import pyshtools as pysh
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import csv
import time
#Compute the PSD from input file which contains x,y,z coo. in each column
t = time.time()
with open('inputfile.txt','r') as f:
lines = list(csv.reader(f, delimiter = ' ', skipinitialspace = True))
x = num.zeros(len(lines));
y = num.zeros(len(lines));
z = num.zeros(len(lines));
for i in range(len(lines)):
x[i] = lines[i][0];
y[i] = lines[i][1];
z[i] = lines[i][2];
dim=len(lines) #number of points of the object
print(dim)
#Next lines are not necessary if file imported is already centered
#calculation of the centroid of the object
Cx=num.mean(x)
Cy=num.mean(y)
Cz=num.mean(z)
#Translation of the object to have the centroid at the origin (Cx = Cy = Cz = 0)
for i in range(dim):
xc[i] = x[i] - Cx
yc[i] = y[i] - Cy
zc[i] = z[i] - Cz
d=num.zeros(dim); #distance from centroid
lat=num.zeros(dim); #latitude
lon=num.zeros(dim); #longitude
#conversion in spherical coordinates
for i in range(dim):
d[i] = math.sqrt(xc[i]**2 + yc[i]**2 + zc[i]**2)
lat[i] = math.degrees(math.acos(zc[i]/d[i]))-90 #SHTools package take -90° < lat < 90°
lon[i] = math.degrees(math.atan2(yc[i], xc[i])) #and 0 < lon < 360°
norm=1; #Geodesy 4-pi normalized harmonics
csphase=1; #do not apply the Condon-Shortley phase factor to the associated Legendre functions
lmax=80; #Maximum number of the spherical harmonic degree l
#Compute the spherical harmonic coefficients cilm
cilm, chi2 = pysh.expand.SHExpandLSQ(d, lat, lon, lmax, [norm, csphase]);
#Calculate the power spectrum
P = pysh.spectralanalysis.spectrum(cilm);
degrees = num.arange(cilm.shape[1]);
elapsed = time.time() - t
print('time elapsed [s]:')
print(elapsed)
'''
To save P in a text file
np.savetxt('A.txt', A, delimiter=' ', fmt='%s')
'''
#Plot power spectrum
fig, ax = plt.subplots(1, 1)
plt.loglog(degrees, P, label='PSD from centroid')
ax.set(xlabel='Spherical harmonic degree', ylabel='Power')
ax.grid()
plt.legend(loc="upper right")
#Plot points of the object
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, c='r', marker='o')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
c4science · Help