From b526e8694503a1f74b4d9a911d0157ffddcfbe48 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Fri, 2 Aug 2024 09:32:57 +1000
Subject: [PATCH 01/17] WIP: propagate VTEM changes everywhere

---
 src/AEMnoNuisanceGradientInversionTools.jl | 11 ++++++-----
 src/VTEM1DInversion.jl                     |  8 +++++---
 2 files changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl
index 4559b1e..f5908b3 100644
--- a/src/AEMnoNuisanceGradientInversionTools.jl
+++ b/src/AEMnoNuisanceGradientInversionTools.jl
@@ -102,8 +102,9 @@ function plotconvandlast(soundings, delr, delz;
             (calcresiduals && doreshist) && plotgausshist(res[i][:], title="Line $(soundings[a].linenum) residuals")
     end
     getphidhist(ϕd; doplot=dophiplot, saveplot=dophiplot, prefix=prefix)
-    if calcresiduals 
-        rall = reduce(hcat, res)
+    if calcresiduals
+        idxgood = isassigned.(Ref(res), 1:length(res))
+        rall = reduce(hcat, res[idxgood])
         doreshist && plotgausshist(vec(rall), title="All Lines residuals")
         map(eachrow(rall)) do r
             n = length(r)
@@ -154,7 +155,7 @@ function plotconvandlasteachline(soundings, σ, ϕd, delr, delz, resid;
     height_ratios = nextra == 1 ? [1,1,2,4,0.25] : [1,1,4,0.25]
     fig, ax = plt.subplots(4+nextra, 1, gridspec_kw=Dict("height_ratios" => height_ratios),
         figsize=figsize)
-    lname = "Line_$(soundings[1].linenum)"*postfix
+    lname = "Line_$(soundings[1].linenum)_"*postfix
     x0, y0 = soundings[1].X, soundings[1].Y
     if isdefined(soundings[1], :zTx)
         zTx = [s.zTx for s in soundings]
@@ -190,10 +191,10 @@ function plotconvandlasteachline(soundings, σ, ϕd, delr, delz, resid;
     [a.tick_params(labelbottom=false) for a in ax[1:irow-1]]
     isnothing(idx) || plotprofile(ax[irow], idx, Z, R)
     # eg = extrema(gridr)
-    isa(yl, Nothing) || ax[3].set_ylim(yl...)
+    isa(yl, Nothing) || ax[irow].set_ylim(yl...)
     ax[irow].set_ylabel("mAHD")
     ax[irow].set_xlabel("Distance m")
-    ax[irow].sharex(ax[2])
+    ax[irow].sharex(ax[irow-1])
     ax[irow].set_xlim(extrema(gridr))
     plotNEWSlabels(soundings, gridr, gridz, [ax[irow]], x0, y0, xend, yend, 
         preferEright=preferEright, preferNright=preferNright; fontsize)
diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl
index 1e64fd7..7c68d78 100644
--- a/src/VTEM1DInversion.jl
+++ b/src/VTEM1DInversion.jl
@@ -74,11 +74,13 @@ function allocateJ(FJt, σ, select, nfixed, nmodel, calcjacobian)
     if calcjacobian && !isempty(select)
         J = FJt'
         J = J[select,nfixed+1:nmodel]
-        Wdiag = 1 ./σ[select]
-        res = similar(Wdiag)
     else    
-        res, J, Wdiag = zeros(0), zeros(0), zeros(0)
+        J, Wdiag = zeros(0), zeros(0), zeros(0)
     end
+    # always return an allocated residuals and W - small price to pay I think
+    # since majority of Jacobian allocations are in aem.F
+    res = similar(σ[select]) 
+    Wdiag = 1 ./σ[select]
     W = sparse(diagm(Wdiag))        
     return res, J, W
 end

From 270fb2de207f3bc5215df4a411e498e1d17e6c88 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Tue, 27 Aug 2024 12:06:41 +1000
Subject: [PATCH 02/17] seg-y stuff

---
 zz_portalcurtains/RDP.jl                      | 46 ++++++++++++++++++-
 zz_portalcurtains/makesegy_from_all_xyzrho.jl | 15 ++++++
 zz_portalcurtains/makesegy_from_xyzrho.jl     | 11 +++++
 zz_portalcurtains/writesegy.jl                | 11 +++++
 4 files changed, 82 insertions(+), 1 deletion(-)
 create mode 100644 zz_portalcurtains/makesegy_from_all_xyzrho.jl
 create mode 100644 zz_portalcurtains/makesegy_from_xyzrho.jl
 create mode 100644 zz_portalcurtains/writesegy.jl

diff --git a/zz_portalcurtains/RDP.jl b/zz_portalcurtains/RDP.jl
index 58e074d..7e5d1eb 100644
--- a/zz_portalcurtains/RDP.jl
+++ b/zz_portalcurtains/RDP.jl
@@ -1,6 +1,7 @@
 module RDP
 using LinearAlgebra, Dates, ArchGDAL, Printf, DataFrames, CSV,
-    PyPlot, Images, FileIO, HiQGA, Interpolations, NearestNeighbors, DelimitedFiles
+    PyPlot, Images, FileIO, HiQGA, Interpolations, NearestNeighbors, DelimitedFiles,
+    SegyIO
 import GeoFormatTypes as GFT
 import GeoDataFrames as GDF
 import HiQGA.transD_GP.LineRegression.getsmoothline, HiQGA.transD_GP.CommonToAll.colstovtk
@@ -50,6 +51,49 @@ function makeextent(lnum, ncols, nrows, gridr, gridz; suffix="")
         lnum, suffix, 0, 0, ncols, -nrows, gridr[1], maximum(gridz), gridr[end], minimum(gridz))
 end    
 
+function XYZ_zmid_gridtoSEGY(σ, X, Y, Z; dr=nothing, zall=nothing, dz=nothing, fname="line", dst_dir="", suffix="",
+            nanval=-6 #=-6 is in log10=#,)
+    img, gridr, gridz, topofine, R = transD_GP.makegrid(σ, X, Y, Z; 
+            dr, zall, dz)
+    img[isnan.(img)] .= nanval       
+    segypath = dst_dir
+    xymid, = transD_GP.getallxyinr(X, Y, dr)
+    xm, ym, = map(i->xymid[i,:], 1:2)
+    rm = transD_GP.CommonToAll.cumulativelinedist(xm, ym)
+    topom = (interpolate((gridr,), topofine, Gridded(Linear())))(rm)
+    block = SeisBlock(Float32.(img))
+    set_header!(block, :dt, dz*1000) # ms
+    set_header!(block, :dtOrig, dz*1000) # ms
+    set_header!(block, :nsOrig, size(img,1))
+    set_header!(block, :SourceX, round.(Int, xm))
+    set_header!(block, :SourceY, round.(Int, ym))
+    set_header!(block, :GroupX, round.(Int, xm))
+    set_header!(block, :GroupY, round.(Int, ym))
+    set_header!(block, :CDPTrace, Array(0:length(xm)-1))
+    set_header!(block, :Tracenumber, Array(0:length(xm)-1))
+    set_header!(block, :Inline3D, 1)
+    set_header!(block, :Crossline3D, Array(1:length(xm)))
+    set_header!(block, :ElevationScalar, round.(Int, topom))
+    set_header!(block, :DelayRecordingTime, -round(Int, maximum(gridz))) # shift to topo as start
+    fname = joinpath(segypath, fname*"_"*suffix)
+    segy_write(fname*".segy", block)
+end    
+
+function writesegyfromxyzrhodir(nlayers::Int; src_dir="", src_epsg=0, dst_dir="", dr=nothing, dz=nothing)
+    @assert !isnothing(src_epsg)
+    lines = transD_GP.getprobabilisticlinesfromdirectory(src_dir)
+    isdir(dst_dir) || mkpath(dst_dir)
+    ioproj = open(joinpath(dst_dir, "0000_projection.txt"), "w")
+    write(ioproj, "EPSG: $src_epsg")
+    close(ioproj)
+    map(lines) do ln
+        @info "doing line $ln"
+        fname = "line_$ln"
+        X, Y, Z, zall, ρlow, ρmid, ρhigh, ρavg, ϕmean, ϕsdev = transD_GP.readxyzrhoϕ(ln, nlayers; pathname=src_dir)
+        [XYZ_zmid_gridtoSEGY(-ρ, X, Y, Z; dr, zall, dz, dst_dir, fname, suffix=str) for (ρ, str) in zip([ρlow, ρmid, ρhigh],["high", "mid", "low"])]
+    end
+end
+
 function writepathextentfiles(line, nrows, ncols, gridr, topo, gridz, X, Y; suffix="", dst_dir="", src_epsg=0)
     geompath = joinpath(dst_dir,"geometry")
     isdir(geompath) || mkpath(geompath)
diff --git a/zz_portalcurtains/makesegy_from_all_xyzrho.jl b/zz_portalcurtains/makesegy_from_all_xyzrho.jl
new file mode 100644
index 0000000..3a07ac3
--- /dev/null
+++ b/zz_portalcurtains/makesegy_from_all_xyzrho.jl
@@ -0,0 +1,15 @@
+using HiQGA, PyPlot
+includet("RDP.jl")
+nlayers = 52
+dr = 250
+dz = 1
+items = [item for item in walkdir("/Users/anray/Documents/work/projects/largeaem/final_01/summaries_all")]
+idx_summary = [basename(it[1]) == "summary" for it in items]
+src_dir = [it[1] for it in items[idx_summary]]
+src_epsg = [28353, 28354, 28351, 28354, 28354, 28354, 28354, 28352, 28352, 28352, 28352]
+map(zip(src_dir, src_epsg)) do (sdir, epsg)
+    parts = split(sdir,"/")
+    fname = parts[end-2]*"_"*parts[end-1]
+    dst_dir = joinpath("/Users/anray/Documents/work/projects/largeaem/final_01/","segy",fname)
+    RDP.writesegyfromxyzrhodir(nlayers; src_dir=sdir, src_epsg=epsg, dst_dir, dr, dz)
+end
diff --git a/zz_portalcurtains/makesegy_from_xyzrho.jl b/zz_portalcurtains/makesegy_from_xyzrho.jl
new file mode 100644
index 0000000..21e985a
--- /dev/null
+++ b/zz_portalcurtains/makesegy_from_xyzrho.jl
@@ -0,0 +1,11 @@
+using HiQGA, PyPlot
+includet("RDP.jl")
+src_epsg = 28354
+nlayers = 52
+src_dir = "/Users/anray/Documents/work/projects/largeaem/final_01/summaries_all/AusAEM_03_ERC/2021_ERC_01/summary"
+dr = 250
+dz = 1
+parts = split(src_dir,"/")
+fname = parts[end-2]*"_"*parts[end-1]
+dst_dir = joinpath("/Users/anray/Documents/work/projects/largeaem/final_01/","segy",fname)
+RDP.writesegyfromxyzrhodir(nlayers; src_dir, src_epsg, dst_dir, dr, dz)
diff --git a/zz_portalcurtains/writesegy.jl b/zz_portalcurtains/writesegy.jl
new file mode 100644
index 0000000..6e9bcb1
--- /dev/null
+++ b/zz_portalcurtains/writesegy.jl
@@ -0,0 +1,11 @@
+using SegyIO, PyPlot, HiQGA
+cd(@__DIR__)
+includet("RDP.jl")
+pathname="/Users/anray/Documents/work/projects/largeaem/final_01/summaries_all/AusAEM_03_ERC/2021_ERC_01/summary"
+nlayers = 52
+line = 1001001
+X, Y, Z, zall, ρlow, ρmid, ρhigh, = transD_GP.readxyzrhoϕ(line, nlayers; pathname)
+dr = 250.
+dz = 2
+RDP.XYZ_zmid_gridtoSEGY(-ρlow, X, Y, Z; dr, zall, dz)# one line
+

From 069572f1d01ac75fe35386b81d839acc560f9e9e Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Tue, 27 Aug 2024 12:39:12 +1000
Subject: [PATCH 03/17] WIP: resid changes as always return an allocated
 residuals and W since majority of Jacobian allocations are in aem.F

---
 src/AEMwithNuisanceGradientInversionTools.jl |  5 +++--
 src/SkyTEM1DInversion.jl                     | 10 ++++++----
 src/TEMPEST1DInversion.jl                    | 10 ++++++----
 src/VTEM1DInversion.jl                       |  6 +++---
 4 files changed, 18 insertions(+), 13 deletions(-)

diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl
index c02ac83..b74923a 100644
--- a/src/AEMwithNuisanceGradientInversionTools.jl
+++ b/src/AEMwithNuisanceGradientInversionTools.jl
@@ -107,7 +107,8 @@ function plotconvandlast(soundings, delr, delz, nufieldnames::Vector{Symbol};
             (calcresiduals && doreshist) && plotgausshist(res[i][:], title="Line $(soundings[a].linenum) residuals")
     end
     getphidhist(ϕd, doplot=dophiplot, saveplot=dophiplot, prefix=prefix)
-    if calcresiduals 
+    if calcresiduals
+        idxgood = isassigned.(Ref(res), 1:length(res))
         rall = reduce(hcat, res)
         doreshist && plotgausshist(vec(rall), title="All Lines residuals")
         map(eachrow(rall)) do r
@@ -159,7 +160,7 @@ function plotconvandlasteachline(soundings, σ, nu, nufieldnames, ϕd, delr, del
     height_ratios = nextra == 1 ? [1,ones(1+nnu)...,2,4, 0.25] : [1,ones(1+nnu)...,4, 0.25]
     fig, ax = plt.subplots(4+nextra+nnu, 1, gridspec_kw=Dict("height_ratios" => height_ratios),
         figsize=figsize)
-    lname = "Line_$(soundings[1].linenum)"*postfix
+    lname = "Line_$(soundings[1].linenum)_"*postfix
     x0, y0 = soundings[1].X, soundings[1].Y
     if isdefined(soundings[1], :z_tx) # make this a symbol key TODO
         zTx = [s.z_tx for s in soundings]
diff --git a/src/SkyTEM1DInversion.jl b/src/SkyTEM1DInversion.jl
index 5bca4ab..7c7b6e8 100644
--- a/src/SkyTEM1DInversion.jl
+++ b/src/SkyTEM1DInversion.jl
@@ -94,11 +94,13 @@ function allocateJ(Flow, Fhigh, σlow, σhigh, selectlow, selecthigh, nfixed, nm
     if calcjacobian && (!isempty(selectlow) || !isempty(selecthigh))
         J = [Flow.dBzdt_J'; Fhigh.dBzdt_J']; 
         J = J[[selectlow;selecthigh],nfixed+1:nmodel]
-        Wdiag = [1 ./σlow[selectlow]; 1 ./σhigh[selecthigh]]
-        res = similar(Wdiag)
     else    
-        res, J, Wdiag = zeros(0), zeros(0), zeros(0)
-    end    
+        J  = zeros(0)
+    end
+    # always return an allocated residuals and W - small price to pay I think
+    # since majority of Jacobian allocations are in aem.F    
+    Wdiag = [1 ./σlow[selectlow]; 1 ./σhigh[selecthigh]]
+    res = similar(Wdiag)
     W = sparse(diagm(Wdiag))
     return res, J, W
 end    
diff --git a/src/TEMPEST1DInversion.jl b/src/TEMPEST1DInversion.jl
index 35a1a3c..786a106 100644
--- a/src/TEMPEST1DInversion.jl
+++ b/src/TEMPEST1DInversion.jl
@@ -132,14 +132,16 @@ function Bfield(;
 end
 
 function allocateJ(ntimes, nfixed, nmodel, calcjacobian, vectorsum, σx, σz)
+    multfac = vectorsum ? 1 : 2
     if calcjacobian
-        multfac = vectorsum ? 1 : 2
         J = zeros(multfac*ntimes, nmodel-nfixed)
-        Wdiag = vectorsum ? ones(multfac*ntimes) : 1 ./[σx; σz]
-        res = similar(Wdiag)
     else
-        J, Wdiag, res = zeros(0), zeros(0), zeros(0)
+        J = zeros(0)
     end
+    # always return an allocated residuals and W - small price to pay I think
+    # since majority of Jacobian allocations are in aem.F
+    Wdiag = vectorsum ? ones(multfac*ntimes) : 1 ./[σx; σz]
+    res = similar(Wdiag)
     W = sparse(diagm(Wdiag))
     res, J, W    
 end
diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl
index 7c68d78..df5ae28 100644
--- a/src/VTEM1DInversion.jl
+++ b/src/VTEM1DInversion.jl
@@ -75,7 +75,7 @@ function allocateJ(FJt, σ, select, nfixed, nmodel, calcjacobian)
         J = FJt'
         J = J[select,nfixed+1:nmodel]
     else    
-        J, Wdiag = zeros(0), zeros(0), zeros(0)
+        J = zeros(0)
     end
     # always return an allocated residuals and W - small price to pay I think
     # since majority of Jacobian allocations are in aem.F
@@ -239,8 +239,8 @@ function plotsoundingdata(d, σ, times, zTx, zRx; figsize=(8,4), fontsize=1)
     ax[1].set_ylabel("time s")
     ax[1].tick_params(labelbottom=false)
     axx = ax[1].twiny()
-    axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "r")
-    axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "--w")
+    axx.semilogy(infnanmean(σ./abs.(d), 1)[:], times, "r")
+    axx.semilogy(infnanmean(σ./abs.(d), 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")

From 4a4b3f98d559e8a2c138465d02e0745967bf0ae3 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Tue, 27 Aug 2024 13:24:26 +1000
Subject: [PATCH 04/17] showplot modifications

---
 examples/tempest/synth/gradientbased/03_run_inversion.jl | 5 -----
 src/AEMwithNuisanceGradientInversionTools.jl             | 2 +-
 src/SkyTEM1DInversion.jl                                 | 6 ++++--
 src/VTEM1DInversion.jl                                   | 8 ++++++--
 4 files changed, 11 insertions(+), 10 deletions(-)
 delete mode 100644 examples/tempest/synth/gradientbased/03_run_inversion.jl

diff --git a/examples/tempest/synth/gradientbased/03_run_inversion.jl b/examples/tempest/synth/gradientbased/03_run_inversion.jl
deleted file mode 100644
index c317a29..0000000
--- a/examples/tempest/synth/gradientbased/03_run_inversion.jl
+++ /dev/null
@@ -1,5 +0,0 @@
-## do the inversion
-m, χ², λ², idx = transD_GP.gradientinv(σstart, σ0, tempest, nstepsmax=20, 
-                            ;λ²min, λ²max, β², ntries,
-                            lo, hi, regtype);
-aem = tempest;
\ No newline at end of file
diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl
index b74923a..f33c4d6 100644
--- a/src/AEMwithNuisanceGradientInversionTools.jl
+++ b/src/AEMwithNuisanceGradientInversionTools.jl
@@ -109,7 +109,7 @@ function plotconvandlast(soundings, delr, delz, nufieldnames::Vector{Symbol};
     getphidhist(ϕd, doplot=dophiplot, saveplot=dophiplot, prefix=prefix)
     if calcresiduals
         idxgood = isassigned.(Ref(res), 1:length(res))
-        rall = reduce(hcat, res)
+        rall = reduce(hcat, res[idxgood])
         doreshist && plotgausshist(vec(rall), title="All Lines residuals")
         map(eachrow(rall)) do r
             n = length(r)
diff --git a/src/SkyTEM1DInversion.jl b/src/SkyTEM1DInversion.jl
index 7c7b6e8..75b617c 100644
--- a/src/SkyTEM1DInversion.jl
+++ b/src/SkyTEM1DInversion.jl
@@ -432,9 +432,11 @@ function makenoisydata!(aem, ρ;
     # for Gauss-Newton
     aem.res, aem.J, aem.W = allocateJ(aem.Flow, aem.Fhigh, aem.σlow, aem.σhigh, 
                     aem.selectlow, aem.selecthigh, aem.nfixed, length(aem.ρ))
-
     
-    showplot && plotmodelfield!(aem, ρ; onesigma, color, alpha, model_lw, forward_lw, figsize, revax)
+    if showplot
+        plotwaveformgates(aem)
+        plotmodelfield!(aem, ρ; onesigma, color, alpha, model_lw, forward_lw, figsize, revax)
+    end
     nothing
 end
 
diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl
index df5ae28..9d6a316 100644
--- a/src/VTEM1DInversion.jl
+++ b/src/VTEM1DInversion.jl
@@ -366,7 +366,7 @@ end
 
 # noisy synthetic model making
 function makenoisydata!(aem, ρ;
-        rseed=123, noisefrac=0.03, σ_halt=nothing, useML=false,
+        rseed=123, noisefrac=0.03, σ_halt=nothing, useML=false, showplot=true,
         onesigma=true, color=nothing, alpha=1, model_lw=1, forward_lw=1, figsize=(8,6), revax=true,
         # σ_halt default assumed in Bfield units of pV
         units=1/pVinv)
@@ -381,7 +381,11 @@ function makenoisydata!(aem, ρ;
     aem.ndata, aem.select = getndata(aem.d)
     aem.res, aem.J, aem.W = allocateJ(aem.F.dBzdt_J, aem.σ, aem.select, 
         aem.nfixed, length(aem.ρ), aem.F.calcjacobian)
-    plotmodelfield!(aem, ρ; onesigma, color, alpha, model_lw, forward_lw, figsize, revax)
+
+    if showplot
+        plotwaveformgates(aem)
+        plotmodelfield!(aem, ρ; onesigma, color, alpha, model_lw, forward_lw, figsize, revax)
+    end
     nothing
 end
 

From 7b5de9f3b2316615a903063a24a35edd88ee30fc Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Thu, 29 Aug 2024 11:50:35 +1000
Subject: [PATCH 05/17] changes for ta always greater than mintime

---
 src/AEM_VMD_HMD.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index 3d229fb..41fd619 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -196,7 +196,11 @@ function checkrampformintime(times, ramp, minresptime, maxtime)
             end
             # just so we know, rta < rtb and rta < t so rta < rtb <= t
             ta = times[itime]-rta 
+            if ta <= minresptime
+                break # since ta must always be greater than minresptime
+            end
             tb = max(times[itime]-rtb, minresptime) # rtb > rta, so make sure this is not zero because integ is in log10...
+            @info ta, tb, itime, times[itime]
             @assert ta>tb # else we're in trouble
             if ta < minta
                 minta = ta
@@ -676,6 +680,9 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
             end
 
             ta = F.times[itime]-rta
+            if ta <= F.minresptime
+                break # since ta must always be greater than minresptime
+            end
             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

From 42a7a4b3a2027c8fcdd9d7b41965a77f87edbe5d Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Thu, 29 Aug 2024 13:23:51 +1000
Subject: [PATCH 06/17] removed @info

---
 src/AEM_VMD_HMD.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index 41fd619..312868e 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -200,7 +200,6 @@ function checkrampformintime(times, ramp, minresptime, maxtime)
                 break # since ta must always be greater than minresptime
             end
             tb = max(times[itime]-rtb, minresptime) # rtb > rta, so make sure this is not zero because integ is in log10...
-            @info ta, tb, itime, times[itime]
             @assert ta>tb # else we're in trouble
             if ta < minta
                 minta = ta

From c3f7ff35fb942c978e371b49477a0efe4f8a77d0 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Thu, 29 Aug 2024 18:36:30 +1000
Subject: [PATCH 07/17] dIdt Rx waveform added for VTEM like

---
 src/AEM_VMD_HMD.jl     | 19 ++++++++++++-------
 src/VTEM1DInversion.jl | 14 +++++++++-----
 2 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index 312868e..cf926dc 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -69,7 +69,8 @@ mutable struct HFieldDHT <: HField
     dBazdt_J
     HFD_r_J
     HTD_r_J_interp     
-    dBrdt_J    
+    dBrdt_J
+    isdIdt
 end
 
 function HFieldDHT(;
@@ -94,7 +95,8 @@ function HFieldDHT(;
       freqlow = 1e-4,
       freqhigh = 1e6,
       minresptime = 1.e-6, # I think responses earlier than this are unstable
-      calcjacobian = false
+      calcjacobian = false,
+      isdIdt = false
   )
     @assert all(freqs .> 0.)
     @assert freqhigh > freqlow
@@ -178,7 +180,7 @@ function HFieldDHT(;
             quadnodes, quadweights, preallocate_ω_Hsc(interptimes, lowpassfcs)..., rxwithinloop, provideddt, doconvramp, useprimary,
             nkᵣeval, interpkᵣ, log10interpkᵣ, log10Filter_base, getradialH, getazimH, 
             calcjacobian, Jtemp, similar(Jtemp), Jac_z, Jac_az, Jac_r, HFD_z_J, HTD_z_J_interp, dBzdt_J,  HFD_az_J, HTD_az_J_interp, dBazdt_J, 
-            HFD_r_J, HTD_r_J_interp, dBrdt_J)
+            HFD_r_J, HTD_r_J_interp, dBrdt_J, isdIdt)
 end
 
 function checkrampformintime(times, ramp, minresptime, maxtime)
@@ -667,10 +669,13 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
     for itime = 1:length(F.times)
         for iramp = 1:size(F.ramp,1)-1
             rta, rtb  = F.ramp[iramp,1], F.ramp[iramp+1,1]
-            dt   = rtb - rta
-            dI   = F.ramp[iramp+1,2] - F.ramp[iramp,2]
-            dIdt = dI/dt
-
+            if !F.isdIdt
+                dt   = rtb - rta
+                dI   = F.ramp[iramp+1,2] - F.ramp[iramp,2]
+                dIdt = dI/dt
+            else
+                dIdt = F.ramp[iramp+1,2]
+            end
             if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time
                 break
             end
diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl
index 9d6a316..6c7089d 100644
--- a/src/VTEM1DInversion.jl
+++ b/src/VTEM1DInversion.jl
@@ -49,7 +49,8 @@ function dBzdt(;times           = [1.],
         modelprimary    = false,
         lowpassfcs      = [],
         freqlow         = 1e-4,
-        freqhigh        = 1e6, 
+        freqhigh        = 1e6,
+        isdIdt          = false
         )
    
     @assert size(σ)  == size(d)
@@ -60,7 +61,8 @@ function dBzdt(;times           = [1.],
         rRx = 0., zRx, 
         lowpassfcs, freqlow, freqhigh,
         calcjacobian, nfreqsperdecade,
-        ntimesperdecade, modelprimary
+        ntimesperdecade, modelprimary,
+        isdIdt
     )
     @assert length(F.thickness) >= length(z)
     # for Gauss-Newton
@@ -403,7 +405,8 @@ function makeoperator(sounding::VTEMsoundingData;
             nfreqsperdecade = 5,
             showgeomplot  = false,
             calcjacobian  = false,
-            plotfield     = false
+            plotfield     = false,
+            isdIdt        = false
             )
     
     zall, znall, zboundaries = setupz(zstart, extendfrac, dz=dz, n=nlayers, showplot=showgeomplot)
@@ -411,7 +414,7 @@ function makeoperator(sounding::VTEMsoundingData;
     ρ[z.>=zstart] .= ρbg
     aem = dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary,
         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)
+        rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, z, ρ, calcjacobian, useML, showgates=plotfield, isdIdt)
     plotfield && plotmodelfield!(aem, log10.(ρ[2:end]))
     aem, zall, znall, zboundaries
 end
@@ -420,8 +423,9 @@ function makeoperator(aem::dBzdt, sounding::VTEMsoundingData)
     ntimesperdecade = gettimesperdec(aem.F.interptimes)
     nfreqsperdecade = gettimesperdec(aem.F.freqs)
     modelprimary = aem.F.useprimary === 1. ? true : false
+    isdIdt = aem.F.isdIdt
     dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, lowpassfcs=sounding.lowpassfcs,
-        times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade,
+        times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, isdIdt,
         rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx,
         z=copy(aem.z), ρ=copy(aem.ρ), 
         aem.F.calcjacobian, aem.useML, showgates=false)

From e405cfed21179e9989754724cf7569faf97b3a68 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Thu, 29 Aug 2024 20:47:45 +1000
Subject: [PATCH 08/17] explanatory notes on ramp choices

---
 src/AEM_VMD_HMD.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index cf926dc..004970b 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -670,10 +670,14 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
         for iramp = 1:size(F.ramp,1)-1
             rta, rtb  = F.ramp[iramp,1], F.ramp[iramp+1,1]
             if !F.isdIdt
+                # choice made: ramp is inbetween rta and rtb if current is given
                 dt   = rtb - rta
                 dI   = F.ramp[iramp+1,2] - F.ramp[iramp,2]
                 dIdt = dI/dt
             else
+                # choice made: ramp[i+1] is result of current changes between t[i] and t[i+1]
+                # makes causal physical sense
+                # NRG seems to prefer this as ramp[1] at t[1]=0 is zero.
                 dIdt = F.ramp[iramp+1,2]
             end
             if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time

From c48e13ad5ec3ab1addb3d0c5411587ce6f1225f9 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Thu, 29 Aug 2024 21:04:48 +1000
Subject: [PATCH 09/17] explanatory notes on ramp choices

---
 src/AEM_VMD_HMD.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index 004970b..f999eb7 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -671,13 +671,16 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
             rta, rtb  = F.ramp[iramp,1], F.ramp[iramp+1,1]
             if !F.isdIdt
                 # choice made: ramp is inbetween rta and rtb if current is given
+                # easy convention when I(t) is provided
                 dt   = rtb - rta
                 dI   = F.ramp[iramp+1,2] - F.ramp[iramp,2]
                 dIdt = dI/dt
             else
                 # choice made: ramp[i+1] is result of current changes between t[i] and t[i+1]
-                # makes causal physical sense
+                # makes causal physical sense, ramp[1] is not used at t[1] = 0
                 # NRG seems to prefer this as ramp[1] at t[1]=0 is zero.
+                # other convention possible is ramp[i] is due to current changes between t[i] and t[i+1]
+                # this will not use ramp[end]
                 dIdt = F.ramp[iramp+1,2]
             end
             if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time

From 814928b1eb1548ccd0bd7663e25a851c1eca87b8 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Fri, 30 Aug 2024 12:57:18 +1000
Subject: [PATCH 10/17] explicit ramp choices

---
 src/AEM_VMD_HMD.jl     | 25 +++++++++++++++++--------
 src/CommonToAll.jl     |  8 +++++++-
 src/VTEM1DInversion.jl | 14 +++++++++-----
 3 files changed, 33 insertions(+), 14 deletions(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index f999eb7..d37584f 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -71,6 +71,7 @@ mutable struct HFieldDHT <: HField
     HTD_r_J_interp     
     dBrdt_J
     isdIdt
+    rampchoice
 end
 
 function HFieldDHT(;
@@ -96,11 +97,15 @@ function HFieldDHT(;
       freqhigh = 1e6,
       minresptime = 1.e-6, # I think responses earlier than this are unstable
       calcjacobian = false,
-      isdIdt = false
+      isdIdt = false,
+      rampchoice = :next, # if using dIdt instead of I a choice has to be made
   )
     @assert all(freqs .> 0.)
     @assert freqhigh > freqlow
     @assert all(diff(times) .> 0)
+    if isdIdt
+        @assert (rampchoice == :previous) | (rampchoice == :next)
+    end
     thickness = zeros(nmax)
     zintfc    = zeros(nmax)
     pz        = zeros(Complex{Float64}, nmax)
@@ -180,7 +185,7 @@ function HFieldDHT(;
             quadnodes, quadweights, preallocate_ω_Hsc(interptimes, lowpassfcs)..., rxwithinloop, provideddt, doconvramp, useprimary,
             nkᵣeval, interpkᵣ, log10interpkᵣ, log10Filter_base, getradialH, getazimH, 
             calcjacobian, Jtemp, similar(Jtemp), Jac_z, Jac_az, Jac_r, HFD_z_J, HTD_z_J_interp, dBzdt_J,  HFD_az_J, HTD_az_J_interp, dBazdt_J, 
-            HFD_r_J, HTD_r_J_interp, dBrdt_J, isdIdt)
+            HFD_r_J, HTD_r_J_interp, dBrdt_J, isdIdt, rampchoice)
 end
 
 function checkrampformintime(times, ramp, minresptime, maxtime)
@@ -676,12 +681,16 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
                 dI   = F.ramp[iramp+1,2] - F.ramp[iramp,2]
                 dIdt = dI/dt
             else
-                # choice made: ramp[i+1] is result of current changes between t[i] and t[i+1]
-                # makes causal physical sense, ramp[1] is not used at t[1] = 0
-                # NRG seems to prefer this as ramp[1] at t[1]=0 is zero.
-                # other convention possible is ramp[i] is due to current changes between t[i] and t[i+1]
-                # this will not use ramp[end]
-                dIdt = F.ramp[iramp+1,2]
+                if F.rampchoice == :next
+                    # choice made: ramp[i+1] is result of current changes between t[i] and t[i+1]
+                    # makes causal physical sense, ramp[1] is not used at t[1] = 0
+                    # NRG seems to prefer this as ramp[1] at t[1]=0 is zero.
+                    dIdt = F.ramp[iramp+1,2]
+                elseif F.rampchoice == :previous
+                    # other convention possible is ramp[i] is due to current changes between t[i] and t[i+1]
+                    # this will not use ramp[end]
+                    dIdt = F.ramp[iramp,2]
+                end
             end
             if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time
                 break
diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl
index d36f4d6..9cf5b90 100644
--- a/src/CommonToAll.jl
+++ b/src/CommonToAll.jl
@@ -15,7 +15,7 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo
         makegrid, whichislast, makesummarygrid, makearray, plotNEWSlabels, 
         plotprofile, gridpoints, splitsoundingsbyline, getsoundingsperline, docontinue, linestartend,
         compatidxwarn, dfn2hdr, getgdfprefix, readlargetextmatrix, pairinteractionplot, flipline, 
-        summaryconductivity, plotsummarygrids1, getVE, writevtkfromsounding, 
+        summaryconductivity, plotsummarygrids1, getVE, writevtkfromsounding, trapezoidrule,
         readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, 
         writeijkfromsounding, nanmean, infmean, nanstd, infstd, infnanmean, infnanstd, 
         kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, 
@@ -898,6 +898,12 @@ function secondderiv(x)
     abs.(sd)
 end
 
+function trapezoidrule(f,t)
+    dt = diff(t)
+    g  = 0.5*(f[1:end-1]+f[2:end])
+    cumsum(g.*dt), 0.5*(t[1:end-1]+t[2:end])
+end
+
 function dodensityestimate(usekde::Bool, data, K::KDEtype, edges)
     if usekde
         kdefunc = kde_(K, data)
diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl
index 6c7089d..7da9da1 100644
--- a/src/VTEM1DInversion.jl
+++ b/src/VTEM1DInversion.jl
@@ -50,7 +50,8 @@ function dBzdt(;times           = [1.],
         lowpassfcs      = [],
         freqlow         = 1e-4,
         freqhigh        = 1e6,
-        isdIdt          = false
+        isdIdt          = false, 
+        rampchoice      = :next
         )
    
     @assert size(σ)  == size(d)
@@ -62,7 +63,7 @@ function dBzdt(;times           = [1.],
         lowpassfcs, freqlow, freqhigh,
         calcjacobian, nfreqsperdecade,
         ntimesperdecade, modelprimary,
-        isdIdt
+        isdIdt, rampchoice
     )
     @assert length(F.thickness) >= length(z)
     # for Gauss-Newton
@@ -406,7 +407,8 @@ function makeoperator(sounding::VTEMsoundingData;
             showgeomplot  = false,
             calcjacobian  = false,
             plotfield     = false,
-            isdIdt        = false
+            isdIdt        = false,
+            rampchoice    = :next,
             )
     
     zall, znall, zboundaries = setupz(zstart, extendfrac, dz=dz, n=nlayers, showplot=showgeomplot)
@@ -414,7 +416,8 @@ function makeoperator(sounding::VTEMsoundingData;
     ρ[z.>=zstart] .= ρbg
     aem = dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary,
         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, isdIdt)
+        rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, z, ρ, calcjacobian, useML, showgates=plotfield, isdIdt,
+        rampchoice)
     plotfield && plotmodelfield!(aem, log10.(ρ[2:end]))
     aem, zall, znall, zboundaries
 end
@@ -424,8 +427,9 @@ function makeoperator(aem::dBzdt, sounding::VTEMsoundingData)
     nfreqsperdecade = gettimesperdec(aem.F.freqs)
     modelprimary = aem.F.useprimary === 1. ? true : false
     isdIdt = aem.F.isdIdt
+    rampchoice = aem.F.rampchoice
     dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, lowpassfcs=sounding.lowpassfcs,
-        times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, isdIdt,
+        times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, isdIdt, rampchoice,
         rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx,
         z=copy(aem.z), ρ=copy(aem.ρ), 
         aem.F.calcjacobian, aem.useML, showgates=false)

From cd330f501add665c3c8dc6043934b381dd1dee9a Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Tue, 3 Sep 2024 17:59:02 +1000
Subject: [PATCH 11/17] WIP finished for dIdt waveforms: -for dIdT,
 rampchoice=:mid seems most sensible -for I=integrate(dIdt), using t[2:end]
 seems best

---
 src/AEM_VMD_HMD.jl     | 30 ++++++++++++++++++++++++++----
 src/CommonToAll.jl     |  4 +++-
 src/VTEM1DInversion.jl |  4 ++--
 3 files changed, 31 insertions(+), 7 deletions(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index d37584f..c92aafb 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -98,13 +98,13 @@ function HFieldDHT(;
       minresptime = 1.e-6, # I think responses earlier than this are unstable
       calcjacobian = false,
       isdIdt = false,
-      rampchoice = :next, # if using dIdt instead of I a choice has to be made
+      rampchoice = :mid, # if using dIdt instead of I a choice has to be made
   )
     @assert all(freqs .> 0.)
     @assert freqhigh > freqlow
     @assert all(diff(times) .> 0)
     if isdIdt
-        @assert (rampchoice == :previous) | (rampchoice == :next)
+        @assert (rampchoice == :previous) | (rampchoice == :next) | (rampchoice == :mid)
     end
     thickness = zeros(nmax)
     zintfc    = zeros(nmax)
@@ -671,9 +671,16 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
     end    
     F.getradialH && fill!(F.dBrdt, 0.)
     F.getazimH && fill!(F.dBazdt, 0.)
+    nramps = size(F.ramp,1)
+    if F.isdIdt && F.rampchoice == :mid
+        nramps += 1
+    end
     for itime = 1:length(F.times)
-        for iramp = 1:size(F.ramp,1)-1
-            rta, rtb  = F.ramp[iramp,1], F.ramp[iramp+1,1]
+        for iramp = 1:nramps-1
+            if F.rampchoice != :mid
+                rta, rtb  = F.ramp[iramp,1], F.ramp[iramp+1,1]
+            end
+            # end
             if !F.isdIdt
                 # choice made: ramp is inbetween rta and rtb if current is given
                 # easy convention when I(t) is provided
@@ -690,8 +697,23 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
                     # other convention possible is ramp[i] is due to current changes between t[i] and t[i+1]
                     # this will not use ramp[end]
                     dIdt = F.ramp[iramp,2]
+                else # mid
+                    if iramp == 1
+                        rta = F.ramp[1,1]
+                        rtb = (F.ramp[1,1]+F.ramp[2,1])/2
+                        dIdt = F.ramp[1,2]
+                    elseif iramp == nramps-1
+                        rta = (F.ramp[end-1,1]+F.ramp[end,1])/2
+                        rtb = F.ramp[end,1]
+                        dIdt = F.ramp[end,2]
+                    else
+                        rta = (F.ramp[iramp-1,1]+F.ramp[iramp,1])/2
+                        rtb = (F.ramp[iramp,1]+F.ramp[iramp+1,1])/2
+                        dIdt = F.ramp[iramp,2]
+                    end
                 end
             end
+            itime ==1 && @info rta, rtb, dIdt
             if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time
                 break
             end
diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl
index 9cf5b90..671cb76 100644
--- a/src/CommonToAll.jl
+++ b/src/CommonToAll.jl
@@ -901,7 +901,9 @@ end
 function trapezoidrule(f,t)
     dt = diff(t)
     g  = 0.5*(f[1:end-1]+f[2:end])
-    cumsum(g.*dt), 0.5*(t[1:end-1]+t[2:end])
+    # integrated, tstart, tmid, tend
+    # I recommend using tend
+    cumsum(g.*dt), t[1:end-1], 0.5*(t[1:end-1]+t[2:end]), t[2:end]
 end
 
 function dodensityestimate(usekde::Bool, data, K::KDEtype, edges)
diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl
index 7da9da1..54154ac 100644
--- a/src/VTEM1DInversion.jl
+++ b/src/VTEM1DInversion.jl
@@ -51,7 +51,7 @@ function dBzdt(;times           = [1.],
         freqlow         = 1e-4,
         freqhigh        = 1e6,
         isdIdt          = false, 
-        rampchoice      = :next
+        rampchoice      = :mid
         )
    
     @assert size(σ)  == size(d)
@@ -408,7 +408,7 @@ function makeoperator(sounding::VTEMsoundingData;
             calcjacobian  = false,
             plotfield     = false,
             isdIdt        = false,
-            rampchoice    = :next,
+            rampchoice    = :mid,
             )
     
     zall, znall, zboundaries = setupz(zstart, extendfrac, dz=dz, n=nlayers, showplot=showgeomplot)

From cb32d5ecc42bc0e53d1790c6fe6df31eb7e3df29 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Wed, 4 Sep 2024 10:54:00 +1000
Subject: [PATCH 12/17] changed checkrampformintime to align with last WIP

---
 src/AEM_VMD_HMD.jl | 30 +++++++++++++++++++++++-------
 1 file changed, 23 insertions(+), 7 deletions(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index c92aafb..5d9758b 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -99,7 +99,7 @@ function HFieldDHT(;
       calcjacobian = false,
       isdIdt = false,
       rampchoice = :mid, # if using dIdt instead of I a choice has to be made
-  )
+      )
     @assert all(freqs .> 0.)
     @assert freqhigh > freqlow
     @assert all(diff(times) .> 0)
@@ -119,7 +119,7 @@ function HFieldDHT(;
     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)
+         mintime, maxtime = checkrampformintime(times, ramp, minresptime, maxtime, rampchoice, isdIdt)
     end
     interptimes = 10 .^(log10(mintime) : 1/ntimesperdecade : log10(maxtime))
     if freqhigh < 3/mintime
@@ -188,13 +188,30 @@ function HFieldDHT(;
             HFD_r_J, HTD_r_J_interp, dBrdt_J, isdIdt, rampchoice)
 end
 
-function checkrampformintime(times, ramp, minresptime, maxtime)
-    # this checks we don't have ultra small tobs - ramp_time_a
+function checkrampformintime(times, ramp, minresptime, maxtime, rampchoice, isdIdt)
+    # this checks we don't have ultra small t_obs - ramp_time_a
     minta = Inf
     maxta = -Inf
+    nramps = size(ramp, 1)
+    if isdIdt && rampchoice == :mid
+        nramps += 1
+    end
     for itime = 1:length(times)
-        for iramp = 1:size(ramp,1)-1
-            rta, rtb  = ramp[iramp,1], ramp[iramp+1,1]
+        for iramp = 1:nramps-1
+            if (rampchoice == :mid) && isdIdt # receiver voltage waveform (dIdt) and :mid
+                if iramp == 1
+                    rta = F.ramp[1,1]
+                    rtb = (F.ramp[1,1]+F.ramp[2,1])/2
+                elseif iramp == nramps-1
+                    rta = (F.ramp[end-1,1]+F.ramp[end,1])/2
+                    rtb = F.ramp[end,1]
+                else
+                    rta = (F.ramp[iramp-1,1]+F.ramp[iramp,1])/2
+                    rtb = (F.ramp[iramp,1]+F.ramp[iramp+1,1])/2
+                end
+            else # when using current (or voltage but not :mid)
+                rta, rtb  = ramp[iramp,1], ramp[iramp+1,1]
+            end
             if rta >= times[itime] # geq instead of gt as we could have an unlcky time
                 break
             end
@@ -713,7 +730,6 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
                     end
                 end
             end
-            itime ==1 && @info rta, rtb, dIdt
             if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time
                 break
             end

From e9edb5576a413d68400854613bf959eecf439deb Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Wed, 4 Sep 2024 15:40:29 +1000
Subject: [PATCH 13/17] scaled RDP n-dimensional added

---
 zz_portalcurtains/RDP.jl | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/zz_portalcurtains/RDP.jl b/zz_portalcurtains/RDP.jl
index 7e5d1eb..7538309 100644
--- a/zz_portalcurtains/RDP.jl
+++ b/zz_portalcurtains/RDP.jl
@@ -35,6 +35,19 @@ function rdpreduce(pointlist, ϵ=1.0)
     keep
 end
 
+function scaledRDP(xin...;ϵ=1.0) # input x,y,z,etc. 
+    scaled = map(xin) do x
+        (x .- minimum(x))/(maximum(x) - minimum(x))
+    end
+    plistscaled = map(zip(scaled...)) do (x)
+         [x...]
+    end
+    scaledpgood = reduce(vcat, RDP.rdpreduce(plistscaled, ϵ)')
+    scaledout = map(zip(eachcol(scaledpgood), xin)) do (s, x) 
+        s*(maximum(x) - minimum(x)) .+ minimum(x)
+    end
+end
+
 function worldcoordinates(;gridr=nothing, gridz=nothing)
     # from https://support.esri.com/en-us/knowledge-base/faq-what-is-the-format-of-the-world-file-used-for-geore-000002860
     A = step(gridr)

From d6e2939cc7f7d209e62ae513a712ed4786c734e8 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Thu, 5 Sep 2024 20:28:12 +1000
Subject: [PATCH 14/17] had forgotten rtb in \!dIdt

---
 src/AEM_VMD_HMD.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index 5d9758b..67c9ed8 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -694,7 +694,7 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
     end
     for itime = 1:length(F.times)
         for iramp = 1:nramps-1
-            if F.rampchoice != :mid
+            if (F.rampchoice != :mid) || (!F.isdIdt)
                 rta, rtb  = F.ramp[iramp,1], F.ramp[iramp+1,1]
             end
             # end

From 84772ed6ca4c1a25af098248232cb51f2cb52f59 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Fri, 6 Sep 2024 18:12:17 +1000
Subject: [PATCH 15/17] slight RDP changes

---
 zz_portalcurtains/RDP.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/zz_portalcurtains/RDP.jl b/zz_portalcurtains/RDP.jl
index 7538309..a76ac9c 100644
--- a/zz_portalcurtains/RDP.jl
+++ b/zz_portalcurtains/RDP.jl
@@ -266,7 +266,7 @@ function writeimageandcolorbar(img::Array, gridr, gridz, line::Int; cmap="turbo"
     end
     figsize = gridr[end]/shrink, abs(diff([extrema(gridz)...])[1])*VE/shrink
     fig, ax = plt.subplots(;figsize)
-    ax.imshow(img; cmap, extent=[gridr[1], gridr[end], gridz[end], gridz[1]])
+    ax.imshow(img; cmap, extent=[gridr[1], gridr[end], gridz[end], gridz[1]], vmin, vmax)
     ax.set_aspect(VE)
     ax.set_aspect("auto")
     ax.axis("off")

From 78bee0d6343d3f14e80c308a56bb0a90c21121ca Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Mon, 9 Sep 2024 12:09:17 +1000
Subject: [PATCH 16/17] F.ramp -> ramp

---
 src/AEM_VMD_HMD.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index 67c9ed8..7eea2c5 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -200,14 +200,14 @@ function checkrampformintime(times, ramp, minresptime, maxtime, rampchoice, isdI
         for iramp = 1:nramps-1
             if (rampchoice == :mid) && isdIdt # receiver voltage waveform (dIdt) and :mid
                 if iramp == 1
-                    rta = F.ramp[1,1]
-                    rtb = (F.ramp[1,1]+F.ramp[2,1])/2
+                    rta = ramp[1,1]
+                    rtb = (ramp[1,1]+ramp[2,1])/2
                 elseif iramp == nramps-1
-                    rta = (F.ramp[end-1,1]+F.ramp[end,1])/2
-                    rtb = F.ramp[end,1]
+                    rta = (ramp[end-1,1]+ramp[end,1])/2
+                    rtb = ramp[end,1]
                 else
-                    rta = (F.ramp[iramp-1,1]+F.ramp[iramp,1])/2
-                    rtb = (F.ramp[iramp,1]+F.ramp[iramp+1,1])/2
+                    rta = (ramp[iramp-1,1]+ramp[iramp,1])/2
+                    rtb = (ramp[iramp,1]+ramp[iramp+1,1])/2
                 end
             else # when using current (or voltage but not :mid)
                 rta, rtb  = ramp[iramp,1], ramp[iramp+1,1]

From 7fc96ae1bb252c499d8914cc11060e798e71cea7 Mon Sep 17 00:00:00 2001
From: Anandaroop Ray <a2ray@ucsd.edu>
Date: Mon, 9 Sep 2024 17:16:56 +1000
Subject: [PATCH 17/17] don't ever model to less than interptimes[1]

---
 src/AEM_VMD_HMD.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl
index 7eea2c5..37e5b46 100644
--- a/src/AEM_VMD_HMD.jl
+++ b/src/AEM_VMD_HMD.jl
@@ -738,10 +738,13 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
             end
 
             ta = F.times[itime]-rta
-            if ta <= F.minresptime
-                break # since ta must always be greater than minresptime
+            if ta < F.interptimes[1]
+                break # since ta must always be greater than min response time modeled
+            end
+            tb = max(F.times[itime]-rtb, F.interptimes[1]) # rtb > rta, so make sure this is not zero because integ is in log10...
+            if tb >= ta # I don't actually expect to get here
+                break # because int is from ta to tb and ta=t-rta tb=t-rtb and we must have ta>tb  
             end
-            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