Skip to content

Commit

Permalink
improved scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
@maul1609
maul1609 committed Oct 20, 2023
1 parent 6abb9b8 commit 4af76ea
Show file tree
Hide file tree
Showing 4 changed files with 53 additions and 34 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,3 +12,5 @@ doxygen
matlab-dev
fortran_testing
*.png
__pycache__

18 changes: 10 additions & 8 deletions python/analyse_output.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,11 @@
import matplotlib
import sys
matplotlib.use('Agg')
import matplotlib.pyplot as plt

import numpy as np
import fourier_wave_number
import normal_modes_compare

username=getpass.getuser()

Expand All @@ -15,23 +18,22 @@
""" do the fourier analysis
"""

u_jet=[5., 10., 20., 30., 40., 50., 60., 70., 80.,\
90., 100., 125., 150., 175., 200., 250., 300., 350.,];
u_jet=[50., 70., 100.];

c_vis=[0.,0.1,0.2,1.0]
c_vis=[0.2]

plt.ion()
plt.figure()
for j in range(len(c_vis)):
fileNames=['/tmp/' + username + '/output' + str(i) + '_' + str(j) + '.nc' \
fileNames=['/tmp/' + username + '/output_' + str(i) + '_' + str(j) + '.nc' \
for i in range(len(u_jet))]

plt.subplot(1,len(c_vis),j)
plt.subplot(1,len(c_vis),j+1)

fourier_wave_number.do_analysis01(fileNames,u_jet);

fourier_wave_number.do_analysis01(fileNames,u_jet[i]);

# same for each value of c_cvis
normal_modes_compare(u_jet[i],1);
normal_modes_compare.normal_modes_compare(u_jet,1);
plt.title('$C_{vis}$=' + str(c_vis[j]))

plt.savefig('/tmp/' + username + '/full_analysis.png')
Expand Down
36 changes: 23 additions & 13 deletions python/fourier_wave_number.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
import sys
matplotlib.use('Agg')

import matplotlib.pyplot as plt

import numpy as np
from netCDF4 import Dataset as NetCDFFile
import scipy as scy
Expand All @@ -15,7 +17,7 @@
if not os.path.exists('/tmp/' + username):
os.mkdir('/tmp/' + username)

def fourier_wave_number(fileName,u_jet):
def fourier_wave_number(fileName):
"""
uses findpeaks and fourier analysis to track the current wave number and
the rotation speed (i.e. using phase information from the fft).
Expand Down Expand Up @@ -60,7 +62,7 @@ def fourier_wave_number(fileName,u_jet):

ind=len(locs);

phs=np.unwrap(np.angle(Y))*180./np.pi;
phs=np.unwrap(np.angle(Y))*180./np.pi; # phase angle in degrees

print('wave number is : ' + str(f[ind]) + '; phase: ' + str(phs[ind]));
phase[j,i]=phs[ind];
Expand All @@ -69,6 +71,10 @@ def fourier_wave_number(fileName,u_jet):

nc.close();
phaseold[j,:]=phase[j,:];

# phase is the phase of the wave through it's cycle, but
# we want the movement of the wave train around the planet,
# so divide by number of waves.
phase[j,:]=np.unwrap(phaseold[j,:]*np.pi/180.)*180./np.pi/(wave_number[j,:]);
rotation_rate[j,:]=-np.diff(phase[j,:],n=1,axis=0) / \
dt_sec*86400.*365.25/360; # full rotations per year
Expand All @@ -90,20 +96,21 @@ def do_analysis01(fileNames,u_jets):
mean_rot=np.zeros((n_files,9));
std_rot=np.zeros((n_files,9));
for j in range(n_files):
(phase,rotation_rate,wave_number)=fourier_wave_number([fileNames[j]],u_jets)
(phase,rotation_rate,wave_number)=fourier_wave_number([fileNames[j]])
for i in range(0,9):
ind1,=np.where(np.diff(wave_number[j,:])==0)
ind,=np.where(wave_number[j,ind1]==(i+1));
print(j)
ind1,=np.where(np.diff(wave_number[0,:])==0)
ind,=np.where(wave_number[0,ind1]==(i+1));
ind = ind1[ind]
if(len(ind)>2):
mean_rot[j,i]=np.mean(rotation_rate[j,ind[1:-2]]);
std_rot[j,i]=np.std(rotation_rate[j,ind[1:-2]]);
mean_rot[j,i]=np.mean(rotation_rate[0,ind[1:-2]]);
std_rot[j,i]=np.std(rotation_rate[0,ind[1:-2]]);
else:
mean_rot[j,i]=np.nan
std_rot[j,i]=np.nan

h=plt.scatter(np.mgrid[1:10],mean_rot[j,:],s=40,c=u_jets[j]*np.ones(9));

plt.clim((u_jets[0],u_jets[-1]))
plt.xlabel('wave number (per full rotation)');
plt.ylabel('mean rotation rate (rotations per year)');
h=plt.colorbar()
Expand All @@ -115,9 +122,10 @@ def do_analysis01(fileNames,u_jets):
fileName=['/tmp/' + username + '/output.nc']
n_files=len(fileName);

u_jets=[50.]
u_jets=[50., 70., 100.]
thisOne=0

(phase,rotation_rate,wave_number)=fourier_wave_number(fileName,u_jets)
(phase,rotation_rate,wave_number)=fourier_wave_number(fileName)

cmap_lev=64;
map=plt.get_cmap(lut=cmap_lev,name='ocean')
Expand All @@ -140,10 +148,12 @@ def do_analysis01(fileNames,u_jets):
mean_rot[j,i]=np.nan
std_rot[j,i]=np.nan

h=plt.scatter(np.mgrid[1:10],mean_rot[j,:],s=40,c=u_jets[j]*np.ones(9));
h=plt.scatter(np.mgrid[1:10],mean_rot[j,:],s=40,c=u_jets[thisOne]*np.ones(9));
plt.xlabel('wave number (per full rotation)');
plt.ylabel('mean rotation rate (rotations per year)');
h=plt.colorbar()
h.set_label('Jet speed (m/s)')
plt.clim((np.min(u_jets),np.max(u_jets)))
if(thisOne==0):
h=plt.colorbar()
h.set_label('Jet speed (m/s)')

plt.savefig('/tmp/' + username + '/fourier_wave_number.png' ,format='png', dpi=300)
31 changes: 18 additions & 13 deletions python/normal_modes_compare.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ def normal_modes_compare(u_jets,flag):
"""

cmap_lev=64;
map=plt.get_cmap(lut=cmap_lev,name='jet')
map=plt.get_cmap(lut=cmap_lev,name='viridis')
u1=np.linspace(u_jets[0],u_jets[-1],cmap_lev);
int1=sci.interp1d(u1,np.mgrid[1:cmap_lev+1],kind='nearest');

Expand All @@ -42,7 +42,7 @@ def normal_modes_compare(u_jets,flag):
lat_low=65; # lowest
re=5.4155760e7; # radius of saturn in this region (due to squashed spheriod)
#re=5.8232e7;
h_jet=1.2; # standard deviation of jet
h_jet=1.; # standard deviation of jet
lat_jet=77; # latitude of the jet
n_ks=36; # calculate the growth factor for this many k-values

Expand Down Expand Up @@ -125,17 +125,20 @@ def normal_modes_compare(u_jets,flag):
if sigmas_max[n-1,2]<1.e-9:
sigmas_max_i[n-1,2]=np.nan


if flag==1:
h=plt.plot(np.mgrid[1:n_ks+1], \
h=[]
# sigma / k is wave speed in m/s
# k is 2*pi/lambda and lambda is x_len/n
# so this is rotations per year
h.append(plt.plot(np.mgrid[1:n_ks+1], \
(-sigmas_max_i[:,0]/(2.*np.pi*np.mgrid[1:n_ks+1]/x_len)) \
*86400.*365.25/x_len);
h=plt.plot(np.mgrid[1:n_ks+1], \
*86400.*365.25/x_len));
h.append(plt.plot(np.mgrid[1:n_ks+1], \
(-sigmas_max_i[:,1]/(2.*np.pi*np.mgrid[1:n_ks+1]/x_len)) \
*86400.*365.25/x_len,'--');
h=plt.plot(np.mgrid[1:n_ks+1], \
*86400.*365.25/x_len,'--'));
h.append(plt.plot(np.mgrid[1:n_ks+1], \
(-sigmas_max_i[:,2]/(2.*np.pi*np.mgrid[1:n_ks+1]/x_len)) \
*86400.*365.25/x_len,':');
*86400.*365.25/x_len,':'));
plt.ylabel('speed of wave (rotations per year)')
plt.xlabel('number of peaks')
elif flag==2:
Expand All @@ -155,12 +158,13 @@ def normal_modes_compare(u_jets,flag):
h=[h, h2];

if ( len(u_jets) > 1):
row=int1(u_jets[i])
row=int(int1(u_jets[i]))
else:
row=1;

for k in range(len(h)):
h[k].set_color(map(row));
h[k][0].set_color(map(row));


if flag==2:
plt.legend(['Largest','2nd largest','3rd largest','sum'])
Expand All @@ -169,8 +173,9 @@ def normal_modes_compare(u_jets,flag):
if __name__=='__main__':
plt.ion()
plt.figure()

normal_modes_compare([50.],1)
u_jets=[50.,70,100.]
# u_jets=[50.]
normal_modes_compare(u_jets,1)

plt.savefig('/tmp/' + username +'/fourier_wave_number.png' ,format='png', dpi=300)

0 comments on commit 4af76ea

Please sign in to comment.