Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/SPHARM/SpharmSpehereEllipsoidGauss.py
e6c5d1a9727emaster

Last Change
Raw File

SpharmSpehereEllipsoidGauss.py
View Options

from matplotlib import pyplot as plt
import numpy as num
import math
import pyshtools as pysh
import random
from random import randint
import time
import csv
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
#This code compute the spherical harmonics coefficients of a sphere of radius 1
t = time.time()
random.seed(10)
a = 1; #parameters of ellipsoid
b = 0.3;
m = 30000; #number of points of the object (sphere in this case)
d = num.ones(m); #value of the function at lat and lon (always equal to 1)
lat = num.zeros(m); #latitude in degrees, initialization
lon = num.zeros(m); #longitude in degrees, initialization
dg = num.zeros(m);#value of the fucntion at lat and lon for gaussian sphere
lmax = 10; #maximum number of spherical harmonic degree l
norm = 1; #Geodesy 4-pi normalized harmonics
csphase = 1; #do not apply the Condon-Shortley phase factor to the associated Legendre functions
#Gaussian noise
noise = num.random.normal(0, 0.05, m);
print(min(noise))
print(max(noise))
xE=num.zeros(m); #initialization of the cartesian coordinates of perfect ellipsoid
yE=num.zeros(m);
zE=num.zeros(m);
dE=num.zeros(m); #distance from centroid of perfect ellipsoid
dEg=num.zeros(m); #distance from centroid of ellipsoid with gaussian noise
dEd=num.zeros(m); #difference of distances between perfect ellipsoid and ellipsoid with noise
for i in range (m):
lat[i]=randint(-90,90); #generate random latitudes and longitudes angles (theta and phi)
lon[i]=randint(0,360);
for i in range (m): #conversion from spherical coordinates into cartesian coo.
xE[i] = a * math.sin(math.radians(lon[i])) * math.cos(math.radians(lat[i]))
yE[i] = b * math.sin(math.radians(lon[i])) * math.sin(math.radians(lat[i]))
zE[i] = math.cos(math.radians(lon[i]))
dE[i] = math.sqrt(xE[i]**2 + yE[i]**2 + zE[i]**2)
for i in range (m):
dg[i] = d[i] + noise[i]
dEg[i] = dE[i] + noise[i]
dEd[i]=dEg[i]-dE[i]
#Compute the spherical harmonic coefficients cilm
cilmg, chi2 = pysh.expand.SHExpandLSQ(dg, lat, lon, lmax, norm, csphase) #sphere with noise
cilmE, chi2 = pysh.expand.SHExpandLSQ(dE, lat, lon, lmax, norm, csphase) #perfect ellipsoid
cilmEg, chi2 = pysh.expand.SHExpandLSQ(dEg, lat, lon, lmax, norm, csphase) #ellipsoid with noise
cilmEd, chi2 = pysh.expand.SHExpandLSQ(dEd, lat, lon, lmax, norm, csphase) #difference ellipsoid
cilmdiff=cilmEg-cilmE; #difference ellispoid considering the spherical harmonic coefficients
#Calculate the power spectrum
Pg = pysh.spectralanalysis.spectrum(cilmg);
degreesg = num.arange(cilmg.shape[1]);
PE = pysh.spectralanalysis.spectrum(cilmE);
degreesE = num.arange(cilmE.shape[1]);
PEg = pysh.spectralanalysis.spectrum(cilmEg);
degreesEg = num.arange(cilmEg.shape[1]);
PEd = pysh.spectralanalysis.spectrum(cilmEd);
degreesEd = num.arange(cilmEd.shape[1]);
Pdiff = pysh.spectralanalysis.spectrum(cilmdiff);
degreesdiff = num.arange(cilmdiff.shape[1]);
#Plot power spectrum
fig, ax = plt.subplots(1, 1)
plt.loglog(degreesg, Pg,'s--',color='green',label='Sphere Gaussian noise')
plt.loglog(degreesE, PE,color='darkorange',label='Ellipsoid')
plt.loglog(degreesEg, PEg,'--',color='red',label='Ellipsoid Gaussian noise')
plt.loglog(degreesEd, PEd,color='blue',label='Ellipsoid difference distance')
plt.loglog(degreesdiff, Pdiff,'--',color='deepskyblue',label='Ellipsoid difference coefficients')
ax.set(xlabel='Spherical harmonic degree', ylabel='Power')
ax.grid()
plt.legend(loc="best")
#Export Powers
num.savetxt('Pg.txt', Pg, delimiter=' ', fmt='%s')
num.savetxt('PE.txt', PE, delimiter=' ', fmt='%s')
num.savetxt('PEg.txt', PEg, delimiter=' ', fmt='%s')
num.savetxt('PEd.txt', PEd, delimiter=' ', fmt='%s')
num.savetxt('Pdiff.txt', Pdiff, delimiter=' ', fmt='%s')
#Plot the random generated points belonging to the ellipsoid
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xE, yE, zE, c='r', marker='o')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
c4science · Help