diff --git a/Project.toml b/Project.toml index e1484fd7..57a5a0b6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "HiQGA" uuid = "01574e87-63b6-4466-925a-334cf7e21b6b" authors = ["richardt94 <29562583+richardt94@users.noreply.github.com>"] -version = "0.3.12" +version = "0.3.13" [deps] Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" diff --git a/examples/VTEM/DarlingParoo/McMC/01_read_data.jl b/examples/VTEM/DarlingParoo/McMC/01_read_data.jl index d784a1fa..d57f8edd 100644 --- a/examples/VTEM/DarlingParoo/McMC/01_read_data.jl +++ b/examples/VTEM/DarlingParoo/McMC/01_read_data.jl @@ -15,6 +15,7 @@ d = [177,221] soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, fid, linenum, frame_height, d, fname_dat, fname_specs_halt, + tx_rx_dz_pass_through = 0.01, # GA-AEM Z up, +ve is rx over tx startfrom = 1, skipevery = 1, dotillsounding = 2, diff --git a/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl b/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl index f6711da6..d0377e5f 100644 --- a/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl +++ b/examples/VTEM/DarlingParoo/gradientbased/01_read_data.jl @@ -15,6 +15,7 @@ d = [177,221] soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, fid, linenum, frame_height, d, fname_dat, fname_specs_halt, + tx_rx_dz_pass_through = 0.01, # GA-AEM Z up, +ve is rx over tx startfrom = 1, skipevery = 50, dotillsounding = nothing, diff --git a/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl b/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl index 7b9a00c2..80accdd8 100644 --- a/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl +++ b/examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl @@ -50,5 +50,5 @@ times = vec(10 .^mean(log10.([ 0.0100851 0.0113498 ]), dims=2)) include("hedgehog.jl") - +lowpassfcs = [] σ_halt = vec([0.238388 0.175366 0.126135 0.091645 0.071083 0.043692 0.030069 0.019933 0.015175 0.010694 0.008489 0.006127 0.005354 0.005058 0.003871 0.003524 0.002759 0.002185 0.001792 0.001505 0.001298 0.001112 0.000946 0.000805 0.000689 0.000620 0.000546 0.000465 0.000428 0.000399 0.000352 0.000333 0.000317 0.000293 0.000302 0.000334 0.000406 0.000465 0.000458 0.000356 0.000292 0.000193 0.000285 0.000195 0.000177]) \ No newline at end of file diff --git a/examples/VTEM/synth/forwardtests/02_compare_ross.jl b/examples/VTEM/synth/forwardtests/02_compare_ross.jl index bac90612..079efc49 100644 --- a/examples/VTEM/synth/forwardtests/02_compare_ross.jl +++ b/examples/VTEM/synth/forwardtests/02_compare_ross.jl @@ -13,11 +13,11 @@ z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed) rRx = 0.0 zRx = -30.0 zTx = -30.0 -freqlow = 1e-3 +freqlow = 1e-4 freqhigh = 1e6 nkᵣeval = 50 ntimesperdecade = 10 -nfreqsperdecade = 6 +nfreqsperdecade = 5 calcjacobian = true include("../waveletapprox/electronics_halt.jl") ## LM operator diff --git a/examples/VTEM/synth/gradientbased/02_set_options.jl b/examples/VTEM/synth/gradientbased/02_set_options.jl index 1c7aadaa..c34b649b 100644 --- a/examples/VTEM/synth/gradientbased/02_set_options.jl +++ b/examples/VTEM/synth/gradientbased/02_set_options.jl @@ -7,4 +7,4 @@ regtype = :R2 # :R1 is first difference regularization, :R2 is two succesive fir lo, hi = -3, 1 # min, max log10 conductivity λ²min, λ²max=-0.5, 7 # min, max of regularization parameter (Tikhonov parameter) ntries = 8 # how many max tries between λ²min, λ²max -β² = 0. # what fraction of roughness penalty to use to enforce a penalty in deviations from background \ No newline at end of file +β² = 0.001 # what fraction of roughness penalty to use to enforce a penalty in deviations from background \ No newline at end of file diff --git a/examples/VTEM/synth/waveletapprox/electronics_halt.jl b/examples/VTEM/synth/waveletapprox/electronics_halt.jl index 7b9a00c2..80accdd8 100644 --- a/examples/VTEM/synth/waveletapprox/electronics_halt.jl +++ b/examples/VTEM/synth/waveletapprox/electronics_halt.jl @@ -50,5 +50,5 @@ times = vec(10 .^mean(log10.([ 0.0100851 0.0113498 ]), dims=2)) include("hedgehog.jl") - +lowpassfcs = [] σ_halt = vec([0.238388 0.175366 0.126135 0.091645 0.071083 0.043692 0.030069 0.019933 0.015175 0.010694 0.008489 0.006127 0.005354 0.005058 0.003871 0.003524 0.002759 0.002185 0.001792 0.001505 0.001298 0.001112 0.000946 0.000805 0.000689 0.000620 0.000546 0.000465 0.000428 0.000399 0.000352 0.000333 0.000317 0.000293 0.000302 0.000334 0.000406 0.000465 0.000458 0.000356 0.000292 0.000193 0.000285 0.000195 0.000177]) \ No newline at end of file diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index d72f088f..3d229fbf 100644 --- a/src/AEM_VMD_HMD.jl +++ b/src/AEM_VMD_HMD.jl @@ -21,6 +21,7 @@ mutable struct HFieldDHT <: HField ramp :: Array{Float64, 2} log10ω :: Array{Float64, 1} interptimes :: Array{Float64, 1} + minresptime HFD_z :: Array{ComplexF64, 1} HFD_r :: Array{ComplexF64, 1} HFD_az :: Array{ComplexF64, 1} @@ -79,7 +80,7 @@ function HFieldDHT(; zTx = -35.0, zRx = -37.5, times = 10 .^LinRange(-6, -1, 50), - ramp = ones(10, 10), + ramp = ones(10, 2), nfreqsperdecade = 5, ntimesperdecade = 10, glegintegorder = 5, @@ -92,6 +93,7 @@ function HFieldDHT(; getazimH = false, freqlow = 1e-4, freqhigh = 1e6, + minresptime = 1.e-6, # I think responses earlier than this are unstable calcjacobian = false ) @assert all(freqs .> 0.) @@ -103,20 +105,29 @@ function HFieldDHT(; ϵᵢ = similar(pz) rTE = zeros(length(pz)-1) rTM = similar(rTE) - if freqhigh < 3/minimum(times) - freqhigh = 3/minimum(times) + + ## all this below for pesky zero time at start instead of shutoff ... + mintime, maxtime = extrema(times) + @assert mintime > 0 # else interpolation in log10 trouble + mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for + maxtime = 10^(log10(maxtime) + 1) # go ahead a decade further + if doconvramp + mintime, maxtime = checkrampformintime(times, ramp, minresptime, maxtime) + end + interptimes = 10 .^(log10(mintime) : 1/ntimesperdecade : log10(maxtime)) + if freqhigh < 3/mintime + freqhigh = 3/mintime end if freqlow > 0.33/maximum(times) freqlow = 0.33/maximum(times) end if isempty(freqs) - freqs = 10 .^(log10(freqlow):1/nfreqsperdecade:log10(freqhigh)) + freqs = 10 .^(log10(freqlow):1/nfreqsperdecade:log10(freqhigh)) end interpkᵣ = 10 .^range(minimum(log10.(Filter_base))-0.5, maximum(log10.(Filter_base))+0.5, length = nkᵣeval) J0_kernel_h, J1_kernel_h, J0_kernel_v, J1_kernel_v = map(x->zeros(ComplexF64, length(interpkᵣ), length(freqs)), 1:4) J01kernelhold = zeros(ComplexF64, length(Filter_base)) - log10ω = log10.(2*pi*freqs) - interptimes = 10 .^(minimum(log10.(times))-1:1/ntimesperdecade:maximum(log10.(times))+1) + log10ω = log10.(2*pi*freqs) HFD_z = zeros(ComplexF64, length(freqs)) # space domain fields in freq HFD_r = zeros(ComplexF64, length(freqs)) # space domain fields in freq HFD_az = zeros(ComplexF64, length(freqs)) # space domain fields in freq @@ -161,7 +172,7 @@ function HFieldDHT(; # Jtemp = calcjacobian ? zeros(ComplexF64, nmax) : zeros(0) useprimary = modelprimary ? one(Float64) : zero(Float64) - HFieldDHT(thickness, pz, ϵᵢ, zintfc, rTE, rTM, zRx, zTx, rTx, rRx, freqs, times, ramp, log10ω, interptimes, + HFieldDHT(thickness, pz, ϵᵢ, zintfc, rTE, rTM, zRx, zTx, rTx, rRx, freqs, times, ramp, log10ω, interptimes, minresptime, HFD_z, HFD_r, HFD_az, HFD_z_interp, HFD_r_interp, HFD_az_interp, HTD_z_interp, HTD_r_interp, HTD_az_interp, dBzdt, dBrdt, dBazdt, J0_kernel_h, J1_kernel_h, J0_kernel_v, J1_kernel_v, J01kernelhold, lowpassfcs, quadnodes, quadweights, preallocate_ω_Hsc(interptimes, lowpassfcs)..., rxwithinloop, provideddt, doconvramp, useprimary, @@ -170,6 +181,40 @@ function HFieldDHT(; HFD_r_J, HTD_r_J_interp, dBrdt_J) end +function checkrampformintime(times, ramp, minresptime, maxtime) + # this checks we don't have ultra small tobs - ramp_time_a + minta = Inf + maxta = -Inf + for itime = 1:length(times) + for iramp = 1:size(ramp,1)-1 + rta, rtb = ramp[iramp,1], ramp[iramp+1,1] + if rta >= times[itime] # geq instead of gt as we could have an unlcky time + break + end + if rtb > times[itime] # end in this interval + rtb = times[itime] + end + # just so we know, rta < rtb and rta < t so rta < rtb <= t + ta = times[itime]-rta + tb = max(times[itime]-rtb, minresptime) # rtb > rta, so make sure this is not zero because integ is in log10... + @assert ta>tb # else we're in trouble + if ta < minta + minta = ta + end + if ta > maxta + maxta = ta + end + end + end + mintime = max(0.5minta, minresptime) + # though I believe mintime = minta always is safe. + # if maxta > maxtime + maxtime = min(1.5maxta, maxtime) + # 1.5 to make sure we clear maxta in the interptimes + # end + mintime, maxtime +end + #update geometry and dependent parameters - necessary for adjusting geometry #e.g. inverting nuisance parameters function update_ZR!(F::HFieldDHT, zTx, zRx, rTx, rRx; onlyprimary=false) @@ -623,7 +668,7 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu dI = F.ramp[iramp+1,2] - F.ramp[iramp,2] dIdt = dI/dt - if rta > F.times[itime] + if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time break end if rtb > F.times[itime] # end in this interval @@ -631,7 +676,7 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu end ta = F.times[itime]-rta - tb = max(F.times[itime]-rtb, 1e-8) # rtb > rta, so make sure this is not zero... + tb = max(F.times[itime]-rtb, F.minresptime)# rtb > rta, so make sure this is not zero because integ is in log10... a, b = log10(ta), log10(tb) x, w = F.quadnodes, F.quadweights F.dBzdt[itime] += (b-a)/2*dot(getrampresponse((b-a)/2*x .+ (a+b)/2, splz), w)*dIdt diff --git a/src/GradientInversion.jl b/src/GradientInversion.jl index a7edb481..087e10d8 100644 --- a/src/GradientInversion.jl +++ b/src/GradientInversion.jl @@ -261,7 +261,12 @@ function gradientinv( m::AbstractVector, istep += 1 istep > nstepsmax && break end - isa(io, Nothing) || begin @info "Finished "*fname; close(io) end + if !isnothing(io) + @info "Finished "*fname + close(io) + nothing + return + end return map(x->x[1:istep], (mnew, χ², λ², oidx)) end @@ -334,7 +339,12 @@ function gradientinv( m::AbstractVector, istep += 1 istep > nstepsmax && break end - isa(io, Nothing) || begin @info "Finished "*fname; close(io) end + if !isnothing(io) + @info "Finished "*fname + close(io) + nothing + return + end # nuisances mess up occam so this is a kludge as sometimes a smoothest root is not found (istep > nstepsmax) && (istep -= 1) return map(x->x[1:istep], (mnew, nunew, χ², χ²nu, λ², oidx)) diff --git a/src/TEMPEST1DInversion.jl b/src/TEMPEST1DInversion.jl index bce1e390..2341b9dd 100644 --- a/src/TEMPEST1DInversion.jl +++ b/src/TEMPEST1DInversion.jl @@ -111,7 +111,7 @@ function Bfield(; F = AEM_VMD_HMD.HFieldDHT(;nmax, calcjacobian, ntimesperdecade = ntimesperdecade, nfreqsperdecade = nfreqsperdecade, - freqlow=1e-5, + # freqlow =1e-5, nkᵣeval = nkᵣeval, times = times, ramp = ramp, diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl index ec73ae75..acce20af 100644 --- a/src/VTEM1DInversion.jl +++ b/src/VTEM1DInversion.jl @@ -34,6 +34,7 @@ function dBzdt(;times = [1.], ramp = [0. 1], rTx = 10., zTx = -30., + zRx = -30.01, # now Z down for use useML = false, z = [-1.], ρ = [-1.], @@ -46,6 +47,9 @@ function dBzdt(;times = [1.], σ = zeros(0), showgates = false, modelprimary = false, + lowpassfcs = [], + freqlow = 1e-4, + freqhigh = 1e6, ) @assert size(σ) == size(d) @@ -53,7 +57,8 @@ function dBzdt(;times = [1.], F = AEM_VMD_HMD.HFieldDHT(; times, ramp, nmax, zTx, rTx, - rRx = 0., zRx = zTx-0.01, + rRx = 0., zRx, + lowpassfcs, freqlow, freqhigh, calcjacobian, nfreqsperdecade, ntimesperdecade, modelprimary ) @@ -100,6 +105,7 @@ mutable struct VTEMsoundingData <: Sounding fid :: Float64 linenum :: Int zTx :: Float64 + zRx :: Float64 rTx :: Float64 lowpassfcs :: Array{Float64, 1} times :: Array{Float64, 1} @@ -109,14 +115,15 @@ mutable struct VTEMsoundingData <: Sounding end returnforwrite(s::VTEMsoundingData) = [s.X, s.Y, s.Z, s.fid, - s.linenum, s.zTx, s.rTx] + s.linenum, s.zTx, s.zRx, s.rTx] function getndata(S::VTEMsoundingData) getndata(S.data)[1] end function VTEMsoundingData(;rRx=nothing, zRx=nothing, zTx=12., - rTx=-12.,lowpassfcs=[], + rTx=-12.,lowpassfcs=[], + tx_rx_dz = -0.01, # not reading but using format Z down -ve is rx above tx times=[1., 2.], ramp=[1 2; 3 4], noise=[1.], data=[1.], sounding_string="sounding", X=nothing, Y=nothing, Z=nothing, @@ -124,13 +131,13 @@ function VTEMsoundingData(;rRx=nothing, zRx=nothing, zTx=12., @assert rTx > 0 @assert zTx < 0 # my coordinate system z down isnothing(rRx) && (rRx = 0.) - isnothing(zRx) && (zRx = zTx-0.1) # place receiver just above tx centre + isnothing(zRx) && (zRx = zTx + tx_rx_dz) !isempty(lowpassfcs) && @assert all(lowpassfcs .> 0) @assert all(diff(times) .>0 ) @assert all(diff(ramp[:,1]) .>0 ) @assert all((noise .>0) .| isnan.(noise)) @assert length(data) == length(noise) - VTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, zTx, rTx, + VTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, zTx, zRx, rTx, lowpassfcs, times, ramp, noise, data) end @@ -152,6 +159,7 @@ function read_survey_files(; Z = -1, fid = -1, linenum = -1, + tx_rx_dz_pass_through = 0.01, # Z up GA-AEM reading convention +ve is rx above tx nanchar = "*") @assert frame_height > 0 @@ -171,8 +179,10 @@ function read_survey_files(; fiducial = soundings[:,fid] whichline = soundings[:,linenum] d = soundings[:,d[1]:d[2]] - zTx = -soundings[:,frame_height] # my coordinate system - + zTx = soundings[:,frame_height] # read in Z up + zRx = -(zTx .+ tx_rx_dz_pass_through) # my coordinate system Z down + zTx = -zTx # my coordinate system Z down + @info "reading $fname_specs_halt" include(fname_specs_halt) @assert size(d, 2) == length(times) @@ -180,15 +190,14 @@ function read_survey_files(; σ_halt[:] .*= units d[:] .*= units σ = sqrt.((multnoise*d).^2 .+ (σ_halt').^2) - - makeqcplots && plotsoundingdata(d, σ, times, zTx; figsize, fontsize) + makeqcplots && plotsoundingdata(d, σ, times, zTx, zRx; figsize, fontsize) nsoundings = size(d, 1) s_array = Array{VTEMsoundingData, 1}(undef, nsoundings) fracdone = 0 for is in 1:nsoundings l, f = Int(whichline[is]), fiducial[is] - s_array[is] = VTEMsoundingData(;zTx=zTx[is], rTx, - times, ramp, noise=σ[is,:], data=d[is,:], + s_array[is] = VTEMsoundingData(;zTx=zTx[is], zRx=zRx[is], rTx, + times, ramp, noise=σ[is,:], data=d[is,:], lowpassfcs, sounding_string="sounding_$(l)_$f", X=easting[is], Y=northing[is], Z=topo[is], fid=f, linenum=l) @@ -201,7 +210,7 @@ function read_survey_files(; return s_array end -function plotsoundingdata(d, σ, times, zTx; figsize=(8,4), fontsize=1) +function plotsoundingdata(d, σ, times, zTx, zRx; figsize=(8,4), fontsize=1) f, ax = plt.subplots(2, 2, figsize=figsize, gridspec_kw=Dict("width_ratios" => [1,0.01])) nsoundings = size(d, 1) plot_d = permutedims(d) @@ -216,9 +225,11 @@ function plotsoundingdata(d, σ, times, zTx; figsize=(8,4), fontsize=1) axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "r") axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "--w") axx.set_xlabel("avg noise fraction") + ax[2].plot(1:nsoundings, zRx, label="zRx") ax[2].plot(1:nsoundings, zTx, label="zTx") ax[2].set_xlabel("sounding #") ax[2].set_ylabel("height m") + ax[2].invert_yaxis() ax[2].sharex(ax[1]) ax[2].legend() ax[2,2].axis("off") @@ -378,8 +389,8 @@ function makeoperator(sounding::VTEMsoundingData; z, ρ, = makezρ(zboundaries; zfixed, ρfixed) ρ[z.>=zstart] .= ρbg aem = dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, - times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, - rTx=sounding.rTx, zTx=sounding.zTx, z, ρ, calcjacobian, useML, showgates=plotfield) + times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, lowpassfcs=sounding.lowpassfcs, + rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, z, ρ, calcjacobian, useML, showgates=plotfield) plotfield && plotmodelfield!(aem, log10.(ρ[2:end])) aem, zall, znall, zboundaries end @@ -388,9 +399,10 @@ function makeoperator(aem::dBzdt, sounding::VTEMsoundingData) ntimesperdecade = gettimesperdec(aem.F.interptimes) nfreqsperdecade = gettimesperdec(aem.F.freqs) modelprimary = aem.F.useprimary === 1. ? true : false - dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, + dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, lowpassfcs=sounding.lowpassfcs, times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, - rTx=sounding.rTx, zTx=sounding.zTx, z=copy(aem.z), ρ=copy(aem.ρ), + rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, + z=copy(aem.z), ρ=copy(aem.ρ), aem.F.calcjacobian, aem.useML, showgates=false) end diff --git a/test/verify_loop_td.jl b/test/verify_loop_td.jl index 4b459405..fe34a74f 100644 --- a/test/verify_loop_td.jl +++ b/test/verify_loop_td.jl @@ -5,15 +5,15 @@ rho = [1e20, 1, ] nmax = 200 ## geometry rTx = 15 -rRx = 0.001 -zRx = -0.001 -zTx = -0.002 +rRx = 0.0 +zRx = -0.01 +zTx = -0.005 doconvramp = false modelprimary = true nkᵣeval = 200 ntimesperdecade = 15 nfreqsperdecade = 15 -times = 10 .^LinRange(-6,1,40) +times = 10 .^LinRange(-5,0,5) ## build operator F = transD_GP.AEM_VMD_HMD.HFieldDHT( ntimesperdecade = ntimesperdecade, diff --git a/test/verify_skytem_td.jl b/test/verify_skytem_td.jl index 268fbb39..391f79ec 100644 --- a/test/verify_skytem_td.jl +++ b/test/verify_skytem_td.jl @@ -3,7 +3,7 @@ include("skytem_response.jl") ## modeling parameters ntimesperdecade = 10 nfreqsperdecade = 5 -freqlow = 1e-3 +freqlow = 1e-4 ## LM operator Flm = transD_GP.AEM_VMD_HMD.HFieldDHT( ntimesperdecade = ntimesperdecade, @@ -16,7 +16,7 @@ Flm = transD_GP.AEM_VMD_HMD.HFieldDHT( rTx = rTx, zRx = zRxLM, freqlow = freqlow - ) + ); ## HM operator Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT( ntimesperdecade = ntimesperdecade, diff --git a/test/verify_tempest_geometry.jl b/test/verify_tempest_geometry.jl index 8ecb2113..9c56cde8 100644 --- a/test/verify_tempest_geometry.jl +++ b/test/verify_tempest_geometry.jl @@ -8,7 +8,7 @@ nmax = 100 zstart = 0. extendfrac, dz = 1.06, 2. -zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40, showplot=true) +zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40, showplot=false) z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed) @info z @@ -66,7 +66,7 @@ tempest = transD_GP.TEMPEST1DInversion.Bfield( addprimary = true #this ensures that the geometry update actually changes everything that needs to be ) # plot before testing -transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest,z,ρ) +transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest,log10.(ρ[2:end])) μ = transD_GP.TEMPEST1DInversion.μ₀ Bx, By, Bz = tempest.Hx[1]*μ*1e15, tempest.Hy[1]*μ*1e15, tempest.Hz[1]*μ*1e15 diff --git a/test/verify_vmd_loop_fd.jl b/test/verify_vmd_loop_fd.jl index 33dbc310..54a471ed 100644 --- a/test/verify_vmd_loop_fd.jl +++ b/test/verify_vmd_loop_fd.jl @@ -13,12 +13,14 @@ rRx = 100. zRx = -0.011 zTx = -0.01 nkᵣeval = 200 +doconvramp = false ## VMD modelprimary = true -Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT( +Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = zTx, rRx = rRx, + doconvramp, freqs = freqs, zRx = zRx, nkᵣeval = nkᵣeval, @@ -26,12 +28,13 @@ Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT( transD_GP.AEM_VMD_HMD.getfieldFD!(Fvmd, zfixed, rho) ## now use a tiny loop radius - VMD approximation is worse as radius gets larger rTx = 0.001 -Floop = transD_GP.AEM_VMD_HMD.HFieldDHT( +Floop = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = zTx, rRx = rRx, rTx = rTx, freqs = freqs, + doconvramp, zRx = zRx, nkᵣeval = nkᵣeval, modelprimary = modelprimary) @@ -46,12 +49,13 @@ ax[1].legend() ax[1].set_title("Compare with W&H Fig 4.2") ## Compare H radial getradialH = true -Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT( +Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = 0, rRx = rRx, freqs = freqs, zRx = 0, + doconvramp, nkᵣeval = nkᵣeval, modelprimary = false, # says W&H for within plane of observation !! getradialH = getradialH) @@ -70,13 +74,14 @@ zTx = -0.01 zRx = -0.02 rRx = 0.001 rTx = 50.0 -Floopin = transD_GP.AEM_VMD_HMD.HFieldDHT( +Floopin = transD_GP.AEM_VMD_HMD.HFieldDHT(; nmax = nmax, zTx = zTx, rRx = rRx, rTx = rTx, freqs = freqs, zRx = zRx, + doconvramp, nkᵣeval = nkᵣeval, modelprimary = modelprimary) transD_GP.AEM_VMD_HMD.getfieldFD!(Floopin, zfixed, rho) diff --git a/test/verify_vmd_td.jl b/test/verify_vmd_td.jl index 7b0cdee0..2a82483f 100644 --- a/test/verify_vmd_td.jl +++ b/test/verify_vmd_td.jl @@ -1,37 +1,38 @@ -using PyPlot, Revise, HiQGA.transD_GP, Random +using PyPlot, HiQGA.transD_GP, Random ## set up -nfreqsperdecade = 20 -ntimesperdecade = 20 +nfreqsperdecade = 10 +ntimesperdecade = 10 # only for interpolation in plots modelprimary = true provideddt = true doconvramp = false -nkᵣeval = 200 -times = 10 .^LinRange(-5,-1.1, 60) -lowpassfcs = [1e6] -freqhigh = 1e6 -freqlow = 1e-8 +nkᵣeval = 120 # for HT +times = 10 .^LinRange(-5,-1.1, 5) +lowpassfcs = [] +freqhigh = 1e5 +freqlow = 1e-7 ## model zfixed = [-1e5, 0] rho = [1e12, 100] nmax = 200 ## geometry rRx = 100. -zRx = 0. +zRx = -0.01 zTx = 0. -modelprimary = true ## -F = transD_GP.AEM_VMD_HMD.HFieldDHT(freqlow = freqlow, - freqhigh = freqhigh, - zTx = zTx, - rRx = rRx, - times = times, - nfreqsperdecade = nfreqsperdecade, - ntimesperdecade = ntimesperdecade, - zRx = zRx, - nkᵣeval = nkᵣeval, - modelprimary = modelprimary, - provideddt = provideddt, - lowpassfcs = lowpassfcs) +F = transD_GP.AEM_VMD_HMD.HFieldDHT(; + freqlow = freqlow, + freqhigh = freqhigh, + zTx = zTx, + rRx = rRx, + times = times, + doconvramp, + nfreqsperdecade = nfreqsperdecade, + ntimesperdecade = ntimesperdecade, + zRx = zRx, + nkᵣeval = nkᵣeval, + modelprimary = modelprimary, + provideddt = provideddt, + lowpassfcs = lowpassfcs); transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho) ## plot FD figure() @@ -60,7 +61,7 @@ title("Compare with W&H Fig 4.4") plt.tight_layout() ## Radial fields figure(figsize=(4,7)) -F.useprimary = 0. +# F.useprimary = 0. F.getradialH = true F.provideddt = true transD_GP.AEM_VMD_HMD.getfieldTD!(F, zfixed, rho)