<jupyter_start><jupyter_text>Rotary spectra - test caseA test case to check the validity of rotary spectra code, using falsified circular velocity data.Based on <NAME>. (1997). Data Analysis Methods in Physical Oceanography. pp.427-432,494-500; <NAME>. (1972) Deep Sea Res. 833-846, and https://pyoceans.github.io/python-oceans/ocfis.html.  Imports<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
from numpy.random import seed
from numpy.random import rand
import scipy as sci
import scipy.signal as sig
import scipy.integrate as integ
import cmath
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>Create velocity dataFalsify circular, horizontal velocity data, in CCW direction.<jupyter_code>r = 40000                                       # radius of circle, m
d = int(1e3)                                    # number of data points
T = int(1e6)                                    # total time for motion, s
t = np.linspace(0,T,d)                          # time range, s
dt = t[d-1]-t[d-2]                              # time step size, s
fs = 1/dt                                       # sampling rate, Hz
w = 0.0003                                      # angular velocity, rad/s
u = -r*w*np.sin(w*t)                            # u velocity of particle
v = r*w*np.cos(w*t)                             # y velocity of particle

M2_f = 0.000022364278                           # tidal frequency for M2, Hz
M2_w = 2*np.pi*M2_f                             # angular frequency for M2, rad/s
M2_u = 0.6*( -r*M2_w*np.sin(M2_w*t) )           # M2 velocity, 
M2_v = 0.4*( r*M2_w*np.cos(M2_w*t) )

seed(1)                                         # generate noise for data
values = rand(d)
noise_w = values*0.0001
noise = np.cos(noise_w*t)

u_tot = (u + M2_u)                              # circular velocity + offset tidal velocities w/ noise
v_tot = (v + M2_v)
u_tot = u_tot+noise
v_tot = v_tot+noise

u_tot = u_tot - np.mean(u_tot)                  # remove mean for spectra
v_tot = v_tot - np.mean(v_tot)<jupyter_output><empty_output><jupyter_text>Rotary spectra (NumPy, not averaged)Defines a function to return the CW, CCW, cross, and quadrature spectra for any u and v velocity vectors.<jupyter_code>def spec_rot(u, v):
    fu, fv = list(map(np.fft.fft, (u, v)))             # individual components Fourier series
    pu = fu * np.conj(fu)                              # auto-spectra of the scalar components
    pv = fv * np.conj(fv)
    puv = fu.real * fv.real + fu.imag * fv.imag        # cross spectrum, not used in the current function version
    quv = -fu.real * fv.imag + fv.real * fu.imag       # quadrature spectrum
    cw = (pu + pv - (2*quv)) / 2                       # rotatory components
    ccw = (pu + pv + (2*quv)) / 2
    F = np.fft.fftfreq(n=d,d=dt)                       # frequency range (two-sided)
    return puv, quv, cw, ccw, F
puv, quv, cw, ccw, f = spec_rot(u_tot,v_tot)           # get rotary components (cw and ccw) and frequency range (f)
half_idx = int((d/2))                                  # discard half spectrum for real data 
cw_real = cw[1:half_idx]*2                             # mult. amplitude by 2 to account for discarded data
ccw_real = ccw[1:half_idx]*2     
f_real = f[1:half_idx]                                 # real frequency range (up to Nyquist)
fig,ax = plt.subplots(1,1,figsize=(6,6))               # plot CCW and CW components
ax.loglog(f_real,ccw_real,label='ccw')
ax.loglog(f_real,cw_real,label='cw',ls='--')
ax.axvline(x=2.23e-5,ls='--',lw=0.7,c='g',label=f'$M_2$: 2.23e-5 Hz')
ax.axvline(x=4.77e-5,ls='--',lw=0.7,c='purple',label='CCW velocity: 4.77e-5 Hz')
ax.set_title('Rotary spectra for $u,v$ (no averaging)')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel(f'Energy density [$(m/s)^2/Hz$]')
ax.legend(loc='lower left',fontsize=9)
plt.show()<jupyter_output><empty_output><jupyter_text>Spectra appear to be correct for falsified purely CCW circular velocity data (2nd peak), M2 tide mostly in $u$ (1st peak), and simulated noise. Just needs to be averaged, and checked for normalisation.  Averaged rotary spectra (SciPy)<jupyter_code>def spec_rot_avg(u, v):
    nps = 256           # nperseg
    win = 'parzen'      # window
    puf, pu = sig.welch(u,fs=fs,window=win,nperseg=nps,return_onesided=False)           # auto-spectrum for u
    pvf, pv = sig.welch(v,fs=fs,window=win,nperseg=nps,return_onesided=False)           # auto-spectrum for v
    cuvf, cuv = sig.csd(v,u,fs=fs,window=win,nperseg=nps,return_onesided=False)         # cross spectra (u,v --> v,u)
    quv = cuv.imag                                      # quadrature spectrum, imaginary part of cross spectra
    cw = ((pu + pv) - (2*quv)) / 8                      # rotatory components
    ccw = ((pu + pv) + (2*quv)) / 8
    F = puf                                             # frequency range (two-sided)
    return cuv, quv, cw, ccw, F
cuv_avg, quv_avg, cw_avg, ccw_avg, f_avg = spec_rot_avg(u_tot,v_tot)  # get rotary components (cw and ccw) and frequency range (f)
half_idx_avg = int((len(cw_avg)/2))                                   # discard half spectrum for real data 
cw_real_avg = cw_avg[1:half_idx_avg]*2                                # mult. amplitude by 2 to account for discarded data
ccw_real_avg = ccw_avg[1:half_idx_avg]*2     
f_real_avg = f_avg[1:half_idx_avg]                                    # real frequency range (up to Nyquist)
# plot CCW and CW components
fig,ax = plt.subplots(1,1,figsize=(6,6))
ax.loglog(f_real_avg,ccw_real_avg,label='ccw')
ax.loglog(f_real_avg,cw_real_avg,label='cw',ls='--')
ax.axvline(x=2.23e-5,ls='--',lw=0.7,c='g',label=f'$M_2$: 2.23e-5 Hz')
ax.axvline(x=4.77e-5,ls='--',lw=0.7,c='purple',label='CCW velocity: 4.77e-5 Hz')
ax.set_title('Rotary spectra for $u,v$ (with averaging)')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel(f'Energy density [$(m/s)^2/Hz$]')
ax.legend(loc='upper right',fontsize=9)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.125923438549362, 'fraction_numerical': 0.027535258562793822, 'mean_word_length': 2.3372549019607844, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 20, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '537370', 'license': 'MIT', 'n_tokens_mistral': 2243, 'n_tokens_neox': 2112, 'n_words': 522, 'path': 'archive/.ipynb_checkpoints/rotary_test-checkpoint.ipynb', 'repo_name': 'kurtisanstey/internal_waves_barkley_canyon', 'size': '5964'}