Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/SPHARM/SpharmCubeGauss.py
e6c5d1a9727emaster

Last Change
Raw File

SpharmCubeGauss.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
import random
from random import randint
#Calculate the PSD of a cube with gaussian noise
mm=10000; #number of points per cube side
Maxcoo=1;
Mincoo=-1;
#gaussian noise (different seed to have non identical cube faces)
random.seed(1)
noise1 = num.random.normal(0, 0.05, mm);
random.seed(2)
noise2 = num.random.normal(0, 0.05, mm);
random.seed(3)
noise3 = num.random.normal(0, 0.05, mm);
random.seed(4)
noise4 = num.random.normal(0, 0.05, mm);
random.seed(5)
noise5 = num.random.normal(0, 0.05, mm);
random.seed(6)
noise6 = num.random.normal(0, 0.05, mm);
#####################
#Perfect Cube
#coordinates of points on face nr. 1
x1=num.zeros(mm);
y1=num.zeros(mm);
z1=num.zeros(mm);
for i in range(mm):
x1[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y1[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z1[i] = Mincoo
#coordinates of points on face nr. 2
x2=num.zeros(mm);
y2=num.zeros(mm);
z2=num.zeros(mm);
for i in range(mm):
x2[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y2[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z2[i] = Maxcoo
#coordinates of points on face nr. 3
x3=num.zeros(mm);
y3=num.zeros(mm);
z3=num.zeros(mm);
for i in range(mm):
x3[i] = Mincoo
y3[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z3[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x4=num.zeros(mm);
y4=num.zeros(mm);
z4=num.zeros(mm);
for i in range(mm):
x4[i] = Maxcoo
y4[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z4[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x5=num.zeros(mm);
y5=num.zeros(mm);
z5=num.zeros(mm);
for i in range(mm):
x5[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y5[i] = Mincoo
z5[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x6=num.zeros(mm);
y6=num.zeros(mm);
z6=num.zeros(mm);
for i in range(mm):
x6[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y6[i] = Maxcoo
z6[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
xc=[];
yc=[];
zc=[];
xc = num.append(xc, x1)
xc = num.append(xc, x2)
xc = num.append(xc, x3)
xc = num.append(xc, x4)
xc = num.append(xc, x5)
xc = num.append(xc, x6)
yc = num.append(yc, y1)
yc = num.append(yc, y2)
yc = num.append(yc, y3)
yc = num.append(yc, y4)
yc = num.append(yc, y5)
yc = num.append(yc, y6)
zc = num.append(zc, z1)
zc = num.append(zc, z2)
zc = num.append(zc, z3)
zc = num.append(zc, z4)
zc = num.append(zc, z5)
zc = num.append(zc, z6)
#################
#Cube with gaussian noise
#coordinates of points on face nr. 1
x1n=num.zeros(mm);
y1n=num.zeros(mm);
z1n=num.zeros(mm);
for i in range(mm):
x1n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y1n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z1n[i] = Mincoo+noise1[i];
#coordinates of points on face nr. 2
x2n=num.zeros(mm);
y2n=num.zeros(mm);
z2n=num.zeros(mm);
for i in range(mm):
x2n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y2n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z2n[i] = Maxcoo+noise2[i];
#coordinates of points on face nr. 3
x3n=num.zeros(mm);
y3n=num.zeros(mm);
z3n=num.zeros(mm);
for i in range(mm):
x3n[i] = Mincoo+noise3[i]
y3n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z3n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x4n=num.zeros(mm);
y4n=num.zeros(mm);
z4n=num.zeros(mm);
for i in range(mm):
x4n[i] = Maxcoo+noise4[i]
y4n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z4n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x5n=num.zeros(mm);
y5n=num.zeros(mm);
z5n=num.zeros(mm);
for i in range(mm):
x5n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y5n[i] = Mincoo+noise5[i]
z5n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x6n=num.zeros(mm);
y6n=num.zeros(mm);
z6n=num.zeros(mm);
for i in range(mm):
x6n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y6n[i] = Maxcoo+noise6[i]
z6n[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
xcn=[];
ycn=[];
zcn=[];
xcn = num.append(xcn, x1n)
xcn = num.append(xcn, x2n)
xcn = num.append(xcn, x3n)
xcn = num.append(xcn, x4n)
xcn = num.append(xcn, x5n)
xcn = num.append(xcn, x6n)
ycn = num.append(ycn, y1n)
ycn = num.append(ycn, y2n)
ycn = num.append(ycn, y3n)
ycn = num.append(ycn, y4n)
ycn = num.append(ycn, y5n)
ycn = num.append(ycn, y6n)
zcn = num.append(zcn, z1n)
zcn = num.append(zcn, z2n)
zcn = num.append(zcn, z3n)
zcn = num.append(zcn, z4n)
zcn = num.append(zcn, z5n)
zcn = num.append(zcn, z6n)
npoints=6*mm;
#To Export coordinates of noisy cube
XCN = num.zeros((npoints,3))
for i in range(npoints):
XCN[i, 0] = xcn[i]
XCN[i, 1] = ycn[i]
XCN[i, 2] = zcn[i]
#num.savetxt('XCN.txt', XCN, delimiter=' ', fmt='%s')
d=num.zeros(npoints); #distance from centroid and surface of perfect cube
dn=num.zeros(npoints); #distance from centroid and surface of noisy cube
lat=num.zeros(npoints); #latitude (perfect cube)
lon=num.zeros(npoints); #longitude (perfect cube)
latn=num.zeros(npoints); #latitude (noisy cube)
lonn=num.zeros(npoints); #longitude (noisy cube)
#conversion in spherical coordinates
for i in range(npoints):
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°
dn[i] = math.sqrt(xcn[i] ** 2 + ycn[i] ** 2 + zcn[i] ** 2)
latn[i] = math.degrees(math.acos(zcn[i] / dn[i])) - 90 # SHTools package take -90° < lat < 90°
lonn[i] = math.degrees(math.atan2(ycn[i], xcn[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=100; #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]); #perfect cube
cilmn, chi2n = pysh.expand.SHExpandLSQ(dn, latn, lonn, lmax, [norm, csphase]); #noisy cube
cilmdiff=cilmn-cilm; #Difference
#Calculate the power spectrum
P = pysh.spectralanalysis.spectrum(cilm);
degrees = num.arange(cilm.shape[1]);
Pn = pysh.spectralanalysis.spectrum(cilmn);
degreesn = num.arange(cilmn.shape[1]);
Pdiff = pysh.spectralanalysis.spectrum(cilmdiff);
degreesdiff = num.arange(cilmdiff.shape[1]);
num.savetxt('P.txt', P, delimiter=' ', fmt='%s')
num.savetxt('Pn.txt', Pn, delimiter=' ', fmt='%s')
num.savetxt('Pdiff.txt', Pdiff, delimiter=' ', fmt='%s')
#Plot power spectrum
fig, ax = plt.subplots(1, 1)
plt.loglog(degrees, P, label='Cube')
plt.loglog(degreesn, Pn, label='Cube noise')
plt.loglog(degreesdiff, Pdiff, label='Difference')
ax.set(xlabel='Spherical harmonic degree', ylabel='Power')
ax.grid()
plt.legend(loc="upper right")
plt.show()
c4science · Help