Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/SPHARM/SpharmSphereNoises.py
e6c5d1a9727emaster

Last Change
Raw File

SpharmSphereNoises.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
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
#Compute the PSD of a perfect sphere and spheres with different degree of noise
#(to show that power with L>0 is not zero when we have noise on a sphere)
m = 10000; #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
lmax = 30; #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
for i in range (m):
lat[i]=randint(-90,90); #generate random latitudes and longitudes angles (theta and phi)
lon[i]=randint(0,360);
#Random Noise
random.seed(10)
d2=num.zeros(m);
d3=num.zeros(m);
d4=num.zeros(m);
for i in range(m):
d2[i] = d[i] - randint(0, 10) / 100 #different degrees of noise
d3[i] = d[i] - randint(0, 10) / 1000
d4[i] = d[i] - randint(0, 10) / 10000
diff4=num.zeros(m);
for i in range(m):
diff4[i]=d[i]-d4[i]
print("LAT")
print(lat)
#Compute the spherical harmonic coefficients cilm
cilm, chi2 = pysh.expand.SHExpandLSQ(d, lat, lon, lmax, norm, csphase) #perfect sphere
cilm2, chi2 = pysh.expand.SHExpandLSQ(d2, lat, lon, lmax, norm, csphase) #spheres with noise
cilm3, chi2 = pysh.expand.SHExpandLSQ(d3, lat, lon, lmax, norm, csphase)
cilm4, chi2 = pysh.expand.SHExpandLSQ(d4, lat, lon, lmax, norm, csphase)
#Calculate the power spectrum
P = pysh.spectralanalysis.spectrum(cilm);
degrees = num.arange(cilm.shape[1]);
P2 = pysh.spectralanalysis.spectrum(cilm2);
degrees2 = num.arange(cilm2.shape[1]);
P3 = pysh.spectralanalysis.spectrum(cilm3);
degrees3 = num.arange(cilm3.shape[1]);
P4 = pysh.spectralanalysis.spectrum(cilm4);
degrees4 = num.arange(cilm4.shape[1]);
#Plot power spectrum
fig, ax = plt.subplots(1, 1)
plt.loglog(degrees, P, label='Perfect sphere')
plt.loglog(degrees2, P2,label='Noise 10^-2')
plt.loglog(degrees3, P3,label='Noise 10^-3')
plt.loglog(degrees4, P4,label='Noise 10^-4')
ax.set(xlabel='Spherical harmonic degree', ylabel='Power')
ax.grid()
plt.legend(loc="best")
plt.show()
c4science · Help