From 4af76eab68aa880184f33493f3bdf9cbc9a8bee6 Mon Sep 17 00:00:00 2001 From: Paul Connolly Date: Sat, 21 Oct 2023 00:55:35 +0100 Subject: [PATCH] improved scripts --- .gitignore | 2 ++ python/analyse_output.py | 18 +++++++++-------- python/fourier_wave_number.py | 36 ++++++++++++++++++++++------------ python/normal_modes_compare.py | 31 +++++++++++++++++------------ 4 files changed, 53 insertions(+), 34 deletions(-) diff --git a/.gitignore b/.gitignore index f9dbda4..3c31d9e 100644 --- a/.gitignore +++ b/.gitignore @@ -12,3 +12,5 @@ doxygen matlab-dev fortran_testing *.png +__pycache__ + diff --git a/python/analyse_output.py b/python/analyse_output.py index 5ccb065..5fcdabc 100644 --- a/python/analyse_output.py +++ b/python/analyse_output.py @@ -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() @@ -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') diff --git a/python/fourier_wave_number.py b/python/fourier_wave_number.py index 8d7fe1e..e66b7d8 100644 --- a/python/fourier_wave_number.py +++ b/python/fourier_wave_number.py @@ -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 @@ -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). @@ -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]; @@ -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 @@ -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() @@ -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') @@ -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) diff --git a/python/normal_modes_compare.py b/python/normal_modes_compare.py index 23df9b0..2bb6be9 100644 --- a/python/normal_modes_compare.py +++ b/python/normal_modes_compare.py @@ -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'); @@ -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 @@ -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: @@ -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']) @@ -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) \ No newline at end of file