some minor bugfixes

src/Acquisition.jl

@@ -14,6 +14,7 @@ import JuMIT.Models
 import JuMIT.Wavelets
 import JuMIT.DSP
 using Distributions
+using DSP # from Julia
 
 """
 Acquisiton has supersources, sources and receivers.
@@ -128,7 +129,7 @@ recorded.
 """
 function Geom_find(acq::Geom; rpos::Array{Float64,1}=nothing, rpos0::Array{Float64,1}=nothing)
 	rpos==nothing ? error("need rpos") : nothing
-	sson = Array(Vector{Int64}, acq.nss);
+	sson = Array{Vector{Int64}}(acq.nss);
 	for iss = 1:acq.nss
 		rvec = [[acq.rz[iss][ir],acq.rx[iss][ir]] for ir=1:acq.nr[iss]]
 		ir=findfirst(rvec, rpos)
@@ -456,42 +457,81 @@ Constructor of `Src`, which is typical for a input model such that
 the model has `nλ` wavelengths.
 
 # Arguments
-* `nλ::Int64=10` : number of wavelenghts in the model
-nλ
-"""
-function Src_fixed_mod(nss::Int64, ns::Int64, nfield::Int64, model::Models.Seismic, nλ::Int64=10)
 
-	x=model.mgrid.x; z=model.mgrid.z
-	# dominant wavelength using model dimensions
+* `nss::Int64` : number of supersources
+* `ns::Int64` : number of source per each supersource
+* `nfield::Int64` :
+
+# Keyword Arguments
+
+* `mod::Models.Seismic` :
+* `nλ::Int64=10` : number of wavelengths in the mod
+* `wav_func::Function=(fqdom, tgrid)->Wavelets.ricker(fqdom,tgrid)` : which wavelet to generate, see Wavelets.jl
+* `tmaxfrac::Float64=1.0` : by default the maximum modelling time is computed using the average velocity and the diagonal distance of the model, use this fraction to increase of reduce the maximum time
+"""
+function Src_fixed_mod(
+		nss::Int64, ns::Int64, nfield::Int64; 
+		mod::Models.Seismic=nothing, 
+		nλ::Int64=10,
+		wav_func::Function=(fqdom, tgrid)->Wavelets.ricker(fqdom,tgrid),
+		tmaxfrac::Float64=1.0
+		)
+
+	x=mod.mgrid.x; z=mod.mgrid.z
+	# dominant wavelength using mod dimensions
 	λdom=mean([(abs(x[end]-x[1])), (abs(z[end]-z[1]))])/real(nλ)
 	# average P velocity
-	vavg=Models.χ([mean(model.χvp)], model.vp0, -1)[1]
+	vavg=Models.χ([mean(mod.χvp)], mod.vp0, -1)[1]
 
 	fqdom = vavg/λdom
 
 	# maximum distance the wave travels
 	d = sqrt((x[1]-x[end]).^2+(z[1]-z[end]).^2)
 
 	# use two-way maximum distance to get tmax
-	tmax=2.*d/vavg
+	tmax=2.*d/vavg .* tmaxfrac
 
 	# choose sampling interval to obey max freq of source wavelet
-	δmin = minimum([model.mgrid.δx, model.mgrid.δz])
-	vpmax = model.vp0[2]
+	δmin = minimum([mod.mgrid.δx, mod.mgrid.δz])
+	vpmax = mod.vp0[2]
 	δt=0.5*δmin/vpmax
 
 	# check if δt is reasonable
 	#(δt > 0.1/fqdom) : error("decrease spatial sampling or nλ")
 
 	tgrid=Grid.M1D(0.0, tmax, δt)
-	wav = Wavelets.ricker(fqdom=fqdom, tgrid=tgrid);
+	wav = wav_func(fqdom, tgrid);
 
 	src=Src_fixed(nss, ns, nfield, wav, tgrid)
 	print(src)
 	# choose ricker waveletes of fdom
 	return src
 end
 
+"""
+Generate band-limited random source signals 
+"""
+function Src_fixed_random(nss::Int64, ns::Int64, nfield::Int64, fmin::Float64, fmax::Float64, tgrid::Grid.M1D, tmax::Float64=tgrid.x[end] )
+	fs = 1/ tgrid.δx;
+	designmethod = Butterworth(4);
+	filtsource = Bandpass(fmin, fmax; fs=fs);
+
+	itmax = indmin(abs.(tgrid.x-tmax))
+	# 20% taper window
+	twin = DSP.taper(ones(itmax),20.) 
+	X = zeros(itmax)
+	wavsrc = [repeat(zeros(tgrid.nx),inner=(1,ns)) for iss=1:nss, ifield=1:nfield] 
+	for ifield in 1:nfield, iss in 1:nss, is in 1:ns
+		X[:] = rand(Uniform(-1.0, 1.0), itmax) .* twin
+		# band limit
+		filt!(X, digitalfilter(filtsource, designmethod), X);
+		wavsrc[iss, ifield][1:itmax,is] = X.*twin
+	end
+	src= Src(nss, fill(ns, nss), nfield, wavsrc, deepcopy(tgrid))
+	print(src)
+	return src
+end
+
 """
 Function that returns Src after time reversal
 """

src/DSP.jl

@@ -7,7 +7,11 @@ using Distributions
 using DSP # from julia
 
 """
-generate random signals in the frequency domain 
+Generate a band-limited 
+random signal of a particular maximum time, `tmax`.
+
+
+
 such that a box car function is applied in the time domain
 """
 function get_random_tmax_signal(;
@@ -18,44 +22,44 @@ function get_random_tmax_signal(;
 	dist::Symbol=:gaussian # distribution type
 	)
 
-# initialize outputs
-S = fill(complex(0.0,0.0),fgrid.nx);
-Sreal = fill(0.0,fgrid.nx);
-s = fill(complex(0.0,0.0),fgrid.nx);
-
-Δf = 1.0 / tmax
-Δf <= fgrid.δx ? error("sampling smaller than grid sampling") :
-Δf >= (fmax-fmin) ? error("need to increase tmax") :
-fvec = [f for f in fmin:Δf:fmax]
-ifvec = fill(0, size(fvec))
-
-println("number of frequencies added to signal:\t", size(fvec,1))
-println("interval between random variable:\t", Δf)
-println("minimum frequency added to signal\t",minimum(fvec))
-println("maximum frequency added to signal\t",maximum(fvec))
-for iff in eachindex(fvec)
-	ifvec[iff] =  indmin((fgrid.x - fvec[iff]).^2.0)
-end
+	# initialize outputs
+	S = fill(complex(0.0,0.0),fgrid.nx);
+	Sreal = fill(0.0,fgrid.nx);
+	s = fill(complex(0.0,0.0),fgrid.nx);
+
+	Δf = 1.0 / tmax
+	Δf <= fgrid.δx ? error("sampling smaller than grid sampling") :
+	Δf >= (fmax-fmin) ? error("need to increase tmax") :
+	fvec = [f for f in fmin:Δf:fmax]
+	ifvec = fill(0, size(fvec))
+
+	println("number of frequencies added to signal:\t", size(fvec,1))
+	println("interval between random variable:\t", Δf)
+	println("minimum frequency added to signal\t",minimum(fvec))
+	println("maximum frequency added to signal\t",maximum(fvec))
+	for iff in eachindex(fvec)
+		ifvec[iff] =  indmin((fgrid.x - fvec[iff]).^2.0)
+	end
 
-if(dist == :gaussian)
-	X = randn(size(fvec));
-elseif(dist == :uniform)
-	X = rand(Uniform(-2.0, 2.0), size(fvec))
-else
-	error("invalid dist")
-end
+	if(dist == :gaussian)
+		X = randn(size(fvec));
+	elseif(dist == :uniform)
+		X = rand(Uniform(-2.0, 2.0), size(fvec))
+	else
+		error("invalid dist")
+	end
 
-for iff in eachindex(fvec)
-	Sreal += X[iff] .* sinc(tmax.*(abs(fgrid.x - fgrid.x[ifvec[iff]])))
-end
+	for iff in eachindex(fvec)
+		Sreal += X[iff] .* sinc(tmax.*(abs(fgrid.x - fgrid.x[ifvec[iff]])))
+	end
 
-S = complex.(Sreal, 0.0);
+	S = complex.(Sreal, 0.0);
 
-# remove mean 
-#S[1] = 0.0; 
-s = ifft(S); 
+	# remove mean 
+	#S[1] = 0.0; 
+	s = ifft(S); 
 
-return S, s
+	return S, s
 
 end
 

src/Gallery.jl

@@ -75,7 +75,19 @@ Gallery of `Seismic` models.
   * `attrib=:seismic_marmousi2_box2` : box for surface seismic studies
 """
 function Seismic(attrib::Symbol, δ::Float64=0.0)
-	if((attrib == :acou_homo1) | (attrib == :acou_homo2))
+	bfrac=0.1; 
+	if((attrib == :acou_homo1))
+		vp0 = [1500., 3500.] # bounds for vp
+		vs0 = [1.0, 1.0] # dummy
+		ρ0 = [1500., 3500.] # density bounds
+		mgrid = M2D(attrib)
+		model= Models.Seismic(vp0, vs0, ρ0,
+		      fill(0.0, (mgrid.nz, mgrid.nx)),
+		      fill(0.0, (mgrid.nz, mgrid.nx)),
+		      fill(0.0, (mgrid.nz, mgrid.nx)),
+		      mgrid)
+
+	elseif((attrib == :acou_homo2))
 		vp0 = [1700., 2300.] # bounds for vp
 		vs0 = [1.0, 1.0] # dummy
 		ρ0 = [1700., 2300.] # density bounds
@@ -91,44 +103,34 @@ function Seismic(attrib::Symbol, δ::Float64=0.0)
 		vs, h= IO.readsu(joinpath(marmousi_folder,"vs_marmousi-ii_0.1.su"))
 		ρ,  h= IO.readsu(joinpath(marmousi_folder,"density_marmousi-ii_0.1.su"))
 		vp .*= 1000.; vs .*= 1000.; #ρ .*=1000
-		bound=0.1; vp0=zeros(2); vs0=zeros(2); ρ0=zeros(2);
-		boundvp=bound*mean(vp); boundvs=bound*mean(vs); boundρ=bound*mean(ρ);
-		vp0[1] = ((minimum(vp) - boundvp)<0.0) ? 0.0 : (minimum(vp) - boundvp)
-		vp0[2] = maximum(vp)+boundvp
-		vs0[1] = ((minimum(vs) - boundvs)<0.0) ? 0.0 : (minimum(vs) - boundvs)
-		vs0[2] = maximum(vs)+boundvs
-		ρ0[1] = ((minimum(ρ) - boundρ)<0.0) ? 0.0 : (minimum(ρ) - boundρ)
-		ρ0[2] = maximum(ρ)+boundρ
+		vp0=Models.bounds(vp,bfrac); 
+		vs0=Models.bounds(vs,bfrac); 
+		ρ0=Models.bounds(ρ, bfrac);
 		mgrid = Grid.M2D(0., 17000., 0., 3500.,size(vp,2),size(vp,1),40)
 		model= Models.Seismic(vp0, vs0, ρ0, Models.χ(vp,vp0,1), Models.χ(vs,vs0,1), Models.χ(ρ,ρ0,1), mgrid)
 	elseif(attrib == :seismic_marmousi2_high_res)
 		vp, h= IO.readsu(joinpath(marmousi_folder,"vp_marmousi-ii.su"))
 		vs, h= IO.readsu(joinpath(marmousi_folder,"vs_marmousi-ii.su"))
 		ρ,  h= IO.readsu(joinpath(marmousi_folder,"density_marmousi-ii.su"))
 		vp .*= 1000.; vs .*= 1000.; #ρ .*=1000
-		bound=0.1; vp0=zeros(2); vs0=zeros(2); ρ0=zeros(2);
-		boundvp=bound*mean(vp); boundvs=bound*mean(vs); boundρ=bound*mean(ρ);
-		vp0[1] = ((minimum(vp) - boundvp)<0.0) ? 0.0 : (minimum(vp) - boundvp)
-		vp0[2] = maximum(vp)+boundvp
-		vs0[1] = ((minimum(vs) - boundvs)<0.0) ? 0.0 : (minimum(vs) - boundvs)
-		vs0[2] = maximum(vs)+boundvs
-		ρ0[1] = ((minimum(ρ) - boundρ)<0.0) ? 0.0 : (minimum(ρ) - boundρ)
-		ρ0[2] = maximum(ρ)+boundρ
+		vp0=Models.bounds(vp,bfrac); 
+		vs0=Models.bounds(vs,bfrac); 
+		ρ0=Models.bounds(ρ, bfrac);
 		mgrid = Grid.M2D(0., 17000., 0., 3500.,size(vp,2),size(vp,1),40)
 		model= Models.Seismic(vp0, vs0, ρ0, Models.χ(vp,vp0,1), Models.χ(vs,vs0,1), Models.χ(ρ,ρ0,1), mgrid)
 
 	elseif(attrib == :seismic_marmousi2_box1)
 		mgrid=Grid.M2D(8500.,9500., 1000., 2000.,5.,5.,40)
 		marm_box1=Models.Seismic_zeros(mgrid)
 		Models.Seismic_interp_spray!(Seismic(:seismic_marmousi2), marm_box1, :interp, :B1)
-		model= marm_box1
+		model= Models.adjust_bounds!(marm_box1, bfrac)
 	elseif(attrib == :seismic_marmousi2_box2)
 		marm_full = Seismic(:seismic_marmousi2)
-		mgrid=Grid.M2D(5000.,13000., marm_full.mgrid.z[1],
+		mgrid=Grid.M2D(6000.,12000., marm_full.mgrid.z[1],
 		 	marm_full.mgrid.z[end] ,5.,5.,40)
 		marm_box2=Models.Seismic_zeros(mgrid)
 		Models.Seismic_interp_spray!(marm_full, marm_box2, :interp, :B1)
-		model= marm_box2
+		model= Models.adjust_bounds!(marm_box2, bfrac)
 	else
 		error("invalid attrib")
 	end
@@ -156,7 +158,7 @@ The sources and receivers are not placed anywhere on the edges of the mesh.
 * `mgrid::Grid.M2D` : a 2-D mesh
 * `attrib::Symbol` : attribute decides output
   * `=:xwell` cross-well acquisition
-  * `=:surf` cross-well; but sources and receivers at unequal depths
+  * `=:surf` surface acquisition
   * `=:vsp` vertical seismic profiling
   * `=:rvsp`  reverse vertical seismic profiling
   * `=:downhole` downhole sources and receivers 
@@ -165,18 +167,20 @@ The sources and receivers are not placed anywhere on the edges of the mesh.
 """
 function Geom(mgrid::Grid.M2D, attrib::Symbol; nss=2, nr=2, rand_flags=[false, false])
 	otx=(0.9*mgrid.x[1]+0.1*mgrid.x[end]); ntx=(0.1*mgrid.x[1]+0.9*mgrid.x[end]);
+	otwx=(0.95*mgrid.x[1]+0.05*mgrid.x[end]); ntwx=(0.05*mgrid.x[1]+0.95*mgrid.x[end]);
 	otz=(0.9*mgrid.z[1]+0.1*mgrid.z[end]); ntz=(0.1*mgrid.z[1]+0.9*mgrid.z[end]);
+	otwz=(0.95*mgrid.z[1]+0.05*mgrid.z[end]); ntwz=(0.05*mgrid.z[1]+0.95*mgrid.z[end]);
 	quatx = (0.75*mgrid.x[1]+0.25*mgrid.x[end]); quatz = (0.75*mgrid.z[1]+0.25*mgrid.z[end]) 
 	tquatx = (0.25*mgrid.x[1]+0.75*mgrid.x[end]); tquatz = (0.25*mgrid.z[1]+0.75*mgrid.z[end]) 
 	halfx = 0.5*(mgrid.x[1]+mgrid.x[end]);	halfz = 0.5*(mgrid.z[1]+mgrid.z[end]);
 	if(attrib == :xwell)
 		geom=Acquisition.Geom_fixed(otz, ntz, otx, otz, ntz, ntx, nss, nr, :vertical, :vertical, rand_flags)
 	elseif(attrib == :surf)
-		geom=Acquisition.Geom_fixed(otx, ntx, otz, otx, ntx, otz, nss, nr, :horizontal, :horizontal, rand_flags)
+		geom=Acquisition.Geom_fixed(otx, ntx, otwz, otx, ntx, otwz, nss, nr, :horizontal, :horizontal, rand_flags)
 	elseif(attrib == :vsp)
-		geom=Acquisition.Geom_fixed(otx, ntx, otz, otz, ntz, otx, nss, nr, :horizontal, :vertical, rand_flags)
+		geom=Acquisition.Geom_fixed(otx, ntx, otwz, otz, ntz, otx, nss, nr, :horizontal, :vertical, rand_flags)
 	elseif(attrib == :rvsp)
-		geom=Acquisition.Geom_fixed(otz, ntz, otx, otx, ntx, otz, nss, nr, :vertical, :horizontal, rand_flags)
+		geom=Acquisition.Geom_fixed(otz, ntz, otx, otx, ntx, otwz, nss, nr, :vertical, :horizontal, rand_flags)
 	elseif(attrib == :downhole)
 		geom=Acquisition.Geom_fixed(quatz, otz, quatx, quatz, otz, quatx, nss, nr, :vertical, :vertical, rand_flags)
 	else
@@ -199,15 +203,15 @@ Gallery of source signals `Src`.
 function Src(attrib::Symbol, nss::Int64=1)
 	if(attrib == :acou_homo1)
 		tgrid = M1D(attrib)
-		wav = Wavelets.ricker(fqdom=10.0, tgrid=tgrid, tpeak=0.25, )
+		wav = Wavelets.ricker(10.0, tgrid, tpeak=0.25, )
 		return Acquisition.Src_fixed(nss, 1, 1, wav, tgrid)
 	elseif(attrib == :acou_homo2)
 		tgrid = M1D(attrib)
-		wav = Wavelets.ricker(fqdom=3.0, tgrid=tgrid, tpeak=0.3, )
+		wav = Wavelets.ricker(3.0, tgrid, tpeak=0.3, )
 		return Acquisition.Src_fixed(nss, 1, 1, wav, tgrid)
 	elseif(attrib == :vecacou_homo1)
 		tgrid = M1D(:acou_homo1)
-		wav = Wavelets.ricker(fqdom=10.0, tgrid=tgrid, tpeak=0.25, )
+		wav = Wavelets.ricker(10.0, tgrid, tpeak=0.25, )
 		return Acquisition.Src_fixed(nss, 1, 3, wav, tgrid)
 	end
 end

src/Interferometry.jl

@@ -35,8 +35,8 @@ function TD_virtual_diff(
 	# central wavelength (use zero for testing)
 	println(string("dominant wavelength in the data:\t",λdom))
 
-	rx = Array(Vector{Float64}, nur); rz = Array(Vector{Float64}, nur);
-	datmat = zeros(2*nt-1, nur, nur, data.nfield);
+	rx = Array{Vector{Float64}}(nur); rz = Array{Vector{Float64}}(nur);
+	datmat = zeros(nur, data.nfield, 2*nt-1, nur);
 	for ifield =1:data.nfield
 		# loop over virtual sources
 		for irs = 1:nur
@@ -61,14 +61,14 @@ function TD_virtual_diff(
 					# stacking over these sources
 					for isson=1:length(sson)
 						if(sson[isson] != [0])
-							datmat[:, ir, irs, ifield] += 
+							datmat[ir, ifield,:,irs] += 
 							xcorr(
 							  datan.d[isson,ifield][:, sson[isson][2]], 
 							  datan.d[isson,ifield][:, sson[isson][1]])
 						end
 					end
 					# normalize depending on the stack
-					nsson != 0 ? datmat[:, ir, irs, ifield] /= nsson : nothing
+					nsson != 0 ? datmat[ ir,  ifield, :, irs] /= nsson : nothing
 				end
 			end
 			if(irvec != [])
@@ -82,7 +82,7 @@ function TD_virtual_diff(
 	sz = [[urpos[1][irs]] for irs in 1:nur];
 
 	# bool array acoording to undef
-	ars = [isdefined(rx,irs) ? true : false for irs=1:nur];
+	ars = [isassigned(rx,irs) ? true : false for irs=1:nur];
 
 	# select only positions that are active
 	rx = rx[ars];	rz = rz[ars];