!!$
!!$
!!$=head1 NAME
!!$
!!$Transf - Transposes and FFTs a dataset
!!$
!!$=head1 SYNOPSIS
!!$
!!$Transf  <in.H >out.H
!!$
!!$TransfM  intag=in.H outtag=out.H
!!$
!!$=head1 INPUT PARAMETERS
!!$
!!$=over 4
!!$
!!$=item f_min - float
!!$
!!$        [0]: Minimum frequency (only needed for forward op)
!!$
!!$=item f_min1 - float
!!$
!!$        [f_min]: Minimum frequency with 100% pass (only needed for forward op)
!!$
!!$=item f_max - float
!!$
!!$        [.5/d1]: Maximum frequency  (only needed for forward op)
!!$
!!$=item f_max1 - float
!!$
!!$        [f_max]: Maximum frequency with 100% pass (only needed for forward op)
!!$
!!$=item maxsize - integer
!!$
!!$        [50]: Maximum megabytes to use
!!$
!!$=item shift2 - integer
!!$
!!$  [0]: How many n=1 axes to introduce or remove after the frequency axis
!!$
!!$
!!$=item verb  - logical
!!$
!!$        [n]: Whether or not to be verbose
!!$
!!$=item  wei  - logical
!!$
!!$        [false]: Do WEI style transforms
!!$
!!$=back
!!$
!!$
!!$=head1 DESCRIPTION
!!$
!!$Transforms a dataset (and its inverse) from (nt,nx,...) to (nx ...,nt)
!!$
!!$=head1 NOTES
!!$
!!$ Forward/Inverse determined from esize (4-forward, 8-inverse)
!!$
!!$ WARNING WARNING WARNING: Unpipable
!!$
!!$=cut
!!>
!!-#------------------------------------------------
!!
!!Author: Robert Clapp, ESMB 463, 7230253
!!
!!Date Created:Fri Sep 22 12:42:13 PDT 2000
!!
!!Purpose:
!!
program Transf 
  use sep
  use sep3d_struct_mod
  use cos_taper_mod
  implicit none
  integer :: i,ib,nt,i1,i2,iin1,iin2
  integer :: nb !block size
  integer :: shift2 !whether or not shift the second axis
  integer :: maxsize !maximummemory
  integer :: ir,ii   !imag and real pointers
  integer :: i_f_min,i_f_max
  integer :: nh !helix coordinate system
  integer :: nf !number of frequencies
  integer :: done,nbs
  logical :: verb,inv
  real    :: df !sampling of frequencies
  real    :: f_min,f_max
  real    :: dt,of,ot,w,f_min1,f_max1
  integer :: n,ndim,esize,np,ierr
  integer :: ng(2),fw(2),jw(2),nw(2),j
  character(len=128) :: label
  real, allocatable :: block_in(:,:),temp(:),block_out(:,:)
  complex,allocatable :: blockc_in(:,:),blockc_out(:,:),block_t(:),scalet(:),tt(:),block_cc(:,:)
  integer :: sign,iverb,nsect,impi,nmpi
  logical :: wei
  real :: pi
  real arg
  type(sep3d) ::  input,output,space
  complex, allocatable :: shift(:)
  character(len=20) :: intag,outtag
  integer, external :: init_loop_calc,do_sep_loop
  integer :: nc(5),nr(5),fc(5),fr(5),jc(5),jr(5),ndimc,ndimr,ns(5),fs(5)
  integer :: imap(5),n_c(2),n_r(2),fdone,ndone,stat
  integer, external :: sep_thread_num



  !----------------------------------------------------------------
  call sep_init()
call sep_add_doc_line("NAME")
call sep_add_doc_line("    Transf - Transposes and FFTs a dataset")
call sep_add_doc_line("")
call sep_add_doc_line("SYNOPSIS")
call sep_add_doc_line("    Transf <in.H >out.H")
call sep_add_doc_line("")
call sep_add_doc_line("    TransfM intag=in.H outtag=out.H")
call sep_add_doc_line("")
call sep_add_doc_line("INPUT PARAMETERS")
call sep_add_doc_line("    f_min - float")
call sep_add_doc_line("                [0]: Minimum frequency (only needed for forward op)")
call sep_add_doc_line("")
call sep_add_doc_line("    f_min1 - float")
call sep_add_doc_line("                [f_min]: Minimum frequency with 100% pass (only needed for forward op)")
call sep_add_doc_line("")
call sep_add_doc_line("    f_max - float")
call sep_add_doc_line("                [.5/d1]: Maximum frequency  (only needed for forward op)")
call sep_add_doc_line("")
call sep_add_doc_line("    f_max1 - float")
call sep_add_doc_line("                [f_max]: Maximum frequency with 100% pass (only needed for forward op)")
call sep_add_doc_line("")
call sep_add_doc_line("    maxsize - integer")
call sep_add_doc_line("                [50]: Maximum megabytes to use")
call sep_add_doc_line("")
call sep_add_doc_line("    shift2 - integer")
call sep_add_doc_line("          [0]: How many n=1 axes to introduce or remove after the frequency axis")
call sep_add_doc_line("")
call sep_add_doc_line("    verb - logical")
call sep_add_doc_line("                [n]: Whether or not to be verbose")
call sep_add_doc_line("")
call sep_add_doc_line("    wei - logical")
call sep_add_doc_line("                [false]: Do WEI style transforms")
call sep_add_doc_line("")
call sep_add_doc_line("DESCRIPTION")
call sep_add_doc_line("    Transforms a dataset (and its inverse) from (nt,nx,...) to (nx ...,nt)")
call sep_add_doc_line("")
call sep_add_doc_line("NOTES")
call sep_add_doc_line("     Forward/Inverse determined from esize (4-forward, 8-inverse)")
call sep_add_doc_line("")
call sep_add_doc_line("     WARNING WARNING WARNING: Unpipable")
call sep_add_doc_line("")

  !----------------------------------------------------------------
  intag="in"
  outtag="out"
  nmpi=1

  call doc(" ")
  call init_sep3d(intag,input,"INPUT")
  call from_param("maxsize",maxsize,50)
  call from_param("shift2",shift2,0)
  call from_param("verb",verb,.false.)
  call from_param("wei",wei,.false.)

  if(input%data_format(1:1)=="F") then
    inv=.false.
    call init_sep3d(output,"SCRATCH","COMPLEX","REGULAR",(/1,1,1,1,1/))
    ndimr=sep3d_ndims(input)
    call init_sep3d(output,space,"SCRATCH")
  else if(input%data_format(1:1)=="C") then
    inv=.true.
    call init_sep3d(output,"SCRATCH","FLOAT","REGULAR",(/1,1,1,1,1/))
    ndimc=sep3d_ndims(input)
    call init_sep3d(input,space,"SCRATCH")
  else
    call seperr("Only know how to handle complex and real data")
  end if


  imap=-1

  !TIme -> frequency
  if(.not. inv) then
    ndimc=sep3d_ndims(input)
    if (wei) ndimc=5
    if (verb .and. wei.and. sep_thread_num()==0) write(0,*) "Using  wei conventions"
    ot=input%o(1);dt=input%d(1);nt=input%n(1)
    nf=.5/dt;df=1./(dt*nt)
    call from_param("f_min",f_min,0.)
    call from_param("f_max",f_max,0.5/dt)
    call from_param("f_min1",f_min1,f_min+df*2)
    call from_param("f_max1",f_max1,f_max-df*3)
    of=0.
    i_f_min=max(1,nint(f_min/df))
    i_f_max=min(int(nt/2),nint(f_max/df))
    nf=i_f_max-i_f_min+1
    f_min = (i_f_min-1)*df
    f_max = (i_f_max-1)*df
    output%n(ndimc)=nf;output%o(ndimc)=f_min
    output%d(ndimc)=df;output%label(ndimc)="Frequency"
    imap(2)=1;
    do i=3,ndimr
     imap(i)=i-1+shift2
    end do
    space%n(ndimc)=nf; space%o(ndimc)=of; space%d(ndimc)=df;
    ndone=product(input%n(2:sep3d_ndims(input)))
  else
    ndimr=ndimc-shift2
    df=input%d(ndimc);of=input%o(ndimc);nf=input%n(ndimc)
    nt=1;do while (nf>nt);nt=nt*2;end do 
    nt=nt*2
    ot=0.
    dt=1./(nt*df)
    output%n(1)=nt; output%o(1)=ot; output%d(1)=dt; output%label(1)="Time"
    imap(1)=2
    do i=2,sep3d_ndims(input)-1
      imap(i)=i+1-shift2
    end do
    if(.not. sep3d_desection(space)) &
     call seperr("trouble desectioning space")
    ndone=product(input%n(1:sep3d_ndims(input)-1))
  end if

  do i=1,ndimc
    if(imap(i)/=-1) then 
      if(.not. inv) then
        space%n(imap(i))=input%n(i)
        space%o(imap(i))=input%o(i)
        space%d(imap(i))=input%d(i)
        space%label(imap(i))=input%label(i)
      end if
      output%n(imap(i))=input%n(i)
      output%o(imap(i))=input%o(i)
      output%d(imap(i))=input%d(i)
      output%label(imap(i))=input%label(i)
    end if
  end do
  call sep3d_set_sep3d(output)
  call sep3d_set_sep3d(space)
   
  nsect=1; if(.not. inv) nsect=nmpi
  call from_param("nsect",nsect,nsect)
  if(nsect>1) then
    if(.not. inv)  then
      if(.not. sep3d_section_tag(output,ndimc,"BLOCK",nsect)) &
         call seperr("trouble sectioning tag")
    else
      if(.not. sep3d_section_tag(space,1,"BLOCK",nsect)) &
         call seperr("trouble sectioning tag")
      if(.not. sep3d_section_tag(output,1,"BLOCK",nsect)) &
         call seperr("trouble sectioning tag")
    end if
  end if

  if (.not. inv) call init_cos_taper(nf,f_min,df,f_min,f_min1,f_max,f_max1)



  jw=1;nw=1;fw=0;ng=(/nh,nf/)

  allocate(scalet(nt))
  if (.not. inv) then
    w=8.*atan(1.) *ot
    do i=1,nt  
      scalet(i)=cexp(cmplx(0.,w*((i-1)*df+of)))
    end do
  else
    w=8.*atan(1.) *of
    do i=1,nt  
      scalet(i)=cexp(cmplx(0.,w*((i-1)*dt+ot)))
    end do
  end if 
!ALLOCATE MEMORY
  if (.not. inv) then
    allocate(temp(nt*2+16))
    allocate(tt(nf))
  else
    allocate(block_t(nt))         !scalet(nt)
  end if 
!INITIALZIE
  if (.not. inv) then
    call RFFTI(nt,temp)
  end if

call sep3d_write_description(outtag,output)


if(.not. inv) then
  if(0/=init_loop_calc(input%ndims,input%n,"MAIN",maxsize*1000*1000/12)) &
    call seperr("trouble initing loop")
else
  if(0/=init_loop_calc(output%ndims,output%n,"MAIN",maxsize*1000*1000/24)) &
    call seperr("trouble initing loop")
end if

 nc=1; nr=1;fc=0; fr=0;jc=1; jr=1;fs=0; ns=1;
 nr(1)=nt;
 nc(ndimc)=nf;
 ns(ndimc)=nf;
 fdone=0;
do while(0==do_sep_loop("MAIN",nr,fr)) !{

   !Figure out my windows
   if(.not. inv) then
     do i=2,ndimr
       if(imap(i)/=-1) then
          nc(imap(i))=nr(i)
          fc(imap(i))=fr(i)
       end if  
     end do
     if(.not. sep3d_local_buffer_size(output,n_c,nc,fc,jc) .or. &
        .not. sep3d_local_buffer_size(input,n_r,nr,fr,jr)) then
       call seperr("trouble getting local buffer sizes")
     end if
     allocate(block_in(n_r(1),product(nr(2:))),block_out(1,1))
     allocate(blockc_in(product(nr(2:)),nf))
     allocate(blockc_out(nc(1),product(nc(2:))))
     call sep3d_grab_headers(intag,input,nh,nwind=nr(2:),fwind=fr(2:))
     if(.not. sep3d_read_data(intag,input,block_in)) &
       call seperr("trouble reading in data")
   else 
     do i=1,ndimc
       if(imap(i)/=-1) then
           nc(i)=nr(imap(i))
           fc(i)=fr(imap(i))
!           nc(imap(i))=nr(i)
!           fc(imap(i))=fr(i)
       end if
     end do
     if(.not. sep3d_local_buffer_size(input,n_c,nc,fc,jc) .or. &
        .not. sep3d_local_buffer_size(output,n_r,nr,fr,jr)) then
       call seperr("trouble getting local buffer sizes")
     end if
     allocate(blockc_in(n_c(1),n_c(2)))
     allocate(blockc_out(n_r(2),nf))
     !read local headers and data
     call sep3d_grab_headers(intag,input,nh,nwind=nc(2:),fwind=fc(2:))
     if(.not. sep3d_read_data(intag,input,blockc_in,nwind=nc(1),fwind=fc(1))) &
       call seperr("trouble reading in data")
     call sep3d_set_window(space,nwind=nc,fwind=fc)
     if(.not. sep3d_transfer_data(input,space,blockc_in,blockc_out))&
       call seperr("trouble transfering data")
     allocate(block_out(n_r(1),product(nr(2:))),block_in(1,1))
   end if
    

   !DO THE TRANSFORM
   if(.not. inv) then!{
     do i=1,size(block_in,2)  !Transform to fourier space
        call RFFTF(nt,block_in(1,i),temp)
     end do
     do i=1,nf !{
       j=i+i_f_min-1
       if(j==1)then
          ir=1; ii=nt;
       else 
         ir=(j-1)*2;ii=ir+1
       end if
       if(wei) then
         do ib=1,n_r(2)
           blockc_in(ib,i)=scalet(j)*cmplx(block_in(ir,ib),-block_in(ii,ib))*cos_taper_filter(i)
         end do
       else
         do ib=1,nbs
!           blockc_in(ib,i)=scalet(j)*cmplx(block_in(ir,ib),-block_in(ii,ib))*cos_taper_filter(i)
           blockc_in(ib,i)=scalet(j)*cmplx(block_in(ir,ib),block_in(ii,ib))*cos_taper_filter(i)
         end do
       end if
     end do !}
   else !}{
     do i=1,n_r(2)!{
        if(wei) then
           block_t(1:nf)=cmplx(real(blockc_out(i,:)), -aimag(blockc_out(i,:)))
        else 
           block_t(1:nf)=cmplx(real(blockc_out(i,:)), -aimag(blockc_out(i,:)))
!           block_t(1:nf)=blockc_out(i,:)
        end if
        block_t(nf+1:nt)=cmplx(0.,0.)
        call ftu(1.,nt,block_t)
        block_out(:,i)=block_t!  *scalet
     end do!}
   end if!}

   !TRANSFER? and write out
   if(.not. inv) then
     call sep3d_set_window(output,nwind=nc, fwind=fc)
     call sep3d_set_window(space,nwind=nc, fwind=fc)
     if(.not. sep3d_transfer_data(space,output,blockc_in,blockc_out))&
       call seperr("trouble transfering data")
     if(.not. sep3d_write_data(outtag,output,blockc_out,nwind=nc, fwind=fc))&
       call seperr("trouble writing data")
   else 
     if(.not. sep3d_write_data(outtag,output,block_out,nwind=nr, fwind=fr))&
       call seperr("trouble writing data")
   end if
   fdone=fdone+n_r(2)
   deallocate(block_in,block_out,blockc_in,blockc_out)
   if(sep_thread_num()==0) write(0,*) "finished ",real(nint(real(fdone)/real(ndone)*1000.)/10.)," %"
end do !}
 

end program
