new version of Plots.seismic

REQUIRE

@@ -1,5 +1,6 @@
 julia 0.6
 DataFrames
+ImageFiltering
 JLD
 CSV
 Distances

src/JuMIT.jl

@@ -23,6 +23,7 @@ include("Analytic.jl")
 include("Fdtd.jl")
 include("FWI.jl")
 include("Interferometry.jl")
+include("SAC.jl")
 include("Plots.jl")
 
 end # module

src/Plots.jl

@@ -6,13 +6,34 @@ using RecipesBase
 using Interpolation
 using Grid
 using Conv
+using ImageFiltering
 import JuMIT.Acquisition
 import JuMIT.Data
 import JuMIT.Models
 
 @userplot Seismic
 
-@recipe function f(p::Seismic; fields=[:vp, :ρ], xlim=nothing, zlim=nothing, geom=nothing) 
+
+"""
+Plot the velocity and density seismic models.
+
+# Arguments
+
+* `model::Models.Seismic` : model that should be plotted
+
+# Keyword Arguments
+
+* `xlim::Vector{Float64}=[model.mgrid.x[1],model.mgrid.x[end]]` : minimum and maximum limits of the second dimension while plotting
+* `zlim::Vector{Float64}=[model.mgrid.z[1],model.mgrid.z[end]]` : minimum and maximum limits of the first dimension while plotting
+* `fields::Vector{Symbol}=[:vp, :ρ]` : fields that are to be plotted, see Models.Seismic_get
+* `contrast_flag=false` : plot only the edges of the model
+
+"""
+
+@recipe function f(p::Seismic; fields=[:vp, :ρ], 
+		   xlim=nothing, zlim=nothing, geom=nothing,
+		   contrast_flag=false
+		   ) 
 
 	model=p.args[1]
 	(xlim===nothing) && (xlim=[model.mgrid.x[1],model.mgrid.x[end]])
@@ -25,55 +46,45 @@ import JuMIT.Models
 	nrow = (model.mgrid.nx <= model.mgrid.nz) ? 1 : length(fields)
 	ncol = (model.mgrid.nx <= model.mgrid.nz) ? length(fields) : 1
 
-	ar=div(model.mgrid.nx,model.mgrid.nz)
-
 	layout := (nrow,ncol)
 	for (i,iff) in enumerate(fields)
+		name=replace(string(iff), "ρ", "rho")
+		name=replace(name, "χ", "contrast\t")
 
 		f0 = Symbol((replace("$(fields[i])", "χ", "")),0)
 		m = Models.Seismic_get(model, fields[i])[izmin:izmax,ixmin:ixmax]
 		mx = model.mgrid.x[ixmin:ixmax]
 		mz = model.mgrid.z[izmin:izmax]
+		if(contrast_flag)
+			mmin=minimum(m)
+			mmax=maximum(m)
+			m=imfilter(m, Kernel.Laplacian())
+			mmmin=minimum(m)
+			mmmax=maximum(m)
+			for j in eachindex(m)
+				m[j]=mmin+((m[j]-mmmin)*inv(mmmax-mmmin))*(mmax-mmin)
+			end
+
+		end
 		@series begin        
 			subplot := i
-			aspect_ratio := ar
+			aspect_ratio := :equal
 			seriestype := :heatmap
-			legend := true
-			xlabel := "x"
-			ylabel := "z"
-			color := :grays
-			yflip := true
-	#		l := :plot
-			title := string(iff)
+			legend --> true
+			xlabel --> "x"
+			ylabel --> "z"
+			color --> :grays
+			title := name
 			xlim := (xlim...)
 			ylim := (zlim...)
-			w := 1
+			#yflip := true
 			mx, mz, m
 		end
 
-		if(!(geom===nothing))
-		end
 	end
 
 end
 
-
-"""
-Plot the velocity and density seismic models.
-
-# Arguments
-
-* `model::Models.Seismic` : model that should be plotted
-
-# Keyword Arguments
-
-* `xlim::Vector{Float64}=[model.mgrid.x[1],model.mgrid.x[end]]` : minimum and maximum limits of the second dimension while plotting
-* `zlim::Vector{Float64}=[model.mgrid.z[1],model.mgrid.z[end]]` : minimum and maximum limits of the first dimension while plotting
-* `fields::Vector{Symbol}=[:vp, :ρ]` : fields that are to be plotted, see Models.Seismic_get
-* `overlay_model=nothing` : use a overlay model
-* `use_bounds=false` : impose bounds from the model or not?
-
-"""
 function Seismic(model::Models.Seismic;
 		 xlim::Vector{Float64}=[model.mgrid.x[1],model.mgrid.x[end]],
 		 zlim::Vector{Float64}=[model.mgrid.z[1],model.mgrid.z[end]],

src/SAC.jl

@@ -0,0 +1,61 @@
+
+module SAC
+
+
+using SAC
+using Interpolation
+using Grid
+using ProgressMeter
+
+
+
+function process(T)
+
+	npts=maximum(getfield.(T,:npts))
+
+	Tpro=SAC.SACtr[]
+
+	@showprogress 1 "Processing..." for ir in eachindex(T)
+		TT=deepcopy(T[ir])
+		if(TT.npts == npts)
+			# normalize
+			SAC.multiply!(TT, inv(vecnorm(TT.t)))
+
+			#SAC.update_headers!(TT)
+			push!(Tpro, TT)
+		end
+	end
+	# sort according to distance
+	dist=getfield.(Tpro, :gcarc)
+	Tpro=Tpro[sortperm(dist)]
+	return Tpro
+end
+
+
+function stackamp(T)
+
+
+	tgrid=Grid.M1D(Float64(T[1].b), Float64(T[1].e), Float64(T[1].delta))
+	fgrid=Grid.M1D_rfft(tgrid)
+
+	ampstack=zeros(fgrid.nx)
+	npts=maximum(getfield.(T,:npts))
+
+	@showprogress 1 "Processing..." for ir in eachindex(T)
+		TT=T[ir]
+		if(TT.npts == npts)
+			amp=abs.(rfft(normalize(Float64.(TT.t))))
+
+			for i in eachindex(amp)
+				ampstack[i]+=amp[i]
+			end
+		end
+	end
+
+	normalize!(ampstack)
+	powwav = (ampstack.^2)
+	powwavdb = 10. * log10.(powwav./maximum(powwav)) # power in decibel after normalizing
+	return powwavdb, fgrid
+end
+
+end