# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## This is an example to run CWI using gPar
#
# This example used a doublet recorded in ILAR
#
# In the working directory, it should have a file with array/station's infromation, such as name, lat, lon. Please see the "ar.ilar" in this example for the element and format
#
# A data folder with name of array's name should be presented in the working directory. In this example, it have a folder named "ILAR" 
#
# In folder "ILAR", it has a txt file of doublet list 'db.ILAR.list'. The name of the file is in format of 'db.NAME.list'. Please looking at the example file for its format.
#
# It also has a folder "Data" to store the waveforms. Waveforms from earthquakes are stored separately in folders named by YYYY.MM.DD.HH.MM.SS. The waveforms are stored in sac format with names of YYYY.MM.DD.HH.NET.STA.CHAN.SAC
#
# For each sac file, the header should have coordinates of the station in STLA and STLO

import gpar

# ## Create array object in gPar
# Parameters: 
#             
#             arrayList: name of array list file
#             mode: "db" or "eq". "db" means dealing with earthquake double
#             beamphase: target phase
#             phase_list" a list of phase would like to save travel time infromation from Taup
#             rstime, retime: beginning and the end of time window to select the reference phase for the initial alignment. 
#                             The final windwo would be [reference phase arrival - rstime, refenrence phase arrival + retime]
#             cstime, cetime: similar to rstime and retime but for target phase
#             winlen: windown length in second for in CWI
#             step: length to move the window each step
#             channel: channel of the seismic data would be used
#             filt:filter parameter to bandpass filter the waveform
#             rphase: str, reference phase
#             domain: 'freq' or 'time', domain of cross-correlation calculation
#             fittype: 'cos' or 'parabola' method to subsample the cross-correlation 
#             threshold: minimum cross-correlation coefficient to select the traces to perform beamforming 
#             cut: minimum traces in an array

ar = gpar.build.createArray(arrayList='ar.ilar',mode='db',beamphase='PKIKP',
                            phase_list=['PKIKP','PP','PcP','PKiKP','ScP','SKP','PKP'],
                            rstime=10, retime=60, cstime=10.0, cetime=3, 
                            winlen=3, step=0.05, channel='SHZ',
                            filt=[1,3,3,True],rphase='PKP',domain='freq',
                            fittype='cos',verb=True,threshold=0.2,cut=7)

# Running CWI for doublets in array ILAR and generate plot

ilar = ar.ARRAY.iloc[0]
for e in ilar.doublet:
    e.beamForming(ilar.geometry, steplen=3,step=0.05,
                  rstime=20, retime=60, 
                  cetime=2,cstime=10,
                  filt=[1,3,3,True],channel='SHZ')
    e.plotBeam()

# Write result to a pickle file

ilar.write('ILAR.coda.pkl')



StarcoderdataJupyterScriptsDedupFiltered