%pylab inline

Populating the interactive namespace from numpy and matplotlib

Problem 1

f=linspace(0,1.1e4,300) #in Hz units
om=2*pi*f

kB=1.3806488e-23; eps0=8.854187e-12; q=1.621e-19 # physical constants
me=9.109e-31; mi=me*1836.152*16 # particle masses, electron, O+
No=1.0e11; Te=1000; Ti=1000 # ionospheric condition
Ce=sqrt(kB*Te/me); Ci=sqrt(kB*Ti/mi) # thermal speeds
he=sqrt(eps0*kB*Te/No/q/q); hi=sqrt(eps0*kB*Ti/No/q/q) # Debye lengths
lam=0.3; k=2*pi/lam # radar Bragg wavelength --- EISCAT ESR (Svalbard)

from scipy.special import dawsn as DawsonI # Dawson integral imported

GordeyevI=(sqrt(pi)*exp(-om*om/2/k/k/Ci/Ci)-2j*DawsonI(om/sqrt(2)/k/Ci))/sqrt(2)/k/Ci
GordeyevE=(sqrt(pi)*exp(-om*om/2/k/k/Ce/Ce)-2j*DawsonI(om/sqrt(2)/k/Ce))/sqrt(2)/k/Ce
sigI=1j*om*eps0*(1-1j*om*GordeyevI)/k/k/hi/hi
sigE=1j*om*eps0*(1-1j*om*GordeyevE)/k/k/he/he

dispersion=1j*om*eps0+sigE+sigI
electronline=abs(1j*om*eps0+sigI)**2*2*real(GordeyevE)/abs(dispersion)**2
ionline=abs(sigE)**2*2*real(GordeyevI)/abs(dispersion)**2
spectrum = ionline + electronline




Te = 2000
Ce=sqrt(kB*Te/me); Ci=sqrt(kB*Ti/mi) # thermal speeds
he=sqrt(eps0*kB*Te/No/q/q); hi=sqrt(eps0*kB*Ti/No/q/q) # Debye lengths


GordeyevI=(sqrt(pi)*exp(-om*om/2/k/k/Ci/Ci)-2j*DawsonI(om/sqrt(2)/k/Ci))/sqrt(2)/k/Ci
GordeyevE=(sqrt(pi)*exp(-om*om/2/k/k/Ce/Ce)-2j*DawsonI(om/sqrt(2)/k/Ce))/sqrt(2)/k/Ce
sigI=1j*om*eps0*(1-1j*om*GordeyevI)/k/k/hi/hi
sigE=1j*om*eps0*(1-1j*om*GordeyevE)/k/k/he/he

dispersion=1j*om*eps0+sigE+sigI
electronline2=abs(1j*om*eps0+sigI)**2*2*real(GordeyevE)/abs(dispersion)**2
ionline2=abs(sigE)**2*2*real(GordeyevI)/abs(dispersion)**2
spectrum2 = ionline2 + electronline2

/usr/local/lib/python2.7/dist-packages/ipython-3.1.0-py2.7.egg/IPython/kernel/__main__.py:19: RuntimeWarning: invalid value encountered in divide
/usr/local/lib/python2.7/dist-packages/ipython-3.1.0-py2.7.egg/IPython/kernel/__main__.py:20: RuntimeWarning: invalid value encountered in divide
/usr/local/lib/python2.7/dist-packages/ipython-3.1.0-py2.7.egg/IPython/kernel/__main__.py:37: RuntimeWarning: invalid value encountered in divide
/usr/local/lib/python2.7/dist-packages/ipython-3.1.0-py2.7.egg/IPython/kernel/__main__.py:38: RuntimeWarning: invalid value encountered in divide

plot(f,spectrum,'k-',label='Te=1000')
plot(f,spectrum2,'k--',label='Te=2000')
legend(loc='best')
xlabel('Freq [Hz]')
ylabel('Spectrum')

<matplotlib.text.Text at 0x7f39c1105590>

Problem 2

f=linspace(0,1.1e4,300) #in Hz units
om=2*pi*f

kB=1.3806488e-23; eps0=8.854187e-12; q=1.621e-19 # physical constants
me=9.109e-31;
mH = me*1836.152
mO = mH*16
fH = 0.1
No=1.0e11; NH = fH*No; NO = (1.0-fH)*No
Te=1000; Ti=1000 # ionospheric condition
Ce=sqrt(kB*Te/me); 
CH=sqrt(kB*Ti/mH) # thermal speeds
CO=sqrt(kB*Ti/mO);
he=sqrt(eps0*kB*Te/No/q/q);
hH=sqrt(eps0*kB*Ti/NH/q/q) # Debye lengths
hO=sqrt(eps0*kB*Ti/NO/q/q)
lam=0.3; k=2*pi/lam # radar Bragg wavelength --- EISCAT ESR (Svalbard)

from scipy.special import dawsn as DawsonI # Dawson integral imported

GordeyevH=(sqrt(pi)*exp(-om*om/2/k/k/CH**2)-2j*DawsonI(om/sqrt(2)/k/CH))/sqrt(2)/k/CH
GordeyevO=(sqrt(pi)*exp(-om*om/2/k/k/CO**2)-2j*DawsonI(om/sqrt(2)/k/CO))/sqrt(2)/k/CO
GordeyevE=(sqrt(pi)*exp(-om*om/2/k/k/Ce/Ce)-2j*DawsonI(om/sqrt(2)/k/Ce))/sqrt(2)/k/Ce
sigH=1j*om*eps0*(1-1j*om*GordeyevH)/k/k/hH**2
sigO=1j*om*eps0*(1-1j*om*GordeyevO)/k/k/hO**2
sigE=1j*om*eps0*(1-1j*om*GordeyevE)/k/k/he/he

dispersion=1j*om*eps0+sigE+sigO+sigH
electronline=abs(1j*om*eps0+sigH+sigO)**2*2*real(GordeyevE)/abs(dispersion)**2
Hline = abs(sigE)**2*2*real(GordeyevH)/abs(dispersion)**2
Oline = abs(sigE)**2*2*real(GordeyevO)/abs(dispersion)**2
spectrum3 = fH*Hline + (1-fH)*Oline + electronline

/usr/local/lib/python2.7/dist-packages/ipython-3.1.0-py2.7.egg/IPython/kernel/__main__.py:29: RuntimeWarning: invalid value encountered in divide
/usr/local/lib/python2.7/dist-packages/ipython-3.1.0-py2.7.egg/IPython/kernel/__main__.py:30: RuntimeWarning: invalid value encountered in divide
/usr/local/lib/python2.7/dist-packages/ipython-3.1.0-py2.7.egg/IPython/kernel/__main__.py:31: RuntimeWarning: invalid value encountered in divide

plot(f,spectrum,'k-',label='O only')
plot(f,spectrum3,'k--',label='H & O')
legend(loc='best')
xlabel('Freq [Hz]')
ylabel('Spectrum')

<matplotlib.text.Text at 0x7f39c0a7c850>

Problem 3

f=linspace(1e1,0.9e7,300) #in Hz units
om=2*pi*f

kB=1.3806488e-23; eps0=8.854187e-12; q=1.621e-19 # physical constants
me=9.109e-31;
mH = me*1836.152
mO = mH*16
fH = 0.5
No=1.0e11; NH = fH*No; NO = (1.0-fH)*No
Te=1000; Ti=1000 # ionospheric condition
Ce=sqrt(kB*Te/me); 
CH=sqrt(kB*Ti/mH) # thermal speeds
CO=sqrt(kB*Ti/mO);
he=sqrt(eps0*kB*Te/No/q/q);
hH=sqrt(eps0*kB*Ti/NH/q/q) # Debye lengths
hO=sqrt(eps0*kB*Ti/NO/q/q)
lam=0.3; k=2*pi/lam # radar Bragg wavelength --- EISCAT ESR (Svalbard)

from scipy.special import dawsn as DawsonI # Dawson integral imported

GordeyevH=(sqrt(pi)*exp(-om*om/2/k/k/CH**2)-2j*DawsonI(om/sqrt(2)/k/CH))/sqrt(2)/k/CH
GordeyevO=(sqrt(pi)*exp(-om*om/2/k/k/CO**2)-2j*DawsonI(om/sqrt(2)/k/CO))/sqrt(2)/k/CO
GordeyevE=(sqrt(pi)*exp(-om*om/2/k/k/Ce/Ce)-2j*DawsonI(om/sqrt(2)/k/Ce))/sqrt(2)/k/Ce
sigH=1j*om*eps0*(1-1j*om*GordeyevH)/k/k/hH**2
sigO=1j*om*eps0*(1-1j*om*GordeyevO)/k/k/hO**2
sigE=1j*om*eps0*(1-1j*om*GordeyevE)/k/k/he/he

dispersion=1j*om*eps0+sigE+sigO+sigH

loglog(f,abs(dispersion),'k-')
mx = dispersion.max()
#ylim((1e-3*mx, 1.1*mx));

The two ion modes are impossible to distinguish. The plasma wave is clearly around 3 MHz. There is a "kink" in the dispersion relation around 40 kHz. I could not find the third mode, but it is also probably very close to 40 kHz.

Problem 5

fn = '/home/bhardin2/ece458/JRO_10min_Valley2015014_11.npz.npz'
npzfile = load(fn)
print npzfile.keys()
dat = npzfile['data_ch0']
hts = npzfile['heights']
t = npzfile['epochtime']

['epochtime', 'DST', 'time_zone', 'heights', 'data_ch0', 'data_ch1']

# Plot of power vs height
Nh = len(hts)
P = zeros(Nh)
for hi in range(Nh):
    Estub = dat[:,hi]
    P[hi] = sum(abs(Estub)**2)
semilogx(P,hts,'k-')
xlabel('Power')
ylabel('Altitude [km]')

<matplotlib.text.Text at 0x7f39b8b7fd90>

# Create spectrum at different heights
n = 80 # length of FFT chunks
M = int(len(t)/n) # number of FFT chunks

S = zeros((Nh,n))

for hi in range(Nh):

    F2 = zeros(n)
    for j in range(M):
        Estub = dat[n*j:n*(j+1),hi]
        F = fft.fft(Estub)/n
        F2 += abs(F)**2
    F2 = F2/M
    F2s = fft.fftshift(F2)
    x = arange(n)
    xs = fftshift(x)
    xs[xs>=n/2] = xs[xs>=n/2] - n

    T = t[1]-t[0]
    df = 1/(T*n)
    f = xs*df
    S[hi,:] = log10(abs(F2s)**2)
      

F,H = meshgrid(f,hts)
pcolormesh(F,H,S,cmap='hot')
axis('tight')
xlabel('Frequency [Hz]')
ylabel('Altitude [km]')
cb = colorbar()
cb.set_label('log10(Power)')
clim((0,10))

None of the spectra seem to be exhibiting the double hump expected from incoherent scatter. This is probably because of the presence of coherent scatter from 150 km echoes and the electrojet. Below we show examples from a few altitudes.

hvec = linspace(100,200,5)

for h in hvec:

    figure()
    n = 80 # length of FFT chunks
    M = int(len(t)/n) # number of FFT chunks
    #M = 50
    #print 'Averaging over %.1f seconds' % (t[M*n]-t[0])
    #h = 146.
    hi = argmin(abs(hts - h)) # Extract index closest to h
    #print 'Height: %.1f km' % hts[hi]

    F2 = zeros(n)
    for j in range(M):
        Estub = dat[n*j:n*(j+1),hi]
        F = fft.fft(Estub)/n
        F2 += abs(F)**2
    F2 = F2/M
    F2s = fft.fftshift(F2)
    x = arange(n)
    xs = fftshift(x)
    xs[xs>=n/2] = xs[xs>=n/2] - n

    T = t[1]-t[0]
    df = 1/(T*n)
    f = xs*df

    semilogy(f,F2s,'k')
    #xlabel('FFT bin')
    xlabel('Frequency [Hz]')
    ylabel('Magnitude-squared')
    axis('tight');
    title('%.1f km altitude' % hts[hi])

npzfile.close()

import subprocess
import shutil
nb_name = 'ECE458_HW9.ipynb' # There's probably a way to get this programmatically
cmd = ['ipython', 'nbconvert', '--to', 'slides', nb_name, \
       '--reveal-prefix', '"http://cdn.jsdelivr.net/reveal.js/2.6.2"',]
#      '--template','output_toggle']
# The --template command points to a custom template which opens and closes the input code cells with a click
subprocess.call(cmd)
html_name = '%s.slides.html' % (nb_name.split('.')[0])
shutil.copy(html_name, '/home/bhardin2/public_html/slideshows/')