ok

.github/workflows/blank.yml

@@ -27,7 +27,7 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
@@ -47,7 +47,7 @@ jobs:
   build_docs:
     runs-on: ubuntu-20.04
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1.9'

.gitignore

@@ -1,3 +1,4 @@
 Manifest.toml
 .vscode
-.DS_Store
+.DS_Store
+RUNS/

Project.toml

@@ -4,13 +4,22 @@ authors = ["tduretz <[email protected]>"]
 version = "0.1.0"
 
 [deps]
+ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MathTeXEngine = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SymbolicNumericIntegration = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

README.md

@@ -1,6 +1,13 @@
-# SeismicQ
-
 <h1> <img src="docs/src/assets/logo.png" alt="SeismicQ.jl" width="50"> SeismicQ </h1>
 
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://tduretz.github.io/SeismicQ/dev/)
 [![Build Status](https://github.com/tduretz/SeismicQ/workflows/CI/badge.svg)](https://github.com/tduretz/SeismicQ/actions)
+
+<h1> <img src="docs/src/assets/visuel_ecocup2.png" alt="SeismicQ.jl" width="500">  </h1>
+
+Result: trace of a 1D elastic wave 
+<h1> <img src="docs/src/assets/traceE.png" alt="SeismicQ.jl" width="500">  </h1>
+
+Result: trace of a 1D visco-elastic wave 
+<h1> <img src="docs/src/assets/traceVE.png" alt="SeismicQ.jl" width="500">  </h1>
+

docs/.gitignore

@@ -1,3 +1,6 @@
 Manifest.toml
-.vscode
-build/
+.vscode/
+.vscode/settings.json
+build/
+runs/
+RUNS/

docs/src/index.md

@@ -2,10 +2,30 @@
 *Tout ce que vous avez toujours voulu savoir sur le Q.*
 ## Package Features
 - Compute source functions
-## Function Documentation
+## Function Documentation: Sources
 ```@docs
 Ricker
+```
+## Function Documentation: Rheology
+For deviatoric rheology, updates take the form of:
+```math  
+\tau = \theta_\mathrm{s} \dot\varepsilon + \chi_\mathrm{s} \tau^\mathrm{old}
+```
+For volumetric rheology, updates take the form of:
+```math  
+P = P^\mathrm{old} + \theta_\mathrm{b} \nabla V + \chi_\mathrm{b} \nabla V^\mathrm{old}
+```
+```@docs
+θs
+χs 
+θb
+χb
+```
+## Function Documentation: Treatment
+```@docs
 Spec
 ComputeQgraph
 Getfreq
+GenAttenuatedRicker
+PlotReceiverGather
 ```

examples/FakeData.jl

@@ -0,0 +1,26 @@
+using SeismicQ
+
+function main()
+
+    # Geophone position [m]
+    listₓ = 0:100:5000;
+
+    # Time domain
+    Δt  = 1e-3
+    Nt  = 2000
+
+    # Central frequency of the source [Hz]
+    𝑓₀  = 10.
+
+    # Velocities
+    Vp = 7000 # m/s
+    Vs = 4000 
+    αp = 2e-4
+    αs = 4e-4 # (π * f)/ (Q * V)
+
+    time_vec,acc_vec=GenAttenuatedRicker(listₓ,Δt,Nt,Vp,Vs,αp,αs,𝑓₀)
+
+    PlotReceiverGather(listₓ,time_vec,acc_vec)
+end
+
+main()

examples/GetQTest.jl

@@ -2,19 +2,30 @@ using SeismicQ,Plots
 
 function main()
 
-#P and S velocities
-Vp=7000
-Vs=4000
-
-SeismicMatrix=GenMatrix(Vp,Vs)
+# Geophone position [m]
+listₓ = 0:100:5000;
 
-#distances from source
+# Time domain
+Δt  = 1e-3
+Nt  = 2000
+
+# Central frequency of the source [Hz]
+𝑓₀  = 10.
+
+# Velocities
+Vp = 7000 # m/s
+Vs = 4000 
+αp = 2e-4
+αs = 4e-4 # (π * f)/ (Q * V)
+
+#Selected distances from source
 d1=3000
 d2=5000
 
-#Number of samples and time increment of traces
-Nt=2000
-dt=1e-3
+TimeVec,AccVec=GenAttenuatedRicker(listₓ,Δt,Nt,Vp,Vs,αp,αs,𝑓₀)
+#SeismicMatrix=GenMatrix(Vp,Vs)
+
+dt=Δt
 
 #min and max frequencies to compute Q
 fqmin=5.0
@@ -28,8 +39,8 @@ println("Trace length is $Tmax s using $Nt samples")
 ind1=Int(trunc(d1/100))
 ind2=Int(trunc(d2/100))
 
-rec1=SeismicMatrix[ind1,2:Nt+1]
-rec2=SeismicMatrix[ind2,2:Nt+1]
+rec1=AccVec[ind1,1:Nt]
+rec2=AccVec[ind2,1:Nt]
 
 println("Trace 1 at offset $d1 m is found at index $ind1")
 println("Trace 2 at offset $d2 m is found at index $ind2")

examples/PlotRicker.jl

@@ -0,0 +1,47 @@
+using SeismicQ, GLMakie, MathTeXEngine
+Makie.update_theme!(fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic)))
+
+function MainRicker2D()
+    # Spatial extent
+    l  = (x = 25, y = 25)
+
+    # Discretization
+    Nc  = (x = 200, y = 200) 
+    Δ   = (x = l.x/Nc.x, y = l.y/Nc.y, z=1.0)
+    X   = (v = (x= LinRange(0,l.x,Nc.x+1)            , y= LinRange(0,l.y,Nc.y+1)),
+           c = (x= LinRange(0-Δ.x/2,l.x+Δ.x/2,Nc.x+2) , y= LinRange(0-Δ.y/2,l.y+Δ.y/2,Nc.y+2)),
+           i = (x= LinRange(0,l.x,Nc.x+1)             , y= LinRange(0-Δ.y/2,l.y+Δ.y/2,Nc.y+2)),
+           j = (x= LinRange(0-Δ.x/2,l.x+Δ.x/2,Nc.x+2) , y= LinRange(0,l.y,Nc.y+1))) 
+
+    # Array sizes
+    szv   = (Nc.x+1, Nc.y+1)
+    szc   = (Nc.x+2, Nc.y+2)
+
+    # Source parameters
+    𝑓₀   = 50   # Central frequency of the source [Hz]
+    t₀   = 1.2/𝑓₀
+    σ₀   = l.x/30
+    x₀   = l.x/2
+    y₀   = l.y/2
+    t    = t₀    
+
+    # Compute Ricker function with 2D spatial support
+    f_ext = (v=zeros(szv)  , c=zeros(szc))
+    x2d   = X.v.x * ones(size( X.v.y))'
+    y2d   = ones(size( X.v.x)) * X.v.y'
+    @. f_ext.v = Ricker.( x2d, x₀, y2d, y₀, t, t₀, 𝑓₀, σ₀)
+
+    # Makie visualisation
+    resol = 1000 
+    f     = Figure(resolution = (l.x/l.y*resol, resol), fontsize=25)
+    ax1   = Axis(f[1, 1], title = L" vz on v grid at $t$ = %$(t) [s]", xlabel = L"$x$ [m]", ylabel = L"$y$ [m]")
+    hm    = GLMakie.heatmap!(ax1, X.v.x, X.v.y, f_ext.v, colormap = :turbo)
+    colsize!(f.layout, 1, Aspect(1, l.x/l.y))
+    GLMakie.Colorbar(f[1, 2], hm, label = "Vz [m/s]", width = 20, labelsize = 25, ticklabelsize = 14 )
+    GLMakie.colgap!(f.layout, 20)
+    # DataInspector(f)
+    display(f)
+
+end
+
+MainRicker2D()