Page MenuHomec4science
Files F87354492

SpharmSpehereEllipsoidGauss.py
No OneTemporary
Actions

File Metadata

Created
Sat, Oct 12, 04:26

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()

Event Timeline

c4science · Help