Program Corr
  !##################################################
  !#                                              
  !#  cross-correlate every trace with every other
  !#     <in.H is complex, and n1 is a power of 2
  !#     >out.H is real n1,n2,n2
  !##################################################
  use sep
  implicit none
  integer :: h,i,j,done,ix,iy,ierr,fcounter,bcounter
  integer :: n1,n2,n3,n4,n5,esize,maxmem,nslice,shot_x,shot_y
  real    :: o1,o2,o3,o4,o5,d1,d2,d3,d4,d5,slicemem
  complex :: shot
  complex, allocatable,dimension(:,:,:)  :: res
  complex, allocatable,dimension(:,:,:)  :: data
  complex, allocatable,dimension(:)  :: shottrace
  logical ::weiout
  call sep_init(SOURCE)
  call from_history('n1',n1);call from_history('n2',n2);call from_history('n3',n3)
  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);if (esize /= 8) call seperr(' esize 8 required ')
  call from_param('maxmem',maxmem,500)
  call from_param('weiout',weiout,.true.)
  call from_param('shot_x',shot_x);  call from_param('shot_y',shot_y)

  if (weiout) then
     ix=1.5+ (shot_x - o1)/d1
     iy=1.5+ (shot_y - o2)/d2
  else
     ix=1.5+ (shot_x - o2)/d2
     iy=1.5+ (shot_y - o3)/d3
  end if
  
  call to_history('n1',n1);call to_history('n2',n2);call to_history('n3',n3)
!;call to_history('n4',n2);call to_history('n5',n3)
  call to_history('o1',o1);call to_history('o2',o2);call to_history('o3',o3)
!;call to_history('o4',o2);call to_history('o5',o3)
  call to_history('d1',d1);call to_history('d2',d2);call to_history('d3',d3)
!;call to_history('d4',d2);call to_history('d5',d3)
  call to_history('esize',8)
  call sep_close()

  maxmem=maxmem*1024*1024
  slicemem=8 * n1*n2  !real 
  if(maxmem/slicemem <1) call seperr('  Increase maxmem, or input less data- less than one slice fits in memory.  ')
  nslice=min(maxmem/slicemem,real(n3))  !integer
  done=0
  write(0,*) '   Maximum memory to use',maxmem/1024/1024,'Mbytes'
  write(0,*) '   Memory per slice     ',max(floor(slicemem)/1024,1),'Kbytes'
  write(0,*) '   Total num slices     ',n3
  write(0,*) '   Slices per block     ',nslice
  write(0,*) '   Memory used          ',nslice*8*n1*n2/1024/1024,'Mbytes'
  if (weiout) then
     bcounter=0
     do while (done < n3)
        write(0,*) n3/nslice-done/nslice,'frequency block(s) to do.'
        bcounter=bcounter+1
        allocate(data(n1,n2,min(n3-done,nslice) ))
        done=done+nslice 
        call sreed('in', data, n1*n2*esize*size(data,3))
        write(0,*) 'Finished reading block',bcounter
        fcounter=0
        do h=1,size(data,3)
           fcounter=fcounter+1
           do i=1,n2
              do j=1,n1
                 shot=conjg(data(ix,iy,h))
                 data(j,i,h)= data(j,i,h) * shot
              end do
           end do
           if (mod(h,size(data,3)/10)==0) write(0,*) 'Finished',100*fcounter/size(data,3),'% of block',bcounter
        end do
        call sep_write(data)
        deallocate(data)
     end do
  else
     allocate (data(n1,n2,n3))
     allocate (shottrace(n1))
     call sseek_block('in',n2*(iy-1)+ix ,esize*n1,0)
     call sreed('in',shottrace,esize*n1)
     call sseek('in',0,0)
     call sep_read(data)
     do h=1,n3
        shot=conjg(shottrace(h))
        do j=1,n2
           do i=1,n1
              data=data(i,j,h)*shot 
           end do
        end do
     end do
     call sep_write(data)
     deallocate(data)
     deallocate(shottrace)
  end if

end Program Corr











