revert to default behaviour

Author: MatthewEai

src/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl

@@ -35,6 +35,7 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(;
         wavelet = WaveletCausalRicker(f=5.0),
         freesurface = false,
         imgcondition = "standard",
+        RTM_weight = 0.5,
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
     # active and passive earth model properties.  The active set is in the model-space
@@ -175,6 +176,7 @@ function JetProp2DAcoIsoDenQ_DEO2_FDTD(;
             nbx_inject = nbx_inject,
             wavelet = wavelet,
             freesurface = freesurface,
+            RTM_weight,
             imgcondition = get(icdict, lowercase(imgcondition), WaveFD.ImagingConditionStandard()),
             nthreads = nthreads,
             reportinterval = reportinterval,
@@ -363,7 +365,11 @@ Defaults for arguments are shown inside square brackets.
 * `imgcondition` ["standard"] Selects the type of imaging condition used. Choose from "standard", "FWI", 
     and "RTM". "FWI" and "RTM" will perform Kz wavenumber filtering prior to the imaging condition
     in order to promote long wavelengths (for FWI), or remove long wavelength backscattered energy (for 
-    RTM). Note the true adjoint only exists for "standard" imaging condition currently.
+    RTM).  "MIX" mixes the FWI imaging condition using the parameter RTM_weight. Note the true adjoint 
+    only exists for "standard" imaging condition currently.
+* `RTM_weight` determines the balance of short wavelengths and long wavelengths in the imaging condition. 
+    A value of 0.0 is equivalent to the FWI imaging condition, a value of 0.5 is equivalent to the standard
+    imaging condtion, and a value of 1.0 is equivalent to the RTM imaging condition. 
 * `nthreads [Sys.CPU_THREADS]` The number of threads to use for OpenMP parallelization of the modeling.
 * `reportinterval [500]` The interval at which information about the propagtion is logged.
 
@@ -727,10 +733,21 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     end
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
+
+        # Temporary model updates if mixed imaging condition used. 
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+        δm_RTM = Dict{String,Array{Float32,2}}()
+        δm_All = Dict{String,Array{Float32,2}}()
+    end
+
     for prop in keys(kwargs[:active_modelset])
         δm_ginsu[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        if kwargs[:imgcondition] ===  WaveFD.ImagingConditionWaveFieldSeparationMIX()
+            δm_RTM[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+            δm_All[prop] = zeros(Float32, nz_ginsu, nx_ginsu)
+        end
     end
-
+    
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # Get receiver interpolation coefficients
@@ -804,10 +821,23 @@ function JopProp2DAcoIsoDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             end
 
             # born accumulation
-            if kwargs[:imgcondition] === WaveFD.ImagingConditionStandard() || kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationFWI() || kwargs[:imgcondition] === WaveFD.ImagingConditionWaveFieldSeparationFWI()
+            if kwargs[:imgcondition] !== WaveFD.ImagingConditionWaveFieldSeparationMIX()
+                # Standard born accumulation for FWI, RTM, or standard imaging condition 
                 cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields)
             else
-                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], kwargs[:imgcondition], δm_ginsu, wavefields,0.0)
+                # Born accumulation for mix imaging condition 
+                imcon =  WaveFD.ImagingConditionStandard()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_All, wavefields)
+
+                imcon = WaveFD.ImagingConditionWaveFieldSeparationRTM()
+                cumtime_im += @elapsed WaveFD.adjointBornAccumulation!(p, kwargs[:modeltype], imcon, δm_RTM, wavefields)
+
+                weightAll = 1.0f0 - kwargs[:RTM_weight]
+                weightShort = 2.0f0*kwargs[:RTM_weight] - 1.0f0
+
+                for prop in keys(kwargs[:active_modelset])
+                    δm_ginsu[prop] = (δm_All[prop]  .* weightAll) .+ (δm_RTM[prop]  .* weightShort)
+                end
             end
         end
     end