program Passivize
  !----------------------------------------------------------------!
  ! convolves transmission panels w/ white noise (0-1.) trace
  !       requires time domain input
  !       outputs x,1,w for use in migration program
  ! param list:
  !             time: time output,no has only pos. freqs
  !             impulse: source function =1. for testing
  !             allout: dont chop off neg. freq.s 
  !             all: make source last entire length of trace
  !             seed: random # generator seed, and source location
  !             pad: pad length for input
  !             noise: writes file of all noise traces used
  !             shift: mandate constant time shift [-10000 default for no constant shift]
  !             maxlen: max length of source function
  !             decimate: subsample factor for freq. axis
  !             


  use sep
  use fft
  implicit none
  complex,allocatable,dimension(:,:)::data,csource
  complex,allocatable,dimension(:)::trace,out
  real,allocatable,dimension(:)::r2,r1,phase,rout
  real,allocatable,dimension(:,:)::in,source
  integer,allocatable,dimension(:)::z

  integer                                 :: n1,n2,n3,n4,esize,t,n1in,n1out,pad
  real                                    :: d1,d2,d3,d4,pi,w,w0,dw,tottime
  real                                    :: o1,o2,o3,o4,shift,maxlen
  integer                                 :: i1,i2,i3,decimate,b,a,seed,c,d,sn1,sn2,sloc,tb,te
  logical ::time,impulse,all,allout,ss

  call sep_init ()
  call from_param("impulse",impulse,.false.);  call from_param('shift',shift,-10000.)
  call from_param("all",all,.false.);    call from_param("allout",allout,.false.);
  call from_param("decimate",decimate,1);call from_param('seed',seed,1500)
  call from_param("time",time,.false.);  call from_param('maxlen',maxlen,3.)
  call from_param('pad',pad,n1in) 

  call from_history ("n1",n1in); call from_history ("d1",d1); call from_history("o1",o1)
  call from_history ("n2",n2); call from_history ("d2",d2); call from_history("o2",o2)
  call from_history ("n3",n3); call from_history ("d3",d3); call from_history("o3",o3)

  if (pad.gt.n1in) then
     n1=pad
  else
     n1=n1in
  end if
  call from_history("esize",esize,8); if (esize==8)call seperr("real data required")
  if(time)then
     esize=4 ; n1out=n1 ; call to_history("d1",d1);call to_history("label1","s")
  else
     esize=8 
     n1out=n1/decimate;d1=d1*decimate
     if (.not.allout)  n1out=n1out/2

     call to_history("d1",1/(n1*d1));call to_history("label1","1/s")
     if(o1.ne.0.) call seperr("non-zero time origin not acceptable!")
  end if
  call to_history("n1",n1out);call to_history("n2",n2);call to_history("n3",n3)
  call to_history("d2",d2);call to_history("d3",d3)
  call to_history("o1",o1);call to_history("o2",o2);call to_history("o3",o3)
  call to_history("esize",esize);call to_history("label2","m");call to_history("label3"," ")

  call to_history("n4", 1);call to_history("o4",1.);call to_history("d4",1.)
  !----------------------------------------------------------------
!  ss=.true.
!  sloc=1.5+ (seed-o2)/d2 - 4 ! wavelet is ~8 traces thick
!  if (ss) then 
!     if(.not.time) call seperr("no dice")
!     ss=.true.
!     call from_aux("wavelet","n1",sn1)
!     call from_aux("wavelet","n2",sn2)
!     call to_history("n1",n1out,"source")
!     call to_history("n2",n2,"source")
!     allocate(source(sn1,sn2));source=0.
!     call sep_read(source,"wavelet")
!     allocate(csource(n1out,n2));csource=(0.0,0.0)
!     csource=cmplx(source,0.0);
!     do i2=1,sn2
!        call ftu(.false.,csource(:,i2))
!     end do
!    
!  end if

  if (exist_file("noise") )then 
     call to_history("esize",esize,"noise");
     call to_history("n1",n1out,"noise");
     call to_history("n2",n2,"noise");
     call to_history("n3",n3,"noise");
     call to_history("n4",1,"noise");
     call to_history("n5",1,"noise");
  end if
  call sep_close()

  !----------------------------------------------------------------
  allocate(in  (n1in       ,n2));in=0.

  allocate(phase(n1out));phase=cmplx(0.,0.)
  allocate(trace(n1out));trace=cmplx(0.,0.)
  allocate(r2   (n1out));r2=1.

  allocate(data(n1,n2))
  allocate(out (n1out));out=cmplx(0.,0.)
  if(time)allocate(rout(n1out))
  allocate(r1(1000));r1=0. !rand.vec. for source len & bulk time shift calcs
  !----------------------------------------------------------------
  pi=acos(-1.)
  w0=-pi/d1
  dw=2*pi/(n1*d1)
  !----------------------------------------------------------------
  shots:do i3=1,n3
     !-------------------------------------------------------------
     !!calculate the length of the random trace to use
     seed=(i3-1)*25+seed  !seed is panel number, which increments by 25
     a=(seed-o2)/25+1
     call random_number(r1)
     !I think 13 samples is the real wavelet length
     if(.not.all) then
        b=max(13 ,min(size(r1),nint(r1(a)*maxlen/d1)))
     else
        b=n1out-1
     end if
     write(0,*) 'Source length:',(b-1)*d1,'seconds.(',b,'samples)'

     call random_seed(size=a)
     allocate(z(a))
     do i1=1,a
        z(i1)=i1*seed*1000
     end do
     call random_seed(put=z)
     call random_number(r2)
     r2(b+1:n1out)=0. ! source function should only last a few seconds
     if (impulse) then
        trace(1)=cmplx(1.,0.)
        write(0,*) 'Test phase shift, no randomness'
     else
        trace=cmplx(r2,0.)
     end if

     !----------------------------------------------------------------
     !!calculate a random bulk time shift for each panel
     !total amount of time shift allowed is length of time record
     tottime=(n1-1)*d1
     if(shift.eq.-10000.) shift=r2(a)*(tottime-maxlen); 
     write(0,*) 'Time shift=',shift,'of',tottime
     tb=1.5 + (shift-o1)/d1 -1
     te=tb+b-1
     do i1=1,n1out
        w=w0+(i1-1)*dw
        phase(i1)=w*shift
     end do
     !trace(2::2)=-trace(2::2)
     call ftu(.false.,trace)
     trace=trace*cmplx(cos(phase),sin(phase))
!     if (ss) then
!        do i2=1,sn2
!           csource(:,i2)=csource(:,i2)!*cmplx(cos(phase),sin(phase))
!        end do
!     end if

     write(0,*) 'sloc',sloc,size(csource,1)
     !     panel:do t=1,4! attack volume as 4 panels in memory
     call sep_read(in)
     data(1:n1in   ,:)=cmplx(in(:,:),0.)
     data(n1in+1:n1,:)=cmplx(0.     ,0.)
     !data(2::2,:)=-data(2::2,:)
!     a=0
!     b=0
     trc:do i2=1,n2
        call ftu(.false.,data(:,i2))
        out=data(1:n1out,i2)*trace

        if(time) then
           call ftu(.true.,out) 
           !out(2::2)=-out(2::2)
           rout=real(out)
           call srite("out",rout,esize*n1out)
        else
           call srite("out",out,esize*n1out)
        end if

!        if(ss) then
!           if ((i2.lt.sloc).or.(i2.ge.sloc+sn2)) then
!              rout=0.
!              call srite("source",rout,esize*n1out)
!              a=a+1
!           elseif (i2.le.sloc+sn2-1) then
!              b=b+1
!              write(0,*) 'seeded trace',i2
!              call ftu(.true.,out)
!              rout(tb:te)=real(out(tb:te))
!              call srite("source",rout,esize*n1out )
!           end if
!        end if
     end do trc
!write(0,*) 'a,b',a,b
     !     end do panel
     deallocate(z)
     if (exist_file("noise"))then 
        if(time) then 
           call ftu(.true.,trace)
           rout=real(trace)
           call srite("noise",rout,esize*n1out)
        else
           call srite("noise",trace,esize*n1out)
        end if
     end if
  end do shots
  if (time) deallocate(rout)
  deallocate(data,out,r2,in,trace,phase)

end program Passivize
