
%config InlineBackend.figure_format='retina'
import sys
from datetime import datetime
from rich.console import Console
console = Console(highlight=False)#,force_jupyter=False,width=999)


from cdasws.cdasws import CdasWs, TimeInterval
from cdasws.datarepresentation import DataRepresentation as dr
from cdasws.doi import get_doi_landing_page_url, get_doi_badge_url, get_doi_hdp_url

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

from IPython.core.display import HTML

today = datetime.now().strftime('%Y-%B-%d')

cdas = CdasWs()

def display_link(
    url: str,
    text: str) -> None:
    display(HTML(f'<a href="{url}" target="_blank">{text}</a>'))

def display_image_link(
    url: str,
    image_url: str) -> None:
    display(HTML(f'<a href="{url}" target="_blank"><img src={image_url}/></a>'))
    

# Injects CSS to target the notebook's standard output stream container
display(HTML('''
<style>
div.output_area pre {
    font-size: 16px !important;
}
</style>
'''))


cdas = CdasWs()


# Fetch all observatory groups
obs_groups = cdas.get_observatory_groups()

# Print the names
print("Observatory -> IDs\n")
for indx, group in enumerate(obs_groups):
    #print("'%s'" %(group['Name']),',', end="")
    print("\033[1;32m%s\033[0m -> %s" %(group['Name'], group['ObservatoryId']))
    #print(f'{group["Name"]}')#"'%s'" %(group['Name']),',', end="")

Observatory -> IDs

ACE -> ['AC', 'ACE', 'OMNI (1AU IP Data)']
AMPTE -> ['AMPTE-CCE']
ARTEMIS -> ['THB', 'THC']
AWE -> ['AWE']
Alouette -> ['Alouette-2']
Apollo -> ['ALSEP']
Arase (ERG) -> ['Arase(ERG)', 'ERG']
Balloons -> ['PMC Turbo', 'bar_1A', 'bar_1B', 'bar_1C', 'bar_1D', 'bar_1G', 'bar_1H', 'bar_1I', 'bar_1J', 'bar_1K', 'bar_1M', 'bar_1N', 'bar_1O', 'bar_1Q', 'bar_1R', 'bar_1S', 'bar_1T', 'bar_1U', 'bar_1V', 'bar_2A', 'bar_2B', 'bar_2C', 'bar_2D', 'bar_2E', 'bar_2F', 'bar_2I', 'bar_2K', 'bar_2L', 'bar_2M', 'bar_2N', 'bar_2O', 'bar_2P', 'bar_2Q', 'bar_2T', 'bar_2W', 'bar_2X', 'bar_2Y', 'bar_3A', 'bar_3B', 'bar_3C', 'bar_3D', 'bar_3E', 'bar_3F', 'bar_3G', 'bar_4A', 'bar_4B', 'bar_4C', 'bar_4D', 'bar_4E', 'bar_4F', 'bar_4G', 'bar_4H', 'bar_5A', 'bar_6A', 'bar_7A']
BioSentinel -> ['BioSentinel']
CNOFS -> ['CNOFS']
CRRES -> ['CRRES Satellite']
Cassini -> ['CASSINI', 'Cassini']
Cluster -> ['C1', 'C2', 'C3', 'C4', 'CL', 'CT']
DMSP -> ['DMSP', 'dmsp', 'dmsp-f12', 'dmsp-f13', 'dmsp-f14', 'dmsp-f15', 'dmsp-f16', 'dmsp-f17', 'dmsp-f18', 'dmsp-f6', 'dmsp-f7', 'dmsp-f8', 'dmsp-f9']
DSCOVR -> ['DSCOVR']
Demonstration and Science Experiments (DSX) -> ['DSX']
Dynamics Explorer -> ['DE', 'DE-1', 'DE-2', 'DE1']
Equator-S -> ['EQ']
FAST -> ['FAST']
Formosat -> ['FS-5', 'ROCSAT-1(FORMOSAT-1)']
GOES -> ['G0', 'G10', 'G12', 'GOES-08', 'GOES-09', 'GOES-10', 'GOES-11', 'GOES-12', 'GOES-13', 'GOES-14', 'GOES-15', 'GOES-16', 'GOES-17', 'GOES-18', 'GOES-19', 'GOES10', 'GOES11', 'GOES12', 'GOES13', 'GOES14', 'GOES15', 'GOES16', 'GOES17', 'GOES18', 'GOES8', 'GOES9', 'GOES_10', 'GOES_11', 'GOES_12', 'GOES_6', 'GOES_7', 'GOES_8', 'GOES_9', 'g11']
GOLD -> ['GOLD']
GPS -> ['GPS']
Genesis -> ['Genesis']
Geosynchronous Investigations -> ['STPSat-6']
Geotail -> ['GEOTAIL']
Ground-Based Investigations -> ['CANOPUS', 'DARN', 'FMI', 'SESAME', 'THG_L2']
Hawkeye -> ['HK']
Helio Ephemeris -> ['BEPICOLOMBO', 'BORRELLY', 'COMETGS', 'COMETHMP', 'DAWN', 'EARTH', 'GALILEO', 'GIACOBINI', 'GIOTTO', 'HALEBOPP', 'HALLEY', 'HYAKUTAKE', 'JUICE', 'JUNO', 'JUPITER', 'MARS', 'MERCURY', 'MOON', 'MSL', 'NEPTUNE', 'PHOBOS2', 'PLUTO', 'PSYCHE', 'ROSETTA', 'SAKIGAKE', 'SATURN', 'SOLAR-1', 'SUISEI', 'URANUS', 'VENUS', 'WILD2']
Helios -> ['HELIOS', 'HELIOS-1', 'HELIOS-2', 'HELIOS1', 'HELIOS2', 'Helios-1', 'Helios-2']
IBEX -> ['IBEX']
ICON -> ['ICON', 'ICON ']
IMAGE -> ['IM', 'IMAGE']
IMAP -> ['IMAP']
IMP (All) -> ['I8', 'IMP-8', 'IMP7', 'IMP8', 'OMNI (1AU IP Data)']
ISEE -> ['ISEE 1', 'ISEE-1', 'ISEE-3', 'ISEE1', 'ISEE2']
ISIS -> ['ISIS-1', 'ISIS-2']
ISS -> ['AWE', 'ISS']
Interball -> ['INTERBALL-AURORAL', 'INTERBALL-GROUND', 'INTERBALL-TAIL']
LANL -> ['L1', 'LANL-97A', 'LANL1989_046', 'LANL1990_095', 'LANL1991_080', 'LANL1994_084', 'Lanl-01A', 'Lanl-02A', 'Lanl-97A', 'Lanl1989_046', 'Lanl1990_095', 'Lanl1991_080', 'Lanl1994_084']
MAVEN -> ['MAVEN']
MESSENGER -> ['MESSENGER']
MMS -> ['MMS1', 'MMS2', 'MMS3', 'MMS4', 'mms1', 'mms2', 'mms3', 'mms4']
Mariner -> ['MARINER2', 'Mariner2']
Mars Global Surveyor (MGS) -> ['MGS']
Mars Science Laboratory (MSL) -> ['MSL']
Moon -> ['BGM1']
NOAA -> ['MetOp1', 'MetOp2', 'MetOp3', 'NOAA05', 'NOAA07', 'NOAA08', 'NOAA10', 'NOAA12', 'NOAA14', 'NOAA15', 'NOAA16', 'NOAA18', 'NOAA19']
New Horizons -> ['NEW_HORIZONS', 'NH']
OMNI (Combined 1AU IP Data; Magnetic and Solar Indices) -> ['OMNI (1AU IP Data)']
POES/MetOp -> ['MetOp1', 'MetOp2', 'MetOp3', 'NOAA15', 'NOAA16', 'NOAA18', 'NOAA19']
Parker Solar Probe (PSP) -> ['PSP', 'PSP_FLD']
Pioneer -> ['PIONEER10', 'PIONEER11', 'Pioneer10', 'Pioneer11', 'Pioneer6', 'Pioneer7', 'PioneerVenus']
Polar -> ['PO', 'POLAR', 'Polar', 'po']
REACH -> ['REACH', 'reach']
SAMPEX -> ['SAMPEX']
SOHO -> ['SOHO']
SOLAR -> ['SOLAR-1']
ST5 -> ['ST5']
STEREO -> ['Ahead', 'Behind', 'STA', 'STB', 'STEREO', 'STEREOA', 'STEREOB', 'sta', 'stb']
Smallsats/Cubesats -> ['CIRBE', 'CSSWE', 'Munin', 'SNOE', 'ac6a', 'ac6b', 'ela', 'elb']
Solar Orbiter -> ['SOLO', 'solo']
Sounding Rockets -> ['Dynamo-2', 'Endurance', 'RENU2', 'TOOWINDY']
THEMIS -> ['THA', 'THB', 'THC', 'THD', 'THE', 'THEMIS', 'THG', 'THG_L1', 'THG_L2']
TIMED -> ['TIMED']
TRACERS -> ['TRACERS1', 'TRACERS2', 'TS1', 'TS2']
TSS-1R -> ['TSS-1R']
TWINS -> ['TWINS-', 'TWINS1', 'TWINS2']
Ulysses -> ['ULYSSES', 'Ulysses']
Van Allen Probes (RBSP) -> ['RBSP', 'RBSP-A', 'RBSP-A-RBSPICE', 'RBSP-B', 'RBSP-B-RBSPICE', 'RBSPA', 'RBSPB', 'rbsp-a', 'rbsp-b', 'rbspa', 'rbspb']
Voyager -> ['VG1', 'VG2', 'VGR-1', 'VGR-2', 'VOYAGER1', 'VOYAGER2', 'Voyager 1', 'Voyager 2', 'Voyager1', 'Voyager2']
Wind -> ['OMNI (1AU IP Data)', 'WI', 'WIND', 'Wi', 'Wind', 'wi']


# Fetch all instrument types
instrument_types = cdas.get_instrument_types()

# print the names
for group in instrument_types:
    print("'%s'" %(group['Name']))

' Particles (space)'
'Activity Indices'
'Dust and Debris'
'Electric Fields (space)'
'Engineering'
'Ephemeris/Attitude/Ancillary'
'Ground-Based HF-Radars'
'Ground-Based Imagers'
'Ground-Based Magnetometers, Riometers, Sounders'
'Ground-Based VLF/ELF/ULF, Photometers'
'Housekeeping'
'Imagers (space)'
'Imaging and Remote Sensing (ITM/Earth)'
'Imaging and Remote Sensing (Magnetosphere/Earth)'
'Imaging and Remote Sensing (Sun)'
'Magnetic Fields (Balloon)'
'Magnetic Fields (space)'
'Particle Flux (space)'
'Particles (space)'
'Plasma and Solar Wind'
'Radio and Plasma Waves (space)'
'Space Weather Monitor'


### user-specified ###
myobservatory=['Ulysses']
myinstrument =['Magnetic Fields (space)']


all_datasets = []
for obs in myobservatory:
    print('for %s' %(obs))
    print('Dataset with the link -> Label')
    datasets = cdas.get_datasets(observatoryGroup=obs,instrumentType=myinstrument)
    if datasets:
        all_datasets.extend(datasets)
        
    lines=[]
    for index, dataset in enumerate(datasets):
        dataset_id = dataset["Id"]
        dataset_label = dataset["Label"]
        dataset_note = dataset['Notes']
        #dataset_note = dataset['DatasetLink']['Url']
        #skt_url = f"https://cdaweb.gsfc.nasa.gov/pub/software/cdawlib/0SKELTABLES/{dataset_id}_00000000_v01.skt"
        #response = requests.get(skt_url)
        #skeleton_table_text = response.text

        #doi_info = dataset["Doi"]
        #
        #if doi_info:
        #    # Generate the landing page metadata URL
        #    metadata_url = get_doi_landing_page_url(doi_info)
        #else:
        #    # Fallback if no DOI info is provided directly in the dataset object
        #    metadata_url = f"https://nasa.gov{dataset_id}"

        #print('\033[1;32m%s\033[0m -> %s, %s' %(dataset_id, dataset_label,dataset_note))#, display_link(metadata_url, 'metadata')))
        #print("'\033[1m%s\033[0m' -> %s" %(dataset_id, dataset_label))#, display_link(metadata_url, 'metadata')))
        #print("%s" %(dataset_note))#, display_link(metadata_url, 'metadata')))
        #print('')
        #display(Markdown('\033[1;32m%s\033[0m -> %s' %(dataset_id, dataset_label)))#, display_link(metadata_url, 'metadata')))
        #print(f'{index}. {dataset_id}, {dataset_label}')
        lines.append("[bold green][link=%s]%s[/link][/] -> %s" %(dataset_note,dataset_id,dataset_label))
    console.print("\n".join(lines),justify='left', soft_wrap=True)

    print('='*100)
    print(' ')

for Ulysses
Dataset with the link -> Label

UY_1MIN_VHM -> Ulysses VHM 1 minute. - A. Balogh (Imperial College, London, UK)                                    
UY_1SEC_VHM -> Ulysses VHM 1 second. - A. Balogh (Imperial College, London, UK)                                    
UY_M1_VHM -> Ulysses VHM 1 hour average. - A. Balogh (Imperial College, London, UK)

====================================================================================================


### user-specified ###
mydataset=['UY_1SEC_VHM'] #the order should be consistent with the list above.


mydois=[]
for index, dataset in enumerate(all_datasets):
    dataset_id = dataset["Id"]
    dataset_label = dataset["Label"]

    if dataset_id in mydataset:  # select this dataset for use below
        doi=dataset['Doi']
        #print(doi)
        mydois.append(dataset['Doi'])
        print(f'\nSelected dataset: {dataset_id} (doi:{doi})')
        print(f"Available Period: {dataset['TimeInterval']}")
        display_image_link(get_doi_landing_page_url(doi), get_doi_badge_url(doi))
        display_link(get_doi_hdp_url(doi), 'HDP Description')
        display_link(dataset['Notes'], 'CDAWeb Notes')
        display(HTML(f'<b>Data Citation</b>: {cdas.get_citation(doi)}.  Accessed on {today}.'))
        
        #variables = cdas.get_variables(doi)
        #var_names = []
        #for index, variable in enumerate(variables):
        #    name = variable['Name']
        #    var_names.append(name)
        #    print("'%s'" %(variable['Name']),"->", variable['LongDescription'])
        
        print('='*100)
        print(' ')

Selected dataset: UY_1SEC_VHM (doi:10.48322/mzqv-2n56)
Available Period: {'Start': '1990-10-25T16:50:45.000Z', 'End': '2009-06-30T18:49:11.000Z'}

====================================================================================================


for indx,mydoi in enumerate(mydois):
    variables = cdas.get_variables(mydoi)
    print('for %s' %(mydataset[indx]))
    
    search_window = TimeInterval("2000-01-01T00:00:00Z", "2005-01-01T00:00:00Z")
    filtered_inventory = cdas.get_inventory(mydataset[indx], timeInterval=search_window)

    for interval in filtered_inventory:
        print(f"Available from: \033[1;32m{interval.start}\033[0m to \033[1;32m{interval.end}\033[0m")
    print(' ')
    var_names = []
    for index, variable in enumerate(variables):
        name = variable['Name']
        var_names.append(name)
        print("\033[1m'%s'\033[0m" %(variable['Name']),"->", variable['LongDescription'])
    #    print(f'{index}. {name}, {variable["LongDescription"]}')
    print('='*100)
    print(' ')

for UY_1SEC_VHM
Available from: 2000-01-01 00:00:01+00:00 to 2000-05-27 23:25:22+00:00
Available from: 2000-05-28 02:14:02+00:00 to 2002-10-30 23:40:44+00:00
Available from: 2002-10-31 10:09:02+00:00 to 2003-03-05 23:50:21+00:00
Available from: 2003-03-06 01:28:30+00:00 to 2003-07-07 22:57:50+00:00
Available from: 2003-07-08 00:59:02+00:00 to 2004-11-21 15:28:09+00:00
Available from: 2004-11-22 15:34:36+00:00 to 2004-12-31 23:59:59+00:00
 
'B_RTN' -> Magnetic field components in RTN coordinates
'B_MAG' -> Magnetic field magnitude
====================================================================================================


# assign parameteprs
#dataset_id = 'AC_H0_SWE'
myvariables  = {}
mystart_time = {}
myend_time   = {}

### user-specified ###
myvariables[mydataset[0]]=['B_RTN']
mystart_time[mydataset[0]] = '2004-02-13T00:00:00Z'
myend_time[mydataset[0]]   = '2004-02-16T00:00:00Z'
#######################


# request the data representation as Spacepy
data={}
for dataset in mydataset:
    print(dataset)
    status, data[dataset] = cdas.get_data(dataset, myvariables[dataset], mystart_time[dataset], myend_time[dataset], dataRepresentation = dr.XARRAY)
    if status['http']['status_code'] == 200:
        print('Donwload Success!')
    else:
        print(f"Error! HTTP Status Code: {status['http']['status_code']}")

### change the name of the first index to "time" (hope the first variable is Epoch-related).
for dataset in mydataset:
    keys=list(data[dataset].indexes.keys())
    data[dataset] = data[dataset].rename({keys[0]: "time"})

UY_1SEC_VHM
Donwload Success!


print(data)

{'UY_1SEC_VHM': <xarray.Dataset>
Dimensions:  (time: 256862, B_RTN_labels: 3)
Coordinates:
  * time     (time) datetime64[ns] 2004-02-13T00:00:00.660000 ... 2004-02-15T...
Dimensions without coordinates: B_RTN_labels
Data variables:
    B_RTN    (time, B_RTN_labels) float32 -0.102 0.39 -0.075 ... -1.54 1.328
Attributes: (12/27)
    Acknowledgement:             ['Please acknowledge the Principal', 'Invest...
    ADID_ref:                    ['NSSD0241']
    Data_type:                   ['1SEC>1 second resolution']
    Data_version:                ['1']
    Descriptor:                  ['VHM>Vector Helium Magnetometer']
    Discipline:                  ['Space Physics>Interplanetary Studies']
    ...                          ...
    Time_resolution:             ['1 second resolution']
    TITLE:                       ['Ulysses VHM 1 second', 'magnetic field com...
    Alt_logical_source:          ['Ulysses_VHM_Magnetic-Field_1sec_cdf']
    spase_DatasetResourceID:     ['spase://NASA/NumericalData/Ulysses/FGM/PT0...
    CDAWEB_PARENTS:              ['uy_1sec_vhm_00000000_v01', 'uy_1sec_vhm_20...
    PARENTS:                     ['uy_1sec_vhm_00000000_v01', 'uy_1sec_vhm_20...}


start1 = datetime(2004, 2,14,12, 0, 0)
end1   = datetime(2004, 2,14,15, 0, 0)

start2 = datetime(2004, 2,14,21, 0, 0)
end2   = datetime(2004, 2,15, 0, 0, 0)

fig,ax=plt.subplots(figsize=(20,5))
ax.plot(data['UY_1SEC_VHM']['time'],data['UY_1SEC_VHM']['B_RTN'][:,0])
ax.plot(data['UY_1SEC_VHM']['time'],data['UY_1SEC_VHM']['B_RTN'][:,1])
ax.plot(data['UY_1SEC_VHM']['time'],data['UY_1SEC_VHM']['B_RTN'][:,2])
ax.set_xlim('2004-02-14T06:00:00Z','2004-02-15T06:00:00Z')
ax.set_ylabel('magnetic fields [nT]')
ax.set_ylim(-3,3)

### add a rectangle patch ###
import matplotlib.dates as mdates
import matplotlib.patches as patches
startplt = mdates.date2num(start1)
endplt   = mdates.date2num(end1)
width = endplt - startplt
rect = patches.Rectangle((startplt, 0), width, 1, color='grey',alpha=.25, clip_on=False,transform=ax.get_xaxis_transform())
#rect = plt.Rectangle((start,0), width, height=12, clip_on=False, edgecolor="none", facecolor="blue", linewidth=3,alpha=0.25)
ax.add_patch(rect)

startplt = mdates.date2num(start2)
endplt   = mdates.date2num(end2)
width = endplt - startplt
rect = patches.Rectangle((startplt, 0), width, 1, color='grey',alpha=.25, clip_on=False,transform=ax.get_xaxis_transform())
#rect = plt.Rectangle((start,0), width, height=12, clip_on=False, edgecolor="none", facecolor="blue", linewidth=3,alpha=0.25)
ax.add_patch(rect)
#############################

plt.show()


import scipy.fft as sf

def PSD1(start, end):
    ### parameters ###
    nt=np.size(data['UY_1SEC_VHM']['time'][:].sel(time=slice(start,end)).values)
    t=np.arange(nt)
    dt=1    #[s]
    lt=dt*nt
    ww = sf.fftshift(sf.fftfreq(nt, d=dt))
    
    ### window function ###
    hun=(0.50-0.50*np.cos(2*np.pi*t/lt))  #hunn window
    hum=(0.54-0.46*np.cos(2*np.pi*t/lt))  #humming window
    wel= 1 - ((t - 0.5 * nt) / (0.5 * nt)) ** 2          #Welch window
    
    ### normalized power spectrum for each components and trace ###
    br=data['UY_1SEC_VHM']['B_RTN'][:,0].sel(time=slice(start, end)).values
    bt=data['UY_1SEC_VHM']['B_RTN'][:,1].sel(time=slice(start, end)).values
    bn=data['UY_1SEC_VHM']['B_RTN'][:,2].sel(time=slice(start, end)).values
    fftbrbr=sf.fftshift(np.abs(sf.fft((br-np.mean(br))*hun))**2)#/nt**2
    fftbtbt=sf.fftshift(np.abs(sf.fft((bt-np.mean(bt))*hun))**2)#/nt**2
    fftbnbn=sf.fftshift(np.abs(sf.fft((bn-np.mean(bn))*hun))**2)#/nt**2
    ffttrc=fftbrbr+fftbtbt+fftbnbn
    
    ### power spectrum density ###
    psd =ffttrc#*lt
    
    ### smoothing (boxcar method) ###
    nsmooth=10  #smoothing parameter 
    weights = np.ones(nsmooth, dtype=float) / nsmooth
    psd_s   = np.convolve(psd, weights, mode='valid')
    ww_s    = np.convolve(ww , weights, mode='valid')
    
    return ww_s, psd_s

ww_su,psd_su=PSD1(start1,end1)
ww_sd,psd_sd=PSD1(start2,end2)

### plotting a figure ###
plt.figure(figsize=(12,8))
#plt.loglog(ww  ,psd  ,'-', label='without smoothing')
plt.loglog(ww_su,psd_su,'-', label='upstream')
plt.loglog(ww_sd,psd_sd,'-', label='downstream')
plt.loglog(ww_su[ww_su>0],1e0*ww_su[ww_su>0]**(-5./3.),'k--', label='$k^{-5/3}$')
plt.legend()
plt.xlabel('frequency [s$^{-1}$]')
plt.ylabel('PSD')
plt.xlim(1e-5,1e0)
plt.ylim(1e-1,1e7)
#plt.yticks(np.logspace(-6,2,5))
plt.grid(True,which='both' ,axis='x')
plt.grid(True,which='major',axis='y')
plt.show()


def PSD2(start, end):
    ### parameters ###
    nt=np.size(data['UY_1SEC_VHM']['time'][:].sel(time=slice(start,end)).values)
    t=np.arange(nt)
    dt=1    #[s]
    lt=dt*nt
    ww = sf.fftshift(sf.fftfreq(nt, d=dt))
    
    ### window function ###
    hun=(0.50-0.50*np.cos(2*np.pi*t/lt))  #hunn window
    hum=(0.54-0.46*np.cos(2*np.pi*t/lt))  #humming window
    wel= 1 - ((t - 0.5 * nt) / (0.5 * nt)) ** 2          #Welch window
    
    ### autocorrelation function ###
    br=data['UY_1SEC_VHM']['B_RTN'][:,0].sel(time=slice(start, end)).values
    bt=data['UY_1SEC_VHM']['B_RTN'][:,1].sel(time=slice(start, end)).values
    bn=data['UY_1SEC_VHM']['B_RTN'][:,2].sel(time=slice(start, end)).values
    brbr=np.correlate((br-np.mean(br))*hun,(br-np.mean(br))*hun,'same')#/nt
    btbt=np.correlate((bt-np.mean(bt))*hun,(bt-np.mean(bt))*hun,'same')#/nt
    bnbn=np.correlate((bn-np.mean(bn))*hun,(bn-np.mean(bn))*hun,'same')#/nt
    plt.figure(figsize=(12,4))
    plt.plot(brbr,label='$B_RB^\prime_R$')
    plt.plot(btbt,label='$B_TB^\prime_T$')
    plt.plot(bnbn,label='$B_NB^\prime_N$')
    plt.legend()
    plt.xlabel('lag [s]')
    plt.ylabel('autocorrelation function')
    plt.show()

    ### powerspectrum density ###
    trc=(brbr[0:nt]+btbt[0:nt]+bnbn[0:nt])
    psd=np.abs(sf.fftshift(sf.fft(trc)))
    
    ### smoothing (boxcar method) ###
    nsmooth=10  #smoothing parameter 
    weights = np.ones(nsmooth, dtype=float) / nsmooth
    psd_s   = np.convolve(psd, weights, mode='valid')
    ww_s    = np.convolve(ww , weights, mode='valid')
    
    return ww_s, psd_s

ww_su,psd_su=PSD2(start1,end1)
ww_sd,psd_sd=PSD2(start2,end2)

### plotting a figure ###
plt.figure(figsize=(12,8))
#plt.loglog(ww  ,psd  ,'-', label='without smoothing')
plt.loglog(ww_su,psd_su,'-', label='upstream')
plt.loglog(ww_sd,psd_sd,'-', label='downstream')
plt.loglog(ww_su[ww_su>0],1e0*ww_su[ww_su>0]**(-5./3.),'k--', label='$k^{-5/3}$')
plt.legend()
plt.xlabel('frequency [s$^{-1}$]')
plt.ylabel('PSD')
plt.xlim(1e-5,1e0)
plt.ylim(1e-1,1e7)
#plt.yticks(np.logspace(-6,2,5))
plt.grid(True,which='both' ,axis='x')
plt.grid(True,which='major',axis='y')
plt.show()

Power Spectral Density (PSD)¶

PSD: First Method¶

PSD: Second Method¶