#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Sun Feb 21 21:21:21 2021 @author: Xiaoting Yang """ ## Rossby_wave_ray_tracing_ala_Hoskins_Karoly_1981 import numpy as np import matplotlib.pyplot as plt from cmath import sqrt # define global variables Omega=2*np.pi/(24*3600.0) # rotational rate of the earth a=6371000.0 # radius of the earth tmax=86400*30 # integration length dt=200.0 # time step nt=int(np.ceil(tmax/dt)) # number of time steps n=2 # integer zonal wavenumber w=Omega/30 # angular velocity phijet=45.0 # jet latitude ujet=0 # mid-latitude jet velocity lon_init=0.0 # initial longitude lat_init=5.0 # initial latitude # define functions to be used def sinh(x): res=0.5*(np.exp(x)-np.exp(-x)) return res def cosh(x): res=0.5*(np.exp(x)+np.exp(-x)) return res def tanh(x): res=np.divide(sinh(x),cosh(x)) return res def sech(x): res=np.divide(1,cosh(x)) return res def phi_to_y(phi): # latitude in degrees to northward distance in meter res=np.divide(1+np.sin(phi*np.pi/180),np.cos(phi*np.pi/180)) res=a*np.log(res) return res def y_to_phi(y): # distance in meter to latitude in degrees res=np.arcsin(tanh(y/a))*180/np.pi return res def ubar_m(y): # the mean flow in which the ray propagates phi=y_to_phi(y) res=a*w*sech(y/a)+ ujet*np.divide(np.exp(-(phi-phijet)**2/(15)**2)\ ,np.cos(phi*np.pi/180)); return res def beta_m(y): # Mercator effective beta of Hoskins and Karoly yjet=phi_to_y(phijet) dely=800*1e3; angmomflw = (Omega/a)*np.multiply(sech(y/a),(2*sech(y/a) \ - (0.1*tanh(y/a)**2) + (0.1*sech(y/a)**2))); f = 1/a**2 + 2/dely**2 - 4*(y-yjet)**2/dely**4 \ - 4*(y-yjet)*tanh(y/a)/(dely**2*a); jet = ujet * f * np.multiply(cosh(y/a), np.exp(-(y-yjet)**2/dely**2)); res = angmomflw + jet; return res def dbeta_dy(y): yjet = phi_to_y(phijet); dely = 800*1e3; angmomflw = -(Omega/(a**2)) * np.multiply(np.multiply(sech(y/a) ,tanh(y/a)), ((sech(y/a)**2/2 + 4*sech(y/a)-(tanh(y/a)**2)/10))); f = 1/a**2 + 2/dely**2 - 4*(y-yjet)**2/dely**4 - 4*(y-yjet)*tanh(y/a)/(dely**2*a); dfdy = -8*(y-yjet)/dely**4 - 4*tanh(y/a)/(dely**2*a) - 4*(y-yjet)*sech(y/a)**2/(a**2*dely**2); jet = ujet * np.multiply(np.multiply((dfdy - f*tanh(y/a)/a - 2*(y-yjet)*f/dely**2), cosh(y/a)), np.exp(-(y-yjet)**2/dely**2)); res = angmomflw+jet return res # Initilize arrays that will save location and wavenumbers of ray along path: x=np.zeros((nt,1),dtype=np.complex128)*np.nan; y=np.zeros((nt,1),dtype=np.complex128)*np.nan; phi_hist=np.zeros((nt,1),dtype=np.complex128)*np.nan; l_hist=np.zeros((nt,1),dtype=np.complex128)*np.nan dldt_hist=np.zeros((nt,1),dtype=np.complex128)*np.nan # specify initial location of ray: x[0]=lon_init*a*np.pi/180; y[0]=phi_to_y(lat_init); k=n/a; # wave number (should this be multiplied by 2pi?) l_hist[0]=sqrt(beta_m(y[0])/ubar_m(y[0])-k**2) # integrate along ray path: for i in range(nt-1): Ks=sqrt(beta_m(y[i])/ubar_m(y[i])) dldt=(k*dbeta_dy(y[i])/Ks**2)+(k*w*sech(y[i]/a)*tanh(y[i]/a)) dldt_hist[i] = 1*dldt l_hist[i+1] = 1*l_hist[i] + dldt*dt l_int = (l_hist[i+1]+l_hist[i])/2 cg_x = ubar_m(y[i])+beta_m(y[i])*(k**2-l_int**2)/(k**2+l_int**2)**2 cg_y = 2*beta_m(y[i])*k*l_int/(k**2+l_int**2)**2 x[i+1] = 1*x[i] + cg_x*dt y[i+1] = 1*y[i] + cg_y*dt phi_hist=y_to_phi(y) # make figure with four panels: Time=np.linspace(0,tmax,nt) plt.figure(figsize=(12,8)) plt.clf() # lon/lat polar view of trajectories # ------------------------------------------------------- ax1=plt.subplot(221); xtemp=np.array(y_to_phi(y)-90.0)*np.sin(x/a) ytemp=np.array(y_to_phi(y)-90.0)*np.cos(x/a) ax1.plot(xtemp,ytemp,linewidth=1.2) pos=ax1.get_position(); pax=plt.axes(pos,projection='polar',facecolor='none') pax.set_yticklabels([]) ax1.set_aspect(1); ax1.axis('Off'); # l and dl/dt as function of time # ------------------------------------------------------- ax2a=plt.subplot(222); color='tab:blue'; ax2a.set_xlabel('Time (days)'); ax2a.set_ylabel('l',color=color); ax2a.plot(Time/86400.0,l_hist,color=color); ax2a.ticklabel_format(style='sci') ax2b=ax2a.twinx(); color='tab:red'; ax2b.set_ylabel('dldt',color=color); ax2b.plot(Time/86400.0,dldt_hist,color=color); ax2b.ticklabel_format(style='sci') # Ubar as function of latitude # ------------------------------------------------------- ax3=plt.subplot(223); yplot=np.linspace(0,1.5*np.pi*a/2,100); phi=y_to_phi(yplot); ubar_plot=np.multiply(ubar_m(yplot),np.cos(phi*np.pi/180)); ax3.plot(phi,ubar_plot) ax3.set_xlabel('Latitude'); ax3.set_ylabel('Ubar'); ax3.set_title(r'$\overline{u}(\phi)$, with $u_{jet}=$'+str(ujet)); K_s_plot=np.sqrt(np.divide(beta_m(y),ubar_m(y))) # latitude of ray and K_s as function of time along path # ------------------------------------------------------- ax4a=plt.subplot(224); color='tab:blue'; ax4a.set_xlabel('Time (days)'); ax4a.set_ylabel('Ray latitude',color=color); ax4a.plot(Time/86400.0,y_to_phi(y),color=color); ax4b=ax4a.twinx(); color='tab:red'; ax4b.set_ylabel(r'$R_E K_s$',color=color); ax4b.plot(Time/86400.0,a*K_s_plot,color=color); # plt.tight_layout() # dont use, causes problems plt.show() #plt.savefig('./Figures/Rossby_ray_tracing_k'+str(n)+'a.png') # second figure: l as a function of latitude: # ------------------------------------------- l_formula=0*l_hist; for i in range(nt): l_formula[i]=sqrt(K_s_plot[i]**2-k**2) plt.figure(figsize=(6,6)); plt.clf(); plt.scatter(phi_hist[::100],l_hist[::100],c=y[::100],cmap='jet'\ ,label="ray tracing; color=latitude"); temp1=min(l_hist.real); temp2=max(l_hist.real); plt.ylim([1.1*temp1-0.1*temp2,1.1*temp2-0.1*temp1]) plt.scatter(phi_hist[::100],l_formula[::100],edgecolors='k',facecolors='None'\ ,label="steady state from $K_s$") plt.xlabel('Latitude'); plt.ylabel('l') plt.legend() plt.title("l vs latitude, ray tracing vs steady calculation") plt.show() #plt.savefig('./Figures/Rossby_ray_tracing_k'+str(n)+'b.png')