module wavelet

  use precision
  use rationnalmod
  use var_glo
  use init
  use fft_1d
  use io
  use integral

  implicit none

contains

  subroutine get_time (t1, t2, r)
    integer, intent(IN) :: t1, t2, r
    integer             :: h, m, s, t

    t = int (real(t2-t1)/real(r))

    if (t<60) then
       print *, t, " secondes"
    else
       s = modulo (t, 60)
       m = t / 60
       if (m<60) then
          print *, m, " m. ", s, " s."
       else
          h = m/60
          m = modulo (m, 60)
          print *, h, " h. ", m, " m. ", s, " s."
       end if
    end if
    
  end subroutine get_time

  !! Returns the the value of the filter
  !! associated to the spline of order "order"
  function wavelet_filter (chi, order)
    real(PR), intent(IN) :: chi
    integer,  intent(IN) :: order
    complex(PR)          :: wavelet_filter, val
    integer              :: i

    val = exp(II*chi*real(ishft (order, -1), PR))
    do i=1, order
       val = val*(1.0d0 + exp(-II*chi))
    end do

    i = ishft (1, order)
    wavelet_filter = val/real(i, PR)

  end function wavelet_filter

  !! Returns the the value of the filter
  !! associated to the dual spline of order "(order, order_dual)"
  function wavelet_filter_dual (chi, order, order_dual)
    real(PR), intent(IN) :: chi
    integer,  intent(IN) :: order, order_dual
    complex(PR)          :: wavelet_filter_dual, val, val_sin
    integer              :: i, j, ordre_tot
    
    
    ordre_tot = order + order_dual
    
    val = dcmplx(0.d0, 0.0d0)

    do i=0, (ordre_tot/2 - 1)
       val_sin = dcmplx(1.0d0, 0.0d0)
       do j=1, 2*i
          val_sin = val_sin*sin(chi/2.0d0)
       end do
       val = val + real(binomial(ordre_tot/2 - 1 + i, i), PR)*val_sin
    end do

    if (btest(order_dual, 0)) then
       val = val*exp (-II*chi/2.0d0)
    end if

    do i=1, order_dual
       val = val*cos(chi/2.0d0)
    end do

    wavelet_filter_dual = val

  end function wavelet_filter_dual

  subroutine wavelet_display_filter_coef (tab)
    type(rationnal), dimension(:), intent(OUT) :: tab
    integer                                    :: i

    do i=1, size(tab, 1)
       call rat_disp (tab(i))
    end do

  end subroutine wavelet_display_filter_coef

  !! Returns the coefficient of the filter
  !! associated to the spline of order "order"
  subroutine wavelet_filter_get_coef (order, tab)
    integer,                       intent(IN)  :: order
    type(rationnal), dimension(:), intent(OUT) :: tab
    integer(print)                             :: puiss
    integer                                    :: i
    
    do i =  -ishft(order, -1), order - ishft(order, -1)
       puiss = ishft(1, order)
       tab (i+ishft(order, -1)+1) = rat_constr (binomial(order, i+ishft(order, -1)), puiss)
    end do

  end subroutine wavelet_filter_get_coef

  !! Returns the coefficient of the filter
  !! associated to the dual spline of order "(order, order_dual)"
  subroutine wavelet_filter_dual_get_coef (order, order_dual, tab)
    integer,                       intent(IN)  :: order, order_dual
    type(rationnal), dimension(:), intent(OUT) :: tab
    type(rationnal)                            :: coef1, coef2, coef3
    integer(print)                             :: puiss
    integer                                    :: j, l, n, order_tot, i
    
    
    order_tot = order + order_dual

    tab = rat_constr (0_print, 1_print)
    
    do j= -ishft(order_dual, -1), order_dual-ishft(order_dual, -1)
       puiss = ishft(1, order_dual)
       coef1 = rat_constr (binomial(order_dual, j+ishft(order_dual, -1)), puiss)
       do n=0, ishft(order_tot, -1)-1
          puiss = ishft(1, 2*n)
          coef2 = rat_constr (binomial(ishft(order_tot, -1)-1+n, n), puiss)
          coef2 = rat_mul (coef1, coef2)
          do l = -n, n
             i = j + l + ishft(order+order_dual, -1)+ishft(order_dual, -1)
             coef3 = rat_constr (binomial(2*n, n+l), 1_print)
             if (btest(l, 0)) then
                coef3 = rat_opp (coef3)
             end if
             coef3 = rat_mul (coef2, coef3)
             tab(i) = rat_add (tab(i), coef3)
          end do
       end do
    end do

  end subroutine wavelet_filter_dual_get_coef

  !! Returns the fourier function of the scale function
  function wavelet_fourier_scale_fct (chi, order)
    real(PR), intent(IN) :: chi
    integer,  intent(IN) :: order
    complex(CPR)         :: val, wavelet_fourier_scale_fct
    integer              :: i

    val = exp (II*chi*real(ishft (order,-1), PR))

    if (chi <> 0.0d0) then
       do i=1, order
          val = val * (1.0d0-exp(-II*chi))/(II*chi)
       end do
    end if

    wavelet_fourier_scale_fct = val

  end function wavelet_fourier_scale_fct

  !! Returns the scale function
  subroutine wavelet_get_scale_fct (order, scale, scale_fct)
    integer,                    intent(IN)  :: order, scale
    complex(CPR), dimension(:), intent(OUT) :: scale_fct
    complex(CPR), dimension(:), allocatable :: Fphi
    complex(CPR), dimension(:), allocatable :: wsave
    real(PR)                                :: z, offset
    integer                                 :: N, i, period

    period = ishft(order, scale)
    N = ishft(Nx, scale)

    offset = real (ishft (order, -1), PR)

    allocate (Fphi(N))
    
    do i=1, N
       z = Pi*real(i-1-N/2, PR)/real(period, PR)
       Fphi(i) = wavelet_filter(z, order)*wavelet_fourier_scale_fct(z, order)
       Fphi(i) = Fphi(i)*exp(-offset*2.0d0*II*z)/real(period, PR)
    end do

    allocate (wsave(4 * N + 15))
    call zffti (N, wsave)

    do i=1, ishft (N, -1)
       scale_fct(i)     = Fphi(i+ishft (N, -1))
       scale_fct(i+ishft (N, -1)) = Fphi(i)
    end do

    call zfftb (N, scale_fct, wsave)

    deallocate (wsave, Fphi)

  end subroutine wavelet_get_scale_fct

  !! Returns the mother wavelet function
  subroutine wavelet_get_mother_wavelet (order, order_dual, scale, mwave)
    integer,                    intent(IN)  :: order, order_dual, scale
    complex(CPR), dimension(:), intent(OUT) :: mwave
    complex(CPR), dimension(:), allocatable :: Fphi
    complex(CPR), dimension(:), allocatable :: wsave
    real(PR)                                :: z1
    integer                                 :: N, i, k, period

    !!period = ishft(order + order_dual - 1, scale)
    period =(order + order_dual - 1)
    N = ishft(Nx, scale)


    allocate (Fphi(N))
    
    do i=1, N
       k = i-1-N/2
       z1 = Pi*real(k, PR)/real(period, PR)
       Fphi(i) = exp(II*z1)
       Fphi(i) = Fphi(i)*conjg(wavelet_filter_dual(z1+Pi, order, order_dual))
       Fphi(i) = Fphi(i)*wavelet_fourier_scale_fct(z1, order)
       Fphi(i) = Fphi(i)*exp(-II*z1)*exp(-II*Pi*real(k, PR))
       !Fphi(i) = conjg(Fphi(i))
    end do

    Fphi = Fphi/real(period, PR)

    allocate (wsave(4 * N + 15))
    call zffti (N, wsave)

    do i=1, ishft (N, -1)
       mwave(i)     = Fphi(i+ishft (N, -1))
       mwave(i+ishft (N, -1)) = Fphi(i)
    end do

    call zfftb (N, mwave, wsave)

    deallocate (wsave, Fphi)

  end subroutine wavelet_get_mother_wavelet

  subroutine wavelet_get_fct_from_base (order, coef, phi, fct)
    integer,                    intent(IN)  :: order
    real(CPR),    dimension(:), intent(IN)  :: coef
    complex(CPR), dimension(:), intent(IN)  :: phi
    complex(CPR), dimension(:), intent(OUT) :: fct
    complex(CPR)                            :: z
    real(PR)                                :: x, y
    integer                                 :: i, j, k, N, period

    N = size (phi, 1)
    period = order*N/Nx

    do i=1, Nx
       x = real((i-1)*order, PR)/real(Nx, PR) - real(ishft (order, -1), PR)
       z = dcmplx (0.d0, 0.d0)
       do j=0, order
          y = 2.0d0*x - real(j-2*ishft (order, -1), PR)
          k = int(y*real(N, PR)/real(period, PR))
          if ((k>=1) .and. (k<=N)) then
             z = z + coef(j+1)*phi(k)
          end if
       end do
       fct(i) = z
    end do

  end subroutine wavelet_get_fct_from_base

  subroutine wavelet_get_fixed_point (N1, N2, coef_filter, scale, fixed_point)
    integer,                       intent(IN)  :: N1, N2, scale
    type(rationnal), dimension(:), intent(IN)  :: coef_filter
    type(rationnal), dimension(:), intent(OUT) :: fixed_point
    type(rationnal), dimension(:), allocatable :: tab, s, coefs
    type(rationnal)                            :: val, val1
    integer                                    :: j, k, l, n, kk
    integer                                    :: k_tab, n_tab, period
    character(len=100)                         :: filename
    real(PR),       dimension(:), allocatable  :: coef

    period = N2-N1
    allocate (s(period))

    ! il faut *2 les coef du filtres
    n = N2 - N1 + 1
    allocate (coefs(n))
    do k = 1, n
       coefs(k) = rat_mul (coef_filter(k), rat_constr (2_print, 1_print))
    end do

    ! calcul de s0

    do k=N1, N2-1
       if ((k>-0.5) .and. (k<0.5)) then
          s(k-N1+1) = rat_constr (1_print, 1_print)
       else
          s(k-N1+1) = rat_constr (0_print, 1_print)
       end if
    end do

    ! calcul des suivants

    do j=1, scale
       period = ishft (period, 1)
       allocate (tab(period))
       do k_tab = 0, period-1
          k = N1*ishft (1, j) + k_tab
          val = rat_constr (0_print, 1_print)

!!$          do n = ((k-N2)/2-1), ((k-N1)/2+1)
!!$             if ((k>=N1+2*n).and.(k<=N2+2*n).and.(n>=N1*ishft(1,j-1)).and.(n<=N2*ishft(1,j-1)-1)) then
!!$                val1 = rat_mul (coefs(k-2*n-N1+1), s(n-N1*ishft(1,j-1)+1))
!!$                val = rat_add (val, val1)
!!$             end if
!!$          end do

          do n_tab=1, size (s)
             n = n_tab-1+ishft(N1,j-1)
             l = k - ishft(n, 1)
             if ((l>=N1) .and. (l<=N2)) then
                val1 = rat_mul (coefs(l-N1+1), s(n_tab))
                val = rat_add (val, val1)
             end if
          end do

          tab(k_tab+1) = val

!!$          if(j==4) then
!!$             print *, "* k :", k, val%num,"/",val%den
!!$          end if

       end do
       deallocate (s)
       allocate (s(period))
       do n = 1, period
          s(n) = tab(n)
       end do
       deallocate (tab)

!!$       if(j==6) then
!!$          allocate(coef(size(s)))
!!$          do n=1, size(s)
!!$             coef(n) = real(s(n)%num)/real(s(n)%den)
!!$          end do
!!$          filename="fct_ech"
!!$          call save_tab(dcmplx(coef),filename);
!!$          deallocate(coef)
!!$       end if

    end do

    fixed_point = s
    
!!$    
!!$    do n=1, Nx
!!$       fixed_point(n) = s((n-1)*(N2-N1)+1)
!!$    end do

    deallocate (coefs, s)

  end subroutine wavelet_get_fixed_point

  subroutine wavelet_get_dual_scale_fct (order, order_dual, j, res)
    integer,                       intent(IN)  :: order, order_dual, j
    type(rationnal), dimension(:), intent(OUT) :: res
    type(rationnal), dimension(:), allocatable :: coef_filter
    integer                                    :: N1, N2, i, scale
    character(len=100)                         :: filename

    scale = 0
    do while (ishft (1, scale) < Nx)
       scale = scale + 1
    end do

    N1 = -(order+order_dual)/2-order_dual/2+1
    N2 = order_dual-order_dual/2+(order+order_dual)/2-1
    allocate (coef_filter(order+2*order_dual-1))

    call wavelet_filter_dual_get_coef (order, order_dual, coef_filter)

    print *, "Coef de phi dual"
    call wavelet_display_filter_coef (coef_filter)

    print *, "Debut fixed"
    call wavelet_get_fixed_point (N1, N2, coef_filter, scale+j, res)
    print *, "Fin fixed"

    filename="fct_ech_dual"
    call save_tab(res,filename);

  end subroutine wavelet_get_dual_scale_fct

  subroutine wavelet_get_scale_fct_2 (order, j, res)
    integer,                       intent(IN)  :: order, j
    type(rationnal), dimension(:), intent(OUT) :: res
    integer                                    :: N1, N2, i, scale

    scale = 0
    do while (ishft (1, scale) < Nx)
       scale = scale + 1
    end do

    N1 = -(order)/2
    N2 = order-order/2

!!$    print *, "Coef de phi"
!!$    call wavelet_display_filter_coef (coef_filter)
    
    call wavelet_get_fixed_point (N1, N2, filter_coef, scale+j, res)

  end subroutine wavelet_get_scale_fct_2

  function wavelet_scale_fct_get_value_r (order, val)
    integer,  intent(IN) :: order
    real(PR), intent(IN) :: val
    real(PR)             :: wavelet_scale_fct_get_value_r
    real(PR)             :: res1, res2
    integer              :: i, j, k

    if ((val < -real(ishft (order, -1), PR)) .or. (val >= real(order-ishft (order, -1), PR))) then
       wavelet_scale_fct_get_value_r = 0.d0
    else
       k = int (val+real(ishft (order, -1), PR)) + 1

       res1 = 0.d0
       do i=1, order
!!$          res2 = real(pol_spline (k, i)%num, PR)/real(pol_spline (k, i)%den, PR)
          res2 = pol_spline_real (k, i)
          do j=1, i-1
             res2 = res2 * val
          end do
          res1 = res1 + res2
       end do

!!$       res1 = real(pol_spline (k, order)%num, PR)/real(pol_spline (k, order)%den, PR)
!!$       res1 = pol_spline_real (k, order)
!!$       do i=1, order-1
!!$          res1 = res1 * val
!!$          res1 = res1 + real(pol_spline (k, -i+order)%num, PR)/real(pol_spline (k, -i+order)%den, PR)
!!$          res1 = res1 + pol_spline_real (k, -i+order)
!!$       end do

       wavelet_scale_fct_get_value_r = res1
    end if

  end function wavelet_scale_fct_get_value_r

  function wavelet_scale_fct_get_value_rat (order, val)
    integer,         intent(IN) :: order
    type(rationnal), intent(IN) :: val
    type(rationnal)             :: wavelet_scale_fct_get_value_rat
    type(rationnal)             :: res1, res2, N1, N2
    integer                     :: i, j, k

    N1 = rat_constr (-int (ishft (order, -1), print), 1_print)
    N2 = rat_constr (int (order - ishft (order, -1), print), 1_print)
    if ((rat_greater (N1, val) == TRUE) .or. (rat_greater_or_equal (val, N2) == TRUE)) then
       wavelet_scale_fct_get_value_rat = rat_constr (0_print, 1_print)
    else
       res1 = rat_constr (0_print, 1_print)
       res2 = rat_sub (val, N1)
       k = res2%num/res2%den + 1
       do i=1, order
          res2 = pol_spline (k, i)
          do j=1, i-1
             res2 = rat_mul (res2, val)
          end do
          res1 = rat_add (res1, res2)
       end do
       wavelet_scale_fct_get_value_rat = res1
    end if
    
  end function wavelet_scale_fct_get_value_rat

  !! Moments

  function moment(a, b, p, fct)
    real(PR),                   intent(IN)  :: a, b
    integer,                    intent(IN)  :: p
    complex(CPR), dimension(:), intent(IN)  :: fct
    complex(CPR), dimension(:), allocatable :: fct2
    complex(CPR)                            :: moment
    complex(CPR)                            :: val
    real(PR)                                :: dx, x, x_p
    integer                                 :: N, i, k

    N = size(fct, 1)
    allocate (fct2(N))
    dx = (b-a)/real(N, PR)

    val = dcmplx (0.d0, 0.0d0)
    do i=0, N-1
       x = a + dx*real(i, PR)
       x_p = 1.0d0
       do k=1, p
          x_p = x_p*x
       end do
!!$       fct2(i+1) = x_p*fct(i+1)
       val = val + x_p*fct(i+1)*dx
    end do

!!$    call int_romberg (fct2, b-a, N, val)
    moment = val

  end function moment
  
  function wavelet_moment_scale_fct (order, moment)
    integer,                     intent(IN) :: order
    integer,                     intent(IN) :: moment
    real(PR)                                :: wavelet_moment_scale_fct, val
    complex(CPR), dimension(:), allocatable :: phi
    real(PR)                                :: x, x_p
    integer                                 :: i, k

    allocate (phi(Nx))

    call wavelet_get_scale_fct (order, 0, phi)

    val = 0.d0
    do i = 1, Nx
       x = real((i-1)*order, PR)/real(Nx, PR)-real(ishft(order, -1), PR)
       x_p = 1.d0
       do k=1, moment
          x_p = x_p*x
       end do
       val = val + x_p*phi(i)*real(order, PR)/real(Nx, PR) 
    end do

    wavelet_moment_scale_fct = val

  end function wavelet_moment_scale_fct

  function wavelet_moment_scale_fct_2 (order, moment)
    integer,                        intent(IN) :: order
    integer,                        intent(IN) :: moment
    type(rationnal), dimension(:), allocatable :: tab_moment
    type(rationnal), dimension(:), allocatable :: coef_filter
    type(rationnal)                            :: wavelet_moment_scale_fct_2
    type(rationnal)                            :: coef1, coef2, coef3
    type(rationnal)                            :: val, val2
    integer                                    :: N1, N2
    integer                                    :: k, l, n, j
    integer(print)                             :: n_puiss

    N1 = -(order)/2
    N2 = order-order/2
    allocate (coef_filter(order+1), tab_moment(moment+1))

    call wavelet_filter_get_coef (order, coef_filter)

    do k=1, order+1
       coef_filter(k) = rat_mul (coef_filter(k), rat_constr(2_print, 1_print))
    end do

    tab_moment(1) = rat_constr (1_print, 1_print)
    print *, ""
    print *, "Momo : ", 0, tab_moment(1)%num, "/", tab_moment(1)%den

    do k = 1, moment
       n_puiss = ishft((ishft(1, k)-1),1)
       coef1 = rat_constr(1_print, n_puiss)
       val = rat_constr(0_print, 1_print)
       do l = 0, k-1
          val2 = rat_constr (binomial(k, l), 1_print)
          coef2 = rat_mul (val2, tab_moment(l+1))
          val2 = rat_constr (0_print, 1_print)
          do n = -ishft(order, -1), order-ishft(order, -1)
             n_puiss = 1_print
             do j = 1, k-l
                n_puiss = n_puiss*n
             end do
             coef3 = rat_constr (n_puiss, 1_print)
             coef3 = rat_mul (coef_filter(n+ishft(order, -1)+1), coef3)
             val2 = rat_add (val2, coef3)
          end do
          val = rat_add(val, rat_mul (val2, coef2))
       end do
       tab_moment(k+1) = rat_mul (val, coef1)
       print *, "Momo : ", k, tab_moment(k+1)%num, "/", tab_moment(k+1)%den
    end do

    wavelet_moment_scale_fct_2 = tab_moment(moment+1)

  end function wavelet_moment_scale_fct_2

  subroutine wavelet_get_moments (N1, N2, coef_filter, moment, tab_moment)
    integer,                       intent(IN)  :: N1, N2
    integer,                       intent(IN)  :: moment
    type(rationnal), dimension(:), intent(IN)  :: coef_filter
    type(rationnal), dimension(:), intent(OUT) :: tab_moment
    type(rationnal), dimension(:), allocatable :: coefs
    type(rationnal)                            :: coef1, coef2, coef3
    type(rationnal)                            :: val, val2
    integer                                    :: k, l, n, j
    integer(print)                             :: n_puiss

    n = size(coef_filter, 1)
    allocate (coefs(n))

    do k=1, n
       coefs(k) = rat_mul (coef_filter(k), rat_constr(2_print, 1_print))
    end do

    tab_moment(1) = rat_constr (1_print, 1_print)

    do k = 1, moment
       n_puiss = ishft((ishft(1, k)-1),1)
       coef1 = rat_constr(1_print, n_puiss)
       val = rat_constr(0_print, 1_print)
       do l = 0, k-1
          val2 = rat_constr (binomial(k, l), 1_print)
          coef2 = rat_mul (val2, tab_moment(l+1))
          val2 = rat_constr (0_print, 1_print)
          do n = N1, N2
             n_puiss = 1_print
             do j = 1, k-l
                n_puiss = n_puiss*n
             end do
             coef3 = rat_constr (n_puiss, 1_print)
             coef3 = rat_mul (coefs(n-N1+1), coef3)
             val2 = rat_add (val2, coef3)
          end do
          val = rat_add(val, rat_mul (val2, coef2))
       end do
       tab_moment(k+1) = rat_mul (val, coef1)
    end do

    deallocate (coefs)

  end subroutine wavelet_get_moments

  subroutine wavelet_moments_scale_fct (order, moment, tab_moment)
    integer,                       intent(IN)  :: order
    integer,                       intent(IN)  :: moment
    type(rationnal), dimension(:), intent(OUT) :: tab_moment
    type(rationnal), dimension(:), allocatable :: coef_filter
    integer                                    :: N1, N2

    N1 = -(order)/2
    N2 = order-order/2
    allocate (coef_filter(order+1))

    call wavelet_filter_get_coef (order, coef_filter)

    call wavelet_get_moments (N1, N2, coef_filter, moment, tab_moment)

  end subroutine wavelet_moments_scale_fct

  subroutine wavelet_moments_dual_scale_fct (order, order_dual, moment, tab_moment)
    integer,                       intent(IN)  :: order, order_dual
    integer,                       intent(IN)  :: moment
    type(rationnal), dimension(:), intent(OUT) :: tab_moment
    type(rationnal), dimension(:), allocatable :: coef_filter
    integer                                    :: N1, N2


    N1 = -(order+order_dual)/2-order_dual/2+1
    N2 = order_dual-order_dual/2+(order+order_dual)/2-1
    allocate (coef_filter(order+2*order_dual-1))

    call wavelet_filter_dual_get_coef (order, order_dual, coef_filter)

    call wavelet_get_moments (N1, N2, coef_filter, moment, tab_moment)

  end subroutine wavelet_moments_dual_scale_fct
  
  function wavelet_moment_mother_wavelet (order, order_dual, moment)
    integer,                     intent(IN) :: order, order_dual
    integer,                     intent(IN) :: moment
    real(PR)                                :: wavelet_moment_mother_wavelet, val
    complex(CPR), dimension(:), allocatable :: psi
    real(PR)                                :: x, x_p
    integer                                 :: i, k

    allocate (psi(Nx))

    call wavelet_get_mother_wavelet (order, order_dual, 0, psi)

    val = 0.d0
    do i = 1, Nx
       x = real((i-1)*order, PR)/real(Nx, PR)-real(ishft(order, -1), PR)
       x_p = 1.d0
       do k=1, moment
          x_p = x_p*x
       end do
       val = val + x_p*psi(i)*real(order, PR)/real(Nx, PR) 
    end do

    wavelet_moment_mother_wavelet = val

  end function wavelet_moment_mother_wavelet

  function wavelet_moment_mother_wavelet_2 (order, order_dual, moment)
    integer,                        intent(IN) :: order, order_dual
    integer,                        intent(IN) :: moment
    type(rationnal), dimension(:), allocatable :: tab_moment
    type(rationnal), dimension(:), allocatable :: coef_filter
    type(rationnal)                            :: wavelet_moment_mother_wavelet_2
    type(rationnal)                            :: coef1, coef2
    type(rationnal)                            :: val, val2
    integer                                    :: N1, N2
    integer                                    :: k, l, n, j
    integer(print)                             :: n_puiss

    N1 = -(order+order_dual)/2-order_dual/2+1
    N2 = order_dual-order_dual/2+(order+order_dual)/2-1
    allocate (coef_filter(order+2*order_dual-1), tab_moment(moment+1))

    call wavelet_filter_dual_get_coef (order, order_dual, coef_filter)

    call wavelet_moments_scale_fct (order, moment, tab_moment)

    val = rat_constr (0_print, 1_print)

    do j=1, order+2*order_dual-1
       n =  - N1 - j
       if (btest (n, 0)) then
          coef1 = coef_filter(j)
       else
          coef1 = rat_opp (coef_filter(j))
       end if
       val2 = rat_constr (0_print, 1_print)
       do l=0, moment
          n_puiss = 1
          do k=1, moment-l
             n_puiss = n_puiss*n
          end do
          n_puiss = n_puiss*binomial(moment,l)
          coef2 = rat_constr (n_puiss, 1_print)
          coef2 = rat_mul (coef2, tab_moment(l+1))
          val2 = rat_add (val2, coef2)
       end do
       val = rat_add (val, rat_mul (coef1, val2))
    end do
    
    n_puiss = ishft (1, moment)
    coef1 = rat_constr (1_print, n_puiss)
    
    val = rat_mul (val, coef1)

    wavelet_moment_mother_wavelet_2 = val

  end function wavelet_moment_mother_wavelet_2

  !! calcul des coef du polynome de Lagrange
  !! aux point 0, ..., l (l=degre du pol)
  
  subroutine get_lagrange_coef (degre, indice, coef)
    integer,                       intent(IN)  :: degre, indice
    type(rationnal), dimension(:), intent(OUT) :: coef
    integer(print), dimension(:),  allocatable :: tab
    integer(print)                             :: den
    integer                                    :: i, j, k

    !! denominateur
    den = 1_print
    do k = 0, degre
       if (k .ne. indice) then
          den = den*(indice-k)
       end if
    end do

    !! coefs du polynome

    allocate (tab(degre+1))
    tab = 0_print
    tab(1) = 1_print
    i = 0
    do k = 0, degre
       if (k .ne. indice) then
          i = i+1
          tab(i+1) = tab(i)
          do j = 2, i
             tab(i-j+2) = -k*tab(i-j+2) + tab(i-j+1)
          end do
          tab(1) = -k*tab(1)
       end if
    end do

    do i = 1, degre+1
       coef(i) = rat_constr (tab(i), den)
    end do

  end subroutine get_lagrange_coef

  !! calcul des coef du polynome de Lagrange
  !! aux point -[l/2], ..., l-[l/2]-1 (l = degre du pol)
  
  subroutine get_lagrange_coef_general (degre, indice, coef)
    integer,                       intent(IN)  :: degre, indice
    type(rationnal), dimension(:), intent(OUT) :: coef
    integer(print), dimension(:),  allocatable :: tab
    integer(print)                             :: den
    integer                                    :: j, k, n
    integer(print)                             :: n_inf, n_sup

    n_inf = -ishft (degre, -1)
    n_sup = degre - ishft (degre, -1) - 1

    !! denominateur
    den = 1_print
    do k = n_inf, n_sup
       if (k .ne. indice) then
          den = den*(indice - k)
       end if
    end do

    !! coefs du polynome

    allocate (tab(degre))
    if (size (tab) .ne. size (coef)) then
       print *, "Erreur dans la taille !! (coef Lagrange general)"
       stop
    end if
    tab = 0_print
    tab(1) = 1_print
    n = 0
    do k = n_inf, n_sup
       if (k .ne. indice) then
          n = n + 1
          tab(n+1) = 1_print
          do j = 2, n
             tab(n-j+2) = tab(n-j+1) - k*tab(n-j+2)
          end do
          tab(1) = -k*tab(1)
       end if
    end do

    do k = 1, degre
       coef(k) = rat_constr (tab(k), den)
    end do

  end subroutine get_lagrange_coef_general

  function wavelet_innerproduct_lagrange (order, fct, degre, scale)
    integer,                       intent(IN)  :: order, degre, scale
    real(PR),        dimension(:), intent(IN)  :: fct
    real(PR)                                   :: wavelet_innerproduct_lagrange
    type(rationnal), dimension(:), allocatable :: lagrange, moments
    type(rationnal)                            :: w
    real(PR)                                   :: c
    integer                                    :: k, n

    allocate (lagrange(degre+1), moments(degre+1))

    call wavelet_moments_scale_fct (order, degre, moments)

    c = 0.d0
    do n=0, degre
       
       call get_lagrange_coef (degre, n, lagrange)
       w = rat_constr (0_print, 1_print)
       do k = 1, degre+1
          w = rat_add (w, rat_mul (lagrange(k), moments(k)))
       end do
       
       c = c + fct(n+1)*real(w%num, PR)/real(w%den, PR)
       
    end do
    
    wavelet_innerproduct_lagrange = c/sqrt(real(ishft (1, scale), PR))

    deallocate (lagrange, moments)

  end function wavelet_innerproduct_lagrange

  function wavelet_ip_dual_lagrange (order, order_dual, fct, degre, scale)
    integer,                       intent(IN)  :: order, order_dual, degre, scale
    real(PR),        dimension(:), intent(IN)  :: fct
    real(PR)                                   :: wavelet_ip_dual_lagrange
    real(PR)                                   :: c
    integer                                    :: n


    c = 0.d0
    do n=0, degre
       c = c + fct(n+1)*dual_weight(n+1)
    end do
    
    wavelet_ip_dual_lagrange = c/normalisation !sqrt(real(ishft (1, scale), PR))

  end function wavelet_ip_dual_lagrange

  !! ATTENTION : il manque le facteur 1/2^(j/2) !!

  function wavelet_ip_dual_lagrange_rat (order, order_dual, fct, degre, scale)
    integer,                       intent(IN)  :: order, order_dual, degre, scale
    type(rationnal), dimension(:), intent(IN)  :: fct
    type(rationnal)                            :: wavelet_ip_dual_lagrange_rat
    type(rationnal)                            :: c
    integer                                    :: n


    c = rat_constr (0_print, 1_print)
    do n=0, degre
       c = rat_add (c , rat_mul (fct(n+1), dual_weight_rat(n+1)))
    end do
    
    wavelet_ip_dual_lagrange_rat = c!/sqrt(real(ishft (1, scale), PR))

  end function wavelet_ip_dual_lagrange_rat

  subroutine wavelet_ipl_test (order, order_dual, degre, j, k)
    integer,                       intent(IN)  :: order, order_dual, degre, j, k
    type(rationnal), dimension(:), allocatable :: phi
    type(rationnal), dimension(:), allocatable :: coef
    real(PR),        dimension(:), allocatable :: fct
    integer                                    :: i, n, scale
    integer                                    :: N1, N2
    character(len=100)                         :: filename

    scale = 0
    do while (ishft (1, scale) < Nx)
       scale = scale + 1
    end do

    N1 = -(order)/2
    N2 = order-order/2

    allocate (phi(ishft (order, scale+j)), fct(ishft (order, scale+j)))
    allocate(coef(order+1))
    
    call wavelet_get_scale_fct_2 (order, j, phi)

    filename="fct_ech"
    call save_tab(phi,filename)

    call wavelet_filter_get_coef (order, coef)
    print *, "Coef de phi"
    call wavelet_display_filter_coef(coef)

    print *, "Fonction : ", size(phi)
    do n=0, degre
       i = 1 + (k+n - N1*ishft(1, j))*Nx
       if (i>size(phi)) then
          fct(n+1) = 0.d0
       else
          fct(n+1) = real(phi(i)%num, PR)/real(phi(i)%den, PR)
       end if
       print *, n+1, i, fct(n+1)
    end do

    z = wavelet_ip_dual_lagrange (order, order_dual, fct, degre, j)

    print *, "Le prod scal vaut : ", z

    deallocate (fct)

  end subroutine wavelet_ipl_test

  !! Calcul d'un produit scalaire

  !! Tests

  subroutine test_wavelet (order, order_dual)
    integer,                       intent(IN)  :: order, order_dual
    integer                                    :: scale, N, i, k
    complex(CPR),    dimension(:), allocatable :: phi, phi_base, f, psi, phi_dual
    real(PR),        dimension(:), allocatable :: coef
    character(len=100)                         :: filename
    real(PR)                                   :: err, N1, N2, N1d, N2d, z
    type(rationnal)                            :: moment_val
    type(rationnal), dimension(:), allocatable :: coef_pol, res
    integer                                    :: degre, indice

    allocate (phi(Nx), psi(Nx), phi_dual(Nx), f(Nx), coef(order+1))

    filename = adjustl(adjustr(main_dir) // adjustl(dir_results)) 
    filename = adjustl(adjustr(filename) // adjustl("ordres.dat"))
    open (unit=99, file=filename, status='replace', action='write')
    write (99, *) order, order_dual
    close (99)

    N1 = -real(ishft(order, -1), PR)
    N2 = real(order-ishft(order, -1), PR)
    N1d = -real(ishft(order_dual, -1), PR) - real(order+order_dual, PR)/2.0d0 + 1
    N2d = real(order_dual-ishft(order_dual, -1), PR) + real(order+order_dual, PR)/2.0d0 - 1


    print *, ""
    print *, "Calcul de la fonction d'echelle avec sa transformee de Fourier"
    scale = 0
    call wavelet_get_scale_fct (order, scale, phi)

    filename="fct_ech"
    call save_tab(phi, filename);

    print *, "Calcul de la fonction d'echelle avec le point fixe"
    scale = 0
    allocate (res(order*Nx))
    call wavelet_get_scale_fct_2 (order, 0, res)

    do i=1, size(phi_dual, 1)
       z = real((i-1)*(size(res, 1)-1), PR)/real(size(phi_dual, 1), PR)
       k = 1+int(z)
       phi_dual(i) = dcmplx(real(res(k)%num, PR)/real(res(k)%den, PR), 0.d0)
    end do
    deallocate (res)

    print *, "Erreur relative entre les 2 methodes : ", rel_err(phi, phi_dual)
    print *, ""

    filename="fct_ech"
    call save_tab(phi_dual, filename);

    !! Test avec la valeur exacte ...

    do i=1, Nx
       phi (i) = wavelet_scale_fct_get_value_r (order, real((i-1)*order, PR)/real(Nx, PR)-real(order/2, PR))
    end do

    print *, "Erreur relative entre les 2 methodes : ", rel_err(phi, phi_dual)
    print *, ""

    filename="fct_ech"
    call save_tab(phi_dual, filename);

    !! Fin test avec la valeur exacte

    do i=0, order
       print *, "Moment ", i, " ", moment(N1, N2, i, phi)
    end do
    print *, ""
    moment_val = wavelet_moment_scale_fct_2(order, order)
    print *, "Moment : ", moment_val%num, "/", moment_val%den

    ! Mother wavelet
    
    print *, "Calcul de l'ondelette mere"
    scale = 0
    call wavelet_get_mother_wavelet (order, order_dual, scale, psi)

    print *, ""
    do i=0, order_dual+4
       print *, "Moment ", i, " ", moment(real(N1-N2d-1, PR)/2.0d0, real(N2-N1d-1, PR)/2.0d0, i, psi)
    end do
    print *, ""

    print *, ""
    do i=0, order_dual+4
       moment_val = wavelet_moment_mother_wavelet_2 (order, order_dual, i)
       print *, "Moment ", i, " ", moment_val%num, "/", moment_val%den
    end do
    print *, ""

    filename="mother_wav"
    call save_tab(psi,filename);

    print *, "Calcul de la fonction d'echelle duale"
    allocate (res((order+2*order_dual-2)*Nx))
    call wavelet_get_dual_scale_fct (order, order_dual, 0, res)

    do i=1, size(phi_dual, 1)
       z = real((i-1)*(size(res, 1)-1), PR)/real(size(phi_dual, 1), PR)
       k = 1+int(z)
       phi_dual(i) = dcmplx(real(res(i)%num, PR)/real(res(i)%den, PR), 0.d0)
    end do

    deallocate (res)

    filename="fct_ech_dual"
    call save_tab(phi_dual,filename);

    scale = 5

    scale = 5
    N = ishft (Nx, scale)
    allocate (phi_base(N))
    call wavelet_get_mother_wavelet (order, order_dual, scale, phi_base)

    filename="fct_ech_rec"
    call save_tab(phi_base,filename);

    
    degre = 3
    allocate (coef_pol(degre+1))
    
    indice = 0
    call get_lagrange_coef (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 0, degre
       call rat_disp (coef_pol(i+1))
    end do
    
    indice = 1
    call get_lagrange_coef (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 0, degre
       call rat_disp (coef_pol(i+1))
    end do
    
    indice = 2
    call get_lagrange_coef (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 0, degre
       call rat_disp (coef_pol(i+1))
    end do
    
    indice = 3
    call get_lagrange_coef (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 0, degre
       call rat_disp (coef_pol(i+1))
    end do

    print *, ""
    print *, " TEST LAGRANGE GENERAL "
    print *, ""

    degre = 4
    deallocate (coef_pol)
    allocate (coef_pol(degre))
    
    indice = -2
    call get_lagrange_coef_general (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 1, degre
       call rat_disp (coef_pol(i))
    end do
    
    indice = -1
    call get_lagrange_coef_general (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 1, degre
       call rat_disp (coef_pol(i))
    end do
    
    indice = 0
    call get_lagrange_coef_general (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 1, degre
       call rat_disp (coef_pol(i))
    end do
    
    indice = 1
    call get_lagrange_coef_general (degre, indice, coef_pol)
    print *, "Indice = ", indice
    do i = 1, degre
       call rat_disp (coef_pol(i))
    end do

    print *, ""
    !call wavelet_ipl_test (order, order_dual, order-2, 1, 3)

    deallocate (phi, phi_dual, phi_base, f, psi, coef, coef_pol)

    print *, ""
    print *, "TEST INTEGRALE SPLINE"
    print *, ""
    print *, "integrale entiere : ", pol_spline_int (order, N1, N2)
    print *, "demi integrale    : ", pol_spline_int (order, N1, 0.5d0*(N1 + N2))
    print *, "demi integrale    : ", pol_spline_int (order, 0.5d0*(N1 + N2), N2)

    do k=1, order
       print *, "Spline    ", pol_spline_real (k, :)
    end do
    print *, ""
    do k=1, order
       print *, "Spline int", pol_spline_real_int (k, :)
    end do

  end subroutine test_wavelet

  !!
  !! Test du produit scalaire
  !!

  function fct_val (x)
    real(PR), intent(IN) :: x
    real(PR)             :: fct_val

    fct_val = exp (-80.0*((x-0.5)*(x-0.5))/2.0d0)
    
  end function fct_val

  subroutine fct_tab(n, tab)
    integer,                intent(IN)  :: n
    real(PR), dimension(:), intent(OUT) :: tab
    integer                             :: i

    do i=1, n
       tab(i) = fct_val(real(i-1, PR)/real(n, PR))
    end do

  end subroutine fct_tab

  subroutine fct_test(order, order_dual, scale, degre)
    integer,               intent(IN)  :: order, order_dual
    integer,               intent(IN)  :: scale, degre
    type(rationnal), dimension(:), allocatable :: phi
    real(PR),dimension(:), allocatable :: fct, ps, tab, tab2
    real(PR)                           :: val
    integer                            :: N1, N2
    integer                            :: k, n, l, dim_tab, dim_ps
    character(len=100)                 :: filename
    integer            :: nb_init, nb_fin, ticker, temps_calcul
    
    N1 = -order/2
    N2 = order - order/2

    call system_clock (COUNT_RATE = ticker)

    dim_tab = Nx*ishft(1, scale)
    allocate (tab(dim_tab), tab2(dim_tab))
    allocate (phi(order*Nx*ishft(1, scale)))
    call wavelet_get_scale_fct_2 (order, scale, phi)
    call fct_tab(dim_tab, tab)

!!$    do n=1, dim_tab
!!$       tab(n) = real(phi(n-N1*dim_tab)%num, PR)/real(phi(n-N1*dim_tab)%den, PR)
!!$    end do

    !! On alloue la fonction 

    allocate (fct(degre+1))

    !! On calcule les produits scalaires qu'il faut.
    !! i.e.   -N2 <= k <= 2^j-N1 (supp(f)=[0,1])

    print *, "Calcul des produits scalaires"
    dim_ps = ishft (1, scale) - N1 + N2 + 1
    allocate (ps(dim_ps))

    call system_clock (COUNT = nb_init)
    do k = -N2, ishft (1, scale) - N1

       do n=0, degre
          fct(n+1) = fct_val(real(k+n, PR)/real(ishft (1, scale), PR))
       end do

       ps (k+N2+1) = wavelet_ip_dual_lagrange (order, order_dual, fct, degre, scale)
    end do
    call system_clock (COUNT = nb_fin)
    temps_calcul = (nb_fin-nb_init)/ticker
    call print_temps_calcul (temps_calcul)

    !! On reconstruit la fonction dans tab2
    print *, "Reconstruction"
    call system_clock (COUNT = nb_init)
    do n = 1, dim_tab

       val = 0.d0
       do k = -N2, ishft (1, scale) - N1
          l = 1-dim_tab*(k+N1)+(n-1)*ishft(1,scale)
          if ((l>=1).and.(l<=size(phi))) then
             val = val + ps (k+N2+1)*real(phi(l)%num, PR)/real(phi(l)%den, PR)
          end if
       end do
       
       tab2(n) = val
       
    end do
    call system_clock (COUNT = nb_fin)
    temps_calcul = (nb_fin-nb_init)/ticker
    call print_temps_calcul (temps_calcul)

    tab2 = tab2*sqrt(real(ishft(1, scale), PR))
    
!!$    call rat_disp (phi(1-N1*Nx*ishft(1, scale)))
!!$    print *, rel_err (tab, tab2)*sqrt(real(ishft(1, scale+3), PR))

    filename="tab1"
    call save_tab(tab,filename);

    filename="tab2"
    call save_tab(tab2,filename);

    print *, "Erreur relative : ", rel_err (tab, tab2)

    deallocate (tab, tab2, ps)

  end subroutine fct_test

  subroutine wavelet_get_ips (order, order_dual, degre, scale, fct, ps)
    integer,                intent(IN)  :: order, order_dual, degre, scale
    real(PR), dimension(:), intent(IN)  :: fct
    real(PR), dimension(:), intent(OUT) :: ps
    real(PR), dimension(:), allocatable :: f
    integer                             :: N1, N2
    integer                             :: i, k, n, j, dim_tab
    
    N1 = -order/2
    N2 = order - order/2

    allocate (f(degre+1))

    do k = -N2, ishft (1, scale) - N1
       do n=0, degre
          j = 1+(k+n)*Nx
          if ((j>=1) .and. (j<=size(fct))) then
             f(n+1) = fct(1+(k+n)*Nx)
          else
             f(n+1) = 0.d0
          end if

       end do
       
       ps (k+N2+1) = wavelet_ip_dual_lagrange (order, order_dual, f, degre, scale)
    end do

    deallocate (f)

  end subroutine wavelet_get_ips

  subroutine wavelet_get_ips_x (order, order_dual, degre, scale, fct, ps)
    integer,                  intent(IN)  :: order, order_dual, degre, scale
    real(PR), dimension(:,:), intent(IN)  :: fct
    real(PR), dimension(:,:), intent(OUT) :: ps
    real(PR), dimension(:),   allocatable :: f
    integer                               :: N1, N2
    integer                               :: i, k, n, j, dim_tab
    
    N1 = -order/2
    N2 = order - order/2

    allocate (f(degre+1))

    ! produit scalaire en x : 

    do i= 1, Nv
       
       do k = -N2, ishft (1, scale) - N1
          do n=0, degre
             j = 1+(k+n)*Nx
             if ((j>=1) .and. (j<=size(fct))) then
                f(n+1) = fct(1+(k+n)*Nx, i)
             else
                f(n+1) = 0.d0
             end if
             
          end do
          
          ps (k+N2+1, i) = wavelet_ip_dual_lagrange (order, order_dual, f, degre, scale)
       end do
    end do

    deallocate (f)

  end subroutine wavelet_get_ips_x

  subroutine wavelet_get_ips_v (order, order_dual, degre, scale, fct, ps)
    integer,                  intent(IN)  :: order, order_dual, degre, scale
    real(PR), dimension(:,:), intent(IN)  :: fct
    real(PR), dimension(:,:), intent(OUT) :: ps
    real(PR), dimension(:),   allocatable :: f
    integer                               :: N1, N2
    integer                               :: i, k, n, j, dim_tab
    
    N1 = -order/2
    N2 = order - order/2

    allocate (f(degre+1))

    ! produit scalaire en x : 

    do i= 1, Nx
       
       do k = -N2, ishft (1, scale) - N1
          do n=0, degre
             j = 1+(k+n)*Nx
             if ((j>=1) .and. (j<=size(fct))) then
                f(n+1) = fct(i, 1+(k+n)*Nx)
             else
                f(n+1) = 0.d0
             end if
             
          end do
          
          ps (i, k+N2+1) = wavelet_ip_dual_lagrange (order, order_dual, f, degre, scale)
       end do
    end do

    deallocate (f)

  end subroutine wavelet_get_ips_v

  subroutine wavelet_get_ips_dual (order, order_dual, degre, scale, fct, ps)
    integer,                intent(IN)  :: order, order_dual, degre, scale
    real(PR), dimension(:), intent(IN)  :: fct
    real(PR), dimension(:), intent(OUT) :: ps
    real(PR), dimension(:), allocatable :: f
    integer                             :: N1d, N2d
    integer                             :: i, k, n, j, dim_tab

    N1d = -(order+order_dual)/2-order_dual/2+1
    N2d = order_dual-order_dual/2+(order+order_dual)/2-1

    allocate (f(degre+1))

    do k = -N2d, ishft (1, scale) - N1d
       do n=0, degre
          j = 1+(k+n)*Nx
          if ((j>=1) .and. (j<=size(fct))) then
             f(n+1) = fct(1+(k+n)*Nx)
          else
             f(n+1) = 0.d0
          end if

       end do
       
       ps (k+N2d+1) = wavelet_ip_dual_lagrange (order, order_dual, f, degre, scale)
    end do

  end subroutine wavelet_get_ips_dual

  subroutine wavelet_rebuild_fct (order, order_dual, scale, phi, ps_in, fct)
    integer,                       intent(IN)  :: order, order_dual, scale
    type(rationnal), dimension(:), intent(IN)  :: phi
    real(PR),        dimension(:), intent(IN)  :: ps_in
    real(PR),        dimension(:), intent(OUT) :: fct
    real(PR)                                   :: val
    integer                                    :: N1, N2
    integer                                    :: n, k, l, dim_tab
    character(len=100)                 :: filename

    N1 = -order/2
    N2 = order - order/2

    dim_tab = Nx*ishft(1, scale)
    do n = 1, dim_tab

       val = 0.d0
       do k = -N2, ishft (1, scale) - N1
          l = 1-dim_tab*(k+N1)+(n-1)*ishft(1,scale)
          if ((l>=1).and.(l<=size(phi))) then
             val = val + ps_in (k+N2+1)*real(phi(l)%num, PR)/real(phi(l)%den, PR)
          end if
       end do
       
       fct(n) = val
       
    end do

    fct = fct*normalisation !sqrt(real(ishft(1, scale), PR))
    
  end subroutine wavelet_rebuild_fct

  !! IMPORTANT
  !! this compute the followiing inner products:
  !! (phi_{fscale, k}(.+v*deltat), phi_dual_{fscale, l})
  !! in tab_ps_advx (v, l, k)
  subroutine wavelet_get_ps_advx ()
    real(PR),        dimension(:), allocatable :: fct
    real(PR)                                   :: v, offset
    integer                                    :: N1, N2, scale
    integer                                    :: k, l, n, j
    integer                                    :: kk, ll

    N1 = -order_bs/2
    N2 = order_bs - order_bs/2

    allocate (fct(degre+1))

    do j = 1, Nv
       !v = real (j-1, PR)/real(Nv, PR)
       offset = -real(ishft (1, fscale-1), PR)*deltat*real (j-1, PR)/real(Nv, PR)
       do l = -N2, ishft (1, fscale) - N1
          !ll = l+N2+1
          do k = -N2, ishft (1, fscale) - N1
             !kk = k+N2+1
             do n=0, degre
                fct (n+1) = wavelet_scale_fct_get_value_r (order_bs, real (l+n-k, PR) + offset)
             end do
             
             tab_ps_advx(j, l+N2+1, k+N2+1) = wavelet_ip_dual_lagrange (order_bs, order_dual_bs, fct, degre, fscale)
          end do
       end do
    end do

    deallocate (fct)

  end subroutine wavelet_get_ps_advx

  !! IMPORTANT
  !! this compute the followiing inner products:
  !! (phi_{fscale, k}(.+v*deltat), phi_dual_{fscale, l})
  !! in tab_ps_advx (v, l, k)
  subroutine wavelet_get_ps_advx_rat ()
    type(rationnal), dimension(:), allocatable :: fct
    type(rationnal)                            :: v, offset, rat_deltat, tmp
    integer                                    :: N1, N2, scale
    integer                                    :: k, l, n, j
    integer                                    :: kk, ll

    N1 = -order_bs/2
    N2 = order_bs - order_bs/2

    rat_deltat = rat_constr (int (1 + int (1.0 - log(deltat)/log(10.0)), print), 1_print)

    allocate (fct(degre+1))

    do j = 1, Nv
       v = rat_constr (int (j-1, print), int (Nv, print))
       offset = rat_constr (int (-ishft (1, fscale), print), 1_print)
       offset = rat_mul (offset, v)
       offset = rat_mul (offset, rat_deltat)
       do l = -N2, ishft (1, fscale) - N1
          ll = l+N2+1
          do k = -N2, ishft (1, fscale) - N1
             kk = k+N2+1
             do n=0, degre
                tmp = rat_constr (int (l+n-k, print), 1_print)
                fct (n+1) = wavelet_scale_fct_get_value_rat (order_bs, rat_add (tmp, offset))
             end do
             
             tmp = wavelet_ip_dual_lagrange_rat (order_bs, order_dual_bs, fct, degre, fscale)
             tab_ps_advx(j, ll, kk) = real (tmp%num, PR)/ (real(tmp%den, PR)*sqrt(real(ishft (1, scale), PR)))
          end do
       end do
    end do

    deallocate (fct)

  end subroutine wavelet_get_ps_advx_rat

  !! IMPORTANT
  !! this compute the followiing inner products:
  !! (phi_{fscale, k}(.+v*deltat), phi_dual_{fscale, l})
  !! in tab_ps_advx (v, l, k)
  subroutine wavelet_get_ps_advx_old (phi)
    type(rationnal), dimension(:), intent(IN)  :: phi
    real(PR),        dimension(:), allocatable :: fct
    real(PR)                                   :: v
    integer                                    :: N1, N2, scale, offset
    integer                                    :: k, l, n, i, j, dim_tab

    N1 = -order_bs/2
    N2 = order_bs - order_bs/2

    dim_tab = Nx*ishft (1, fscale)
    allocate (fct(degre+1))

    do j = 1, Nv
       !v = real (j-1, PR)/real(Nv, PR)
       offset = int (real ((j-1)*Nx*ishft (1, 2*fscale-1), PR)*deltat/real(Nv, PR))
       do l = -N2, ishft (1, fscale) - N1
          
          do k = -N2, ishft (1, fscale) - N1
             fct = 0.d0
             do n=0, degre
                !i = 1 + (l+n-k-N1)*dim_tab - int (real(Nx*ishft (1, 2*fscale-1), PR)*v*deltat)
                i = 1 + (l+n-k-N1)*dim_tab - offset
                if ((i>=1) .and. (i<= dim_tab*order_bs)) then
                   fct (n+1) = real(phi(i)%num, PR)/real(phi(i)%den, PR)
                end if
             end do

             tab_ps_advx(j, l+N2+1, k+N2+1) = wavelet_ip_dual_lagrange (order_bs, order_dual_bs, fct, degre, fscale)
          end do
          
       end do

    end do

    deallocate (fct)

  end subroutine wavelet_get_ps_advx_old

  subroutine wavelet_advection_x (ps_in, v,  ps_out)
    real(PR), dimension(:), intent(IN)  :: ps_in
    integer,                intent(IN)  :: v
    real(PR), dimension(:), intent(OUT) :: ps_out
    real(PR)                            :: val, ps_val
    integer                             :: N1, N2, k, l

    N1 = -order_bs/2
    N2 = order_bs - order_bs/2

    do l = -N2, ishft (1, fscale) - N1
       val = 0.d0
       do k = -N2, ishft (1, fscale) - N1
          ps_val = tab_ps_advx(v, l+N2+1, k+N2+1)
          val = val + ps_val*ps_in(k+N2+1)
       end do
       ps_out (l+N2+1) = val
    end do

    ps_out = ps_out*normalisation !sqrt (real (ishft (1, fscale), PR))

  end subroutine wavelet_advection_x

  subroutine wavelet_advection_v (ps_in, v, E,  ps_out)
    real(PR), dimension(:), intent(IN)  :: ps_in
    integer,                intent(IN)  :: v
    real(PR), dimension(:), intent(IN)  :: E
    real(PR), dimension(:), intent(OUT) :: ps_out
    real(PR)                            :: val, ps_val
    integer                             :: N1, N2, k, l

    N1 = -order_bs/2
    N2 = order_bs - order_bs/2

    !!! A CHANGER !!!

    do l = -N2, ishft (1, fscale) - N1
       val = 0.d0
       do k = -N2, ishft (1, fscale) - N1
          ps_val = tab_ps_advx(v, l+N2+1, k+N2+1)
          val = val + ps_val*ps_in(k+N2+1)
       end do
       ps_out (l+N2+1) = val
    end do

    ps_out = ps_out*normalisation !sqrt (real (ishft (1, fscale), PR))

  end subroutine wavelet_advection_v

  subroutine wavelet_advection_x_old (phi, ps_in, v, deltat, scale, degre, order, order_dual, ps_out)
    type(rationnal), dimension(:), intent(IN) :: phi
    real(PR), dimension(:), intent(IN)  :: ps_in
    real(PR),               intent(IN)  :: v, deltat
    integer,                intent(IN)  :: scale, degre, order, order_dual
    real(PR), dimension(:), intent(OUT) :: ps_out
    real(PR), dimension(:), allocatable :: fct, ps2
    real(PR)                            :: val, ps_val
    integer                             :: N1, N2
    integer                             :: k, l, n, i, dim_tab, dim_ps

    N1 = -order/2
    N2 = order - order/2

    dim_tab = Nx*ishft (1, scale)
    allocate (fct(degre+1))

    do l = -N2, ishft (1, scale) - N1
    
       val = 0.d0
       do k = -N2, ishft (1, scale) - N1
          fct = 0.d0
          do n=0, degre
             i = 1 + (l+n-k-N1)*dim_tab - int (real(Nx*ishft (1, 2*scale-1), PR)*v*deltat)
             if ((i>=1) .and. (i<= dim_tab*order)) then
                fct (n+1) = real(phi(i)%num, PR)/real(phi(i)%den, PR)
             end if
          end do
          
          ps_val = wavelet_ip_dual_lagrange (order, order_dual, fct, degre, scale)
          val = val + ps_val*ps_in(k+N2+1)
          
       end do
       
       ps_out (l+N2+1) = val
    end do
    ps_out = ps_out*normalisation !sqrt (real (ishft (1, scale), PR))

    deallocate (fct)

  end subroutine wavelet_advection_x_old

  subroutine wavelet_write_data (order)
    integer, intent(IN)  :: order
    integer            :: nb_init, nb_fin, ticker, temps_calcul
    character(len=100) :: filename

    call system_clock (COUNT_RATE = ticker)

    print *, " + Construction de Phi x"
    filename="phi_x_128_3"
    call system_clock (COUNT = nb_init)
    call wavelet_get_scale_fct_2 (order, fscale, fct_phi_x)
    call save_tab (fct_phi_x, filename)
    call system_clock (COUNT = nb_fin)
    temps_calcul = (nb_fin-nb_init)/ticker
    call print_temps_calcul (temps_calcul)

!!$    print *, " + Construction de Phi v"
!!$    filename="phi_v_128_3"
!!$    call system_clock (COUNT = nb_init)
!!$    call wavelet_get_scale_fct_2 (order, fscale, fct_phi_v)
!!$    call save_tab (fct_phi_v, filename)
!!$    call system_clock (COUNT = nb_fin)
!!$    temps_calcul = (nb_fin-nb_init)/ticker
!!$    call print_temps_calcul (temps_calcul)

    print *, " + Calculs des produits scalaires"
    filename="ps_128_3_3_8_0.01"
    call system_clock (COUNT = nb_init)
!    call wavelet_get_ps_advx_old (fct_phi_x)
    call wavelet_get_ps_advx ()
!    call wavelet_get_ps_advx_rat ()
!!$    call save_tab (tab_ps_advx, filename)
    call system_clock (COUNT = nb_fin)
    temps_calcul = (nb_fin-nb_init)/ticker
    call print_temps_calcul (temps_calcul)

  end subroutine wavelet_write_data

  subroutine wavelet_read_data (order)
    integer, intent(IN)  :: order
    character(len=5)     :: str
    character(len=100)   :: filename, file
    integer              :: i, j, k
    integer              :: taille_ps
     
    file = "phi_x_128_3"
    filename = adjustl(adjustr(data_dir) // adjustl(file))
    filename = adjustl(adjustr(filename) // adjustl(".dat"))

    open (unit=99, file=filename, action='READ')

    do i=1, order*Nx*ishft (1, fscale)
       read(99,*) fct_phi_x(i)%num, fct_phi_x(i)%den
    end do

    close(99)
     
    file = "phi_v_128_3"
    filename = adjustl(adjustr(data_dir) // adjustl(file))
    filename = adjustl(adjustr(filename) // adjustl(".dat"))

    open (unit=99, file=filename, action='READ')

    do i=1, order*Nv*ishft (1, fscale)
       read(99,*) fct_phi_v(i)%num, fct_phi_v(i)%den
    end do

    close(99)
     
    file = "ps_128_3_3_8_0.01"
    filename = adjustl(adjustr(data_dir) // adjustl(file))
    filename = adjustl(adjustr(filename) // adjustl(".dat"))

    open (unit=99, file=filename, action='READ')

    taille_ps = ishft (1, fscale) + order_bs + 1

    do i=1, Nv
       do j=1, taille_ps
          do k=1, taille_ps
             read(99,*)  tab_ps_advx(i, j, k)
          end do
       end do
    end do

    close(99)
    
  end subroutine wavelet_read_data

  subroutine advx_test(order, order_dual)
    integer,                intent(IN)  :: order, order_dual
    real(PR), dimension(:), allocatable :: fct, fct2, ps_in, ps_out
    type(rationnal)                     :: w
    real(PR)                            :: v, val, x
    integer                             :: N1, N2
    integer                             :: dim_ps
    integer                             :: n, k
    character(len=100)                  :: filename
    integer                             :: t1, d1, t2, d2

    N1 = -order/2
    N2 = order - order/2
    
    print *, ""
    print *, "TEST ADVECTION EN x"
    print *, ""

    print *, " + Calcul des moments"
    call wavelet_moments_dual_scale_fct (order, order_dual, degre, moments)

    print *, " + Calcul des coef de Lagrange"
    do n=0, degre
       call get_lagrange_coef (degre, n, lagrange(n+1,:))
    end do

    print *, " + Calcul des poids"
    do n=0, degre
       w = rat_constr (0_print, 1_print)
       do k = 1, degre+1
          w = rat_add (w, rat_mul (lagrange(n+1, k), moments(k)))
       end do
       dual_weight_rat(n+1) = w
       dual_weight(n+1) = real(w%num, PR)/real(w%den, PR)
    end do

    call wavelet_write_data (order)

!!$    print *, " + Construction de Phi"
!!$    call system_clock (t1, d1)
!!$    allocate (phi(order*Nx*ishft(1, fscale)))
!!$    call wavelet_get_scale_fct_2 (order, fscale, phi)
!!$    call system_clock (t2, d2)
!!$    call get_time (t1, t2, d1)
!!$
!!$    print *, " + Calculs des produits scalaires"
!!$    call system_clock (t1, d1)
!!$    filename="ps_128_3_3_0.01"
!!$    call wavelet_get_ps_advx (phi, degre)
!!$    call save_tab_real_3d (tab_ps_advx, filename)
!!$    !call read_ps (filename)
!!$    call system_clock (t2, d2)
!!$    call get_time (t1, t2, d1)
    
    dim_ps = ishft (1, fscale) - N1 + N2 + 1
    allocate (ps_in(dim_ps), ps_out(dim_ps))

    !! creation de la fonction initiale

    print *, " + Construction de la fonction initiale"
    allocate (fct(Nx*ishft(1, fscale)), fct2(Nx*ishft(1, fscale)))
    !call fct_test_advx (fct)
    call fct_test_advx2 (fct)

    filename="tab1"
    call save_tab(fct,filename);
    
    !! calcul des ps
    call system_clock (t1, d1)
    call wavelet_get_ips (order, order_dual, degre, fscale, fct, ps_in)
    call system_clock (t2, d2)
    call get_time (t1, t2, d1)

    n = 1
    v = 0.5d0
    Nt=20
    print *, " + Debut de la boucle"
    call system_clock (t1, d1)
    do while (n<=Nt)
       call wavelet_advection_x (ps_in, Nv/2,  ps_out)
       ps_in = ps_out
       n = n+1
       print *, "    n = ", n
    end do
    call system_clock (t2, d2)
    print *, "temps : ", t1, d1, t2, d2
    call get_time (t1, t2, d1)

    !! reconstruction
    print *, " + Reconstruction de la fonction"
    call wavelet_rebuild_fct (order, order_dual, fscale, fct_phi_x, ps_in, fct2)

    filename="tab2"
    call save_tab(fct2,filename);

    call fct_test_advx3 (0.55d0, fct)

    print *, rel_err (fct, fct2)*sqrt(real(ishft(1, fscale+3), PR))

    deallocate (fct, fct2, ps_in, ps_out)

  end subroutine advx_test

  !! Debut du vrai code

  subroutine w_advx_demi_pas_temps (degre, f, g)
    integer,                  intent(IN)  :: degre
    real(PR), dimension(:,:), intent(IN)  :: f
    real(PR), dimension(:,:), intent(OUT) :: g
    integer                               :: j
    
    do j = 1, Nv
       call wavelet_get_ips (order_bs, order_dual_bs, degre, fscale, f(:,j), ps_provi1)
       call wavelet_advection_x (ps_provi1, j,  ps_provi2)
       call wavelet_rebuild_fct (order_bs, order_dual_bs, fscale, fct_phi_x, ps_provi2, g(:,j))
    end do

  end subroutine w_advx_demi_pas_temps

  subroutine wavelet_vp_1d (order, order_dual)
    integer, intent(IN)  :: order, order_dual
    real(PR), dimension(:,:), allocatable :: fct, fct2
    real(PR), dimension(:,:), allocatable :: ps_in_x, ps_out_x
    real(PR), dimension(:,:), allocatable :: ps_in_v, ps_out_v
    type(rationnal)      :: w
    integer              :: test_save
    real(PR)             :: temps
    character(len=100)   :: filename
    character(len=4)     :: num_str
    integer              :: nb_init, nb_fin, ticker
    integer              :: N1, N2
    integer              :: dim_ps
    integer              :: n, k, temps_calcul
    integer                             :: t1, d1, t2, d2

    ! initialisation de quelques varaibles
    call system_clock (COUNT_RATE = ticker)

    N1 = -order/2
    N2 = order - order/2
    dim_ps = ishft (1, fscale) - N1 + N2 + 1

    test_save = Nt/save_nbr
    filename = "save"

    n     = 0
    temps = 0.d0
    
    ! Allocations
    allocate (ps_in_x(dim_ps, Nv), ps_out_x(dim_ps, Nv))
    allocate (ps_in_v(Nx, dim_ps), ps_out_v(Nx, dim_ps))
    
    print *, ""
    print *, "TEST ADVECTION"
    print *, ""

    print *, " + Calcul des moments"
    call wavelet_moments_dual_scale_fct (order, order_dual, degre, moments)

    print *, " + Calcul des coef de Lagrange"
    do n=0, degre
       call get_lagrange_coef (degre, n, lagrange(n+1,:))
    end do

    print *, " + Calcul des poids"
    do n=0, degre
       w = rat_constr (0_print, 1_print)
       do k = 1, degre+1
          w = rat_add (w, rat_mul (lagrange(n+1, k), moments(k)))
       end do
       dual_weight_rat(n+1) = w
       dual_weight(n+1) = real(w%num, PR)/real(w%den, PR)
    end do

    !!! Sauvegardes preliminaires des donnees

    call wavelet_write_data (order)

    !!! Chargement preliminaire des donnees

    !call wavelet_read_data (degre)
  
    print *, " + Construction de la fonction initiale"
    call fct_test_landau(Kx, 0.005d0, f_initial)
    !call poisson_1d (f_initial, E_initial)
    
    !! preparation of f

  end subroutine wavelet_vp_1d

end module wavelet
