  subroutine Rot(cube) 
    
    complex,allocatable					:: k_cube(:,:,:),k_temp(:),k_row(:)
    complex,allocatable					:: temp_z(:),temp_x_hdip(:,:),k_row(:)

    type (rotcig_param_type) rotcig_par
    integer status

    logical inside_k_hdip
    integer i_z,i_mx,i_hdip,i_k_hdip
    real geol_dip,off_dip,inp_dip,shift
    real*8 tan_geol_dip,tan_off_dip,tan_inp_dip,tan_out_dip
    real inv_d_khdip,k_hdip_inp,k_hdip_out,w_hdip

    inv_d_khdip=1./d_khdip

! zeroes k_cube
    k_cube=cmplx(0.,0.)
    k_temp=cmplx(0.,0.)
    k_row=cmplx(0.,0.)

  ! FFT along z
    do i_mx=1,kx%n
      do i_hdip=1,kh%n
        temp_z(1:kz%n)=cmplx(cube(1:kz%n,i_mx,i_hdip),0.)
        call fft_1d_data(fft_nz_id,FFT_BACK,kz%n,temp_z(:))
        do i_z=1,int(kz%n/2)+1
          k_cube(i_z,i_mx,i_hdip)=temp_z(i_z)
        end do
      end do 
    end do

  ! FFT along mx and hdip
    do i_z=1,int(kz%n/2)+1
      temp_x_hdip(:,:)=k_cube(i_z,:,:)
      call fft_2d_data(fft_nmx_nhdip_id,FFT_BACK,kx%n,kh%n,temp_x_hdip(:,:))
      k_cube(i_z,:,:)=temp_x_hdip(:,:)
    end do 

    do i_mx=1,n_mx_pad
      do i_z=1,int(kz%n/2)+1
        if(i_z /= 1) then
          tan_geol_dip=-kx%all(i_mx)/kz%all(i_z)
        else
          tan_geol_dip=sign_tan*sign(huge(tan_geol_dip),-kx%all(i_mx))
        end if

        do i_hdip=1,int(kh%n/2)+1
          k_temp(i_hdip+int(kh%n/2)-1)=k_cube(i_z,i_mx,i_hdip)
        end do
        do i_hdip=int(kh%n/2)+2,kh%n
          k_temp(i_hdip-int(kh%n/2)-1)=k_cube(i_z,i_mx,i_hdip)
        end do

        do i_hdip=1,kh%n
          shift=k_hdip(i_hdip)*hdip_0
          k_temp(i_hdip)=k_temp(i_hdip)*cexp(cmplx(0,shift))
        end do

          do i_hdip=1,kh%n
            if(i_z /= 1) then
              tan_off_dip=-k_hdip(i_hdip)/kz%all(i_z)
! Notice opposite sign in adding tangents
              if (abs(1+sign_tan*tan_off_dip*tan_geol_dip)>tiny(tan_off_dip)) &
                   then
                tan_inp_dip= &
                  (tan_geol_dip - sign_tan*tan_off_dip)/&
                  (1+sign_tan*tan_off_dip*tan_geol_dip)
  
                k_hdip_inp=-kz%all(i_z)*tan_inp_dip
              else
                k_hdip_inp=sign(kh%all(i_hdip),-kx%all(i_mx))
              end if
            else
                k_hdip_inp=sign(kh%all(i_hdip),-kx%all(i_mx))
            end if

            i_k_hdip=floor(k_hdip_inp*inv_d_khdip+float(int(kh%n/2)))
            inside_k_hdip=((i_k_hdip >= 1) .and.(i_k_hdip < kh%n))
  
            w_hdip=(k_hdip_inp*inv_d_khdip)+int(kh%n/2) - i_k_hdip
  
            if(inside_k_hdip) then
              k_row(i_hdip)= &
                           (1.-w_hdip)*k_temp(i_k_hdip  ) + &
                           (   w_hdip)*k_temp(i_k_hdip+1)
            else
              if(i_k_hdip == kh%n) then
                k_row(i_hdip)= &
                             (1.-w_hdip)*k_temp(i_k_hdip  ) + &
                             (   w_hdip)*k_temp(1         )
              else
                if(i_k_hdip == 0) then
                  k_row(i_hdip)= &
                             (1.-w_hdip)*k_temp(kh%n) + &
                             (   w_hdip)*k_temp(1         )
                end if
              end if
            end if
  
          end do 
        end if

        do i_hdip=1,kh%n
          shift=-kh%all(i_hdip)*hdip_0
          k_row(i_hdip)=k_row(i_hdip)*cexp(cmplx(0,shift))
        end do

        do i_hdip=1,int(kh%n/2)+1
          k_cube(i_z,i_mx,i_hdip)=k_row(i_hdip+int(kh%n/2)-1)
        end do
        do i_hdip=int(kh%n/2)+2,kh%n
          k_cube(i_z,i_mx,i_hdip)=k_row(i_hdip-int(kh%n/2)-1)
        end do

      end do 
    end do 


! INV FFT along mx and hdip

    do i_z=1,int(kz%n/2)+1
      temp_x_hdip(:,:)=k_cube(i_z,:,:)
      call fft_2d_data(fft_nmx_nhdip_id,FFT_FORW,n_mx_pad,kh%n,temp_x_hdip(:,:))

       k_cube(i_z,:,:)=temp_x_hdip(:,:)
    end do 

! INV FFT along z
    do i_mx=1,kx%n
      do i_hdip=1,kh%n
        do i_z=1,int(kz%n/2)+1
          temp_z(i_z)=k_cube(i_z,i_mx,i_hdip)
        end do
        do i_z=int(kz%n/2)+2,kz%n
          temp_z(i_z)= conjg(k_cube(kz%n-i_z+2,i_mx,i_hdip))
        end do

        call fft_1d_data(fft_nz_id,FFT_FORW,kz%n,temp_z(:))
        cube(1:kz%n,i_mx,i_hdip)=real(temp_z(1:kz%n))

      end do 
    end do 

  return	

  end subroutine Rot
