Program NewCorr
  !##################################################
  !#                                              
  !#  cross-correlate every trace with every other
  !#     <in.H is complex, and n3 is freq, n1 is traces
  !#     >out.H bigger
  !##################################################
  use sep
  use flags
  use timers
!  use omp_lib
!  use fft

  implicit none
  integer :: j1,i1,j2,i2,i3,j3,a,b,c,d,i,j,k1,k2,k3,ompi,ompn
  integer :: n1,n2,n3,n4,n5,esize,n1out
  real    :: o1,o2,o3,d1,d2,d3,maxf,maxt,mem
  real   , allocatable,dimension(:,:,:)  :: data
  complex, allocatable,dimension(:,:,:)  :: res,cdata,half
  complex, allocatable,dimension(:,:)  :: shottrace
  logical:: mig,time,center
  integer*8,pointer:: planF(:),planI(:)
  integer, external :: omp_get_thread_num
  integer, external :: omp_get_num_threads
  real:: startR,finishR,startW,finishW,startFFT,finishFFT, startMult, finishMult

  startR=0;finishR=0;startW=0;finishW=0;startFFT=0;finishFFT=0; startMult=0; finishMult=0
  call sep_init(SOURCE)
  call from_history('n1',n1);call from_history('n2',n2);!call from_history('n3',n3)
n3=1
  call from_history('o1',o1);call from_history('o2',o2);!call from_history('o3',o3)
  call from_history('d1',d1);call from_history('d2',d2);!call from_history('d3',d3)
  call from_history('esize',esize,4)

  call from_param('mig',mig,.false.);call from_param('center',center,.false.)
  call from_param('time',time,.false.); 

  call from_param('maxf',maxf,1./(2.*d1)); call from_param('maxt',maxt,4.)
  if(mig) then 
     center=.false. ; time=.false.
  end if
  n1out=n1;

  if(time) then
     mig=.false.
     call to_history('esize',4)
  else
     call to_history('esize',8)
  end if

  a=1; b=n2
  c=1; d=n3

  write(0,*) 'time,mig,center',time , mig, center
  if(center) then 
     a=n2/2 +1
     b=a 
     write(0,*) "Only constructing the center shot"
     c=n3/2 +1
     d=c
  end if
  write(0,*) 'a,b,c,d',a,b,c,d
  if (mig) then
     call to_history('n1',n2);call to_history('o1',o2);call to_history('d1',d2)
     call to_history('n2',n3);call to_history('o2',o3);call to_history('d2',d3)
     if (esize.eq.4) then
        write(0,*) 'Migrate ready: time in, freq out'
        n1out=min(n1,int(1.5+maxf*d1*n1));if(mod(n1out,2).ne.0) n1out=n1out+1
        call to_history('n3',n1out);call to_history('o3',o1);call to_history('d3',1/(n1*d1))
     else
        write(0,*) 'Migrate ready: freq in , freq out'
        n1out=min(n1,int(1.5+maxf/d1));
        call to_history('n3',n1out);call to_history('o3',o1);call to_history('d3',d1)
     end if
     call to_history('n4',b  );call to_history('o4',o2);call to_history('d4',d2)
     call to_history('n5',d  );call to_history('o5',o3);call to_history('d5',d3)
  else
     if (.not.time.and.esize.eq.4) then !freq out, time in
        n1out=min(n1,int(1.5+maxf*d1*n1))
        write(0,*)'freq out, time in'
        call to_history('n1',n1out);call to_history('o1',o1);call to_history('d1',1/(n1*d1))
     else if(.not.time.and.esize.eq.8) then !freq in/out
        n1out=min(n1,int(1.5+maxf/d1))
        write(0,*)'freq in/out'
        call to_history('n1',n1out);call to_history('o1',o1);call to_history('d1',d1)
     else if(esize.eq.8) then ! time out, freq in
        n1out=min(n1,int(1.5+maxt*n1*d1)) ; 
        d1=1/(n1*d1); 
        call to_history('n1',n1out);call to_history('o1',o1);call to_history('d1',d1)
        write(0,*)' freq in, time out',d1,n1,n1out
     else                 ! time in/out
        n1out=1.5+maxt/d1 ; write(0,*) 'time in/out'
        call to_history('n1',n1out);call to_history('o1',o1);call to_history('d1',d1)
     end if
     call to_history('n2',n2);call to_history('o2',o2);call to_history('d2',d2)
     call to_history('n3',n3);call to_history('o3',o2);call to_history('d3',d3)
     call to_history('n4',b-a+1);call to_history('o4',o2);call to_history('d4',d2)
     call to_history('n5',d-c+1);call to_history('o5',o3);call to_history('d5',d3)
     call to_history('label1','s');call to_history('label2','gx(m)');call to_history('label3','gy(m)')
     call to_history('label4','sx(m)');call to_history('label5','sy(m)')
  end if


  write(0,*) 'n3,n2,n1',n3,n2,n1
  write(0,*) 'd3,d2,d1',d3,d2,d1
  write(0,*) 'esize',esize
  
  allocate(cdata(n1,n2,n3))
  mem=mem+8*size(cdata)
  write(0,*) 'allocate cdata',mem/1024/1024
  allocate(res(n1,n2,n3));res=(0.,0.)
  mem=mem+8*size(res)
  write(0,*) 'allocate res',mem/1024/1024
  if(esize.eq.4) then
     allocate(data(n1,n2,n3))
     mem=mem+4*size(data)
     write(0,*) 'allocate data',mem/1024/1024

     call start_timer(startR)
     call sep_read(data)
write(0,*)' data',data(:,1,1)
     call end_timer(finishR,startR)
     do k3=1,n3
        do k2=1,n2
           do k1=1,n1
              res(k1,k2,k3)=cmplx(data(k1,k2,k3),0.)
           end do
        end do
     end do
write(0,*) 'cdata',res(:,1,1)
write(0,*)'fuck1',sum(cabs(res(:,1,1)))
     if (.not.time) then
        deallocate(data)
        mem=mem-4*size(data)
        write(0,*) 'deallocate data',mem/1024/1024
     end if
     allocate (planF(n2*n3))
     mem=mem+8*size(planF)
     write(0,*) 'allocate planF',mem/1024/1024
     call start_timer(startFFT)
write(0,*)' n3,n2,n1',n3,n2,n1
write(0,*)'size(planF)',size(planF)
write(0,*)'size(res,1),size(res,2),size(res,3)',size(res,1),size(res,2),size(res,3)
write(0,*)'size(data,1),size(data,2),size(data,3)',size(data,1),size(data,2),size(data,3)
write(0,*)'size(cdata,1),size(cdata,2),size(cdata,3)',size(cdata,1),size(cdata,2),size(cdata,3)
     do j=1,n3
        do i=1,n2
           k1=i+(j-1)*n2
!write(0,*) 'k1,i,j',k1,i,j
           call sfftw_plan_dft_1d(planF(k1),n1,&
           res(:,i,j),cdata(:,i,j),FFTW_FORWARD,FFTW_PATIENT)        
        end do
     end do

!call seperr("asdf")
!!$OMP PARALLEL DO PRIVATE(i)
     do i=1,n2*n3
write(0,*) 'i',i
        call sfftw_execute(planF(i))
     end do
!!$OMP END PARALLEL DO 
write(0,*)'fuck',sum(cabs(cdata(:,1,1))),sum(cabs(cdata(:,2,1)))
!write(0,*)' cdata'
!do i=1,n1
!write(0,*) i,cdata(i,1,1)
!end do
     call end_timer(finishFFT,startFFT)
     esize=8
     write(0,*) 'Finished data fft'
  else
     call sep_read(cdata)
  endif
  call sep_close()
write(0,*) 'asdf3'
!$OMP PARALLEL PRIVATE(ompi)
    ompn = omp_get_num_threads()  
!$OMP END PARALLEL 
write(0,*) 'ompn',ompn
    allocate(shottrace(n1,ompn))
    mem=mem+8*size(shottrace)
    write(0,*) 'allocate shottrace',mem/1024/1024
    shottrace=(0.,0.)

!    write(0,*) 'n1',n1,'n1out',n1out
!    allocate(half(n1out,n2,n3))
!    mem=mem+8*size(half)
!    write(0,*) 'allocate half',mem/1024/1024



    if(time) then
       allocate(planI(n2*n3))
       mem=mem+8*size(planI)
       write(0,*) 'allocate planI',mem/1024/1024

       do j=1,n3
          do i=1,n2
             k1=i+(j-1)*n2
             call sfftw_plan_dft_1d(planI(k1),n1,&
             res(1,i,j),res(1,i,j),FFTW_BACKWARD,FFTW_PATIENT)
          end do
       end do
    end if

!    write(0,'(8x,a,f10.3,a)')' *       Time reading       :', &
!      finishR,' (s)'
!    write(0,'(8x,a,f10.3,a)')' *       Time FFTing       :', &
!      finishFFT,' (s)'
write(0,*) 'before looping n3,n2,n1:',n3,n2,n1
  if(.not.mig) then
     shoty:do i3=c,d
        shotx:do i2=a,b
           finishMult=0.
           call start_timer(startMult)
           res=cmplx(0.,0.)
           shottrace=cmplx(0.,0.)
write(0,*)' sum(shotdata)',sum(cabs(cdata(:,i2,i3)))
!$OMP PARALLEL PRIVATE(ompi)
           ompi = omp_get_thread_num()+1    
           do k1=1,n1
              shottrace(k1,ompi)=conjg(cdata(k1,i2,i3))              
           end do
write(0,*) 'shot',sum(cabs(shottrace(:,ompi))),ompi
!$OMP DO PRIVATE(j2,j1,j3)
           y:do j3=1,n3
              x:do j2=1,n2
                 t:do j1=1,n1
                    res(j1,j2,j3)= shottrace(j1,ompi)*cdata(j1,j2,j3)
                 end do t
              end do x
           end do y
!$OMP END DO
!$OMP END PARALLEL
           call end_timer(finishMult,startMult)
           if(time) then
              do k3=1,n3
                 do k2=1,n2
                    j3=n2*(k3-1)+k2
!write(0,*) 'j3,k3,k2',j3,k3,k2
                    call sfftw_execute(planI(j3))
                    do k1=1,n1
                       data(k1,k2,k3)=real(res(k1,k2,k3))
                    end do
                 end do
              end do
              call srite('out',data,4*n1out*n2*n3)
           else
              do k3=1,n3
                 do k2=1,n2
                    do k1=1,n1out
                       half(k1,k2,k3)=res(k1,k2,k3)
                    end do
                 end do
              end do

              call start_timer(startW)
              call srite('out',half,8*n1out*n2*n3)
              call end_timer(finishW,startW)
           end if
!write(0,*) '         shot number:' ,i2
!           if (mod(i2,10)==0) then
!              write(0,*) 'Finished panel/cube',(i3-1)*i2+i2,'of',b*d
!              write(0,'(8x,a,f10.3,a)')' *       Time writing       :', &
!                  finishW,' (s)'
!              write(0,'(8x,a,f10.3,a)')' *       Time multing       :', &
!                  finishMult,' (s)'
!           end if
        end do shotx
     end do shoty
      write(0,'(8x,a,f10.3,a)')' *  Total Time writing       :', &
         finishW,' (s)'
      write(0,'(8x,a,f10.3,a)')' *  Total Time multing       :', &
         finishMult,' (s)'
  else
!    allocate(res(n2,n3,n1));res=(0.,0.)
!    shotym:do i3=c,d
!       shotxm:do i2=a,b
!          shottrace(:,ompi)=conjg(cdata(:,i2,i3))
!!$OMP PARALLEL DO PRIVATE(j3,j2,j1)
!          ym:do j3=1,n3
!             xm:do j2=1,n2
!                w:do j1=1,n1out
!                   res(j2,j3,j1)= shottrace(j1,ompi)*cdata(j1,j2,j3)
!                end do w
!             end do xm
!          end do ym
!!$OMP END PARALLEL DO
!          call srite('out',res(:,:,1:n1out) ,8*n1out*n2*n3)
!          if (mod(i2,10)==0)write(0,*) 'Finished panel',(i3-1)*i2+i2,'of',b*d
!       end do shotxm
!    end do shotym
 end if
write(0,*) 'deallocate'
  if(time) deallocate(data)
  deallocate(cdata,shottrace,res,half)
call auxclose("out")
call auxclose("in")
write(0,*) 'program complete'
end Program NewCorr
