!!
!! The routine mutes data based on pairs t-x (or angle) described in 
!! an input file.
!!
  module mute_mod
!!***************
  implicit none
  integer,private 				:: n1,n2,n3 
  real,private					:: o1,o2,d1,d2 
  real,private,dimension(:,:),allocatable	:: mask2
  real,private,dimension(:),allocatable		:: pad1
  contains

  subroutine mute_init(n1_in,o1_in,d1_in,n2_in,o2_in,d2_in,n3_in,mask2_in)
  integer 				:: n1_in,n2_in,n3_in 
  real					:: o1_in,o2_in,d1_in,d2_in
  real,dimension(n1_in,n2_in)		:: mask2_in 
  n1=n1_in
  n2=n2_in
  n3=n3_in
  o1=o1_in
  o2=o2_in
  d1=d1_in
  d2=d2_in
  allocate(mask2(n1,n2))
  mask2=mask2_in
  end subroutine

  subroutine mute(npe,npi,emx,emt,imx,imt,mask1)
  real					:: time, rnc
  integer, intent(in)			:: npe,npi
  integer				:: if2,if3,n,i2
  integer				:: i1,i3
  integer 				:: mube,mufe,mubi,mufi
  real,dimension(:), intent(in)		:: emt,emx,imt,imx
  real,dimension(:),allocatable		:: pad1
  real,dimension(:), intent(out)	:: mask1
  integer				:: iloc(2),il
allocate(pad1(n1*n2*n3))
  !!
  !! External mute
  !!
do i1=1,n1
  time=o1+(i1-1)*d1
  if (npe.ne.0) then
	  if (time <= emt(1)) then
		if2=int(emx(1)/d2)
	  elseif(time >= emt(npe)) then 
		if2=int(emx(npe)/d2)
	  else
		iloc=minloc(emt, dim=1,MASK=emt>time);il=iloc(1)
		rnc=(emx(il)-emx(il-1))/(emt(il)-emt(il-1))
		if2=int((rnc*(time-emt(il))+emx(il))/d2)
	  end if
  else
	  if2=int(n2/2)
  end if

  !!
  !! Internal mute
  !!
  if (npi.ne.0) then
	  if (time >= imt(npi)) then
		if3=int(imx(npi)/d2)
	  elseif((time >= imt(1)).and.(time < imt(npi))) then 
		iloc=minloc(imt, dim=1,MASK=imt<time);il=iloc(1)
		rnc=(imx(il)-imx(il-1))/(imt(il)-imt(il-1))
		if3=int((rnc*(time-imt(il))+imx(il))/d2)
	  end if
  end if
  mube=int(n2/2)-if2+1; mufe=int(n2/2)+if2
  mubi=int(n2/2)-if3+1; mufi=int(n2/2)+if3+1
  do i2=1,n2
	if ((i2<mube).or.(i2>mufe)) then
		mask2(i1,i2)=0.0
	elseif ((i2>=mubi).and.(i2<=mufi).and.(time>imt(1))) then 
		mask2(i1,i2)=0.0
	elseif ((i2>=mubi).and.(i2<=mufi).and.(time<=imt(1))) then 
		mask2(i1,i2)=1.
	else
		mask2(i1,i2)=1.
	end if
  end do
end do
pad1=0.
mask1=reshape(mask2,(/n1*n2*n3/),pad=pad1)
do i3=1,n3
	mask1((i3-1)*n1*n2+1:i3*n1*n2)=mask1(1:n1*n2)
end do
end subroutine

subroutine mute_close()
    deallocate(pad1)
end subroutine

end module
