module py2f90
    implicit none
    contains
        subroutine transferMat2For(matrix,n1,n2,matrix2)
        implicit none
        integer,intent(in)::n1,n2
        real(8),intent(inout)::matrix(n1,n2)
        real(8),intent(out)::matrix2(n1,n2)
        integer::i,j
        ! loop
        do i=1,n1
            do j=1,n2
                matrix(i,j) = matrix(i,j) + 1.0 
                write(*,"('Row:',i4,1x,'Col:',i4,1x,'Value:',1x,F5.2)")i,j,matrix(i,j)
            enddo
        enddo

        matrix2 = matrix
        end subroutine
end module

Recursive function phi(n,y) result(out)
    implicit none 
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: y
    real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
    integer :: n
    complex(kind = DP) :: out

    if(n<0) then 
        out = 1.0_DP
    else if(n ==0) then
        out = PI**(-0.25_DP)*exp(-0.5_DP*y**2)
    else if(n ==1) then
        out = PI**(-0.25_DP)*sqrt(2.0_DP)*y*exp(-0.5_DP*y**2)
    else if(n>=2) then
        out = sqrt(2.0_DP/(n))*y*phi(n-1,y)-sqrt((n-1.0_DP)/n)*phi(n-2,y)
    end if
end function


subroutine phi_collo_array2(y_collo,n_order,phi_ij)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer, intent(in) :: n_order
    real(kind = DP), intent(in) :: y_collo(n_order)
    complex(kind = DP),intent(out) :: phi_ij(n_order,n_order)
    integer :: i,j
    real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)

    do i = 1,n_order
        ! print *,i
        do j = 1,n_order
            if(j-1<0) then
                return
            else if(j-1 ==0) then
                phi_ij(i,j) = PI**(-0.25_DP)*exp(-0.5_DP*y_collo(i)**2)
            else if(j-1 ==1) then
                phi_ij(i,j) = PI**(-0.25_DP)*sqrt(2.0_DP)*y_collo(i)*exp(-0.5_DP*y_collo(i)**2)
            else if(j-1>=2) then
                phi_ij(i,j) = sqrt(2.0_DP/(j-1))*y_collo(i)*phi_ij(i,j-1)-sqrt((j-1-1.0_DP)/(j-1))*phi_ij(i,j-2)
            end if
        end do
    end do
end subroutine

subroutine phi_collo_array(y_collo,n_order,phi_ij)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer, intent(in) :: n_order
    real(kind = DP), intent(in) :: y_collo(n_order)
    complex(kind = DP),intent(out) :: phi_ij(n_order,n_order)
    integer :: i,j
    real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
    integer(kind = 8) :: tmp(n_order-2)

    do i = 1,n_order
        ! print *,i
        do j = 1,n_order
            if(j-1<0) then
                return
            else if(j-1 ==0) then
                phi_ij(i,j) = PI**(-0.25_DP)
            else if(j-1 ==1) then
                phi_ij(i,j) = PI**(-0.25_DP)*sqrt(2.0_DP)*y_collo(i)
            else if(j-1 ==2) then
                phi_ij(i,j) = sqrt(2.0_DP/(j-1))*y_collo(i)*phi_ij(i,j-1)-sqrt((j-1-1.0_DP)/(j-1))*phi_ij(i,j-2)
                tmp(j-2) = int8(log(phi_ij(i,j)))
                ! print *,tmp(j-2)
                phi_ij(i,j) = exp(-tmp(j-2)+ .0_DP)*phi_ij(i,j)
            else if(j-1>=3) then
                phi_ij(i,j) = sqrt(2.0_DP/(j-1))*y_collo(i)*phi_ij(i,j-1)-sqrt((j-1-1.0_DP)/(j-1))*phi_ij(i,j-2) &
                * exp(-tmp(j-3)+.0_DP)
                tmp(j-2) = int8(log(phi_ij(i,j)))
                ! print *,tmp(j-2)
                phi_ij(i,j) = exp(-tmp(j-2)+ .0_DP)*phi_ij(i,j)
            end if
        end do
        do j = 1,n_order
            if(j-1<0) then
                return
            else if(j-1 <2) then
                phi_ij(i,j) = phi_ij(i,j) *exp(-0.5_DP*y_collo(i)**2)
            else if(j-1>=2) then
                phi_ij(i,j) = phi_ij(i,j) *exp(-0.5_DP*y_collo(i)**2+sum(tmp(1:j-2)))
            end if
        end do
    end do
end subroutine

subroutine phi_grid_array(y,ny,n_order,phi_ij)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer, intent(in) :: n_order,ny
    real(kind = DP), intent(in) :: y(ny)
    complex(kind = DP),intent(out) :: phi_ij(ny,n_order)
    real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
    integer :: i,j
    integer(kind = 8) :: tmp(n_order-2)
    do i = 1,ny
        ! print *,i
        do j = 1,n_order
            if(j-1<0) then
                return
            else if(j-1 ==0) then
                phi_ij(i,j) = PI**(-0.25_DP)
            else if(j-1 ==1) then
                phi_ij(i,j) = PI**(-0.25_DP)*sqrt(2.0_DP)*y(i)
            else if(j-1 ==2) then
                phi_ij(i,j) = sqrt(2.0_DP/(j-1))*y(i)*phi_ij(i,j-1)-sqrt((j-1-1.0_DP)/(j-1))*phi_ij(i,j-2)
                tmp(j-2) = int8(log(phi_ij(i,j)))
                ! print *,tmp(j-2)
                phi_ij(i,j) = exp(-tmp(j-2)+ .0_DP)*phi_ij(i,j)
            else if(j-1>=3) then
                phi_ij(i,j) = sqrt(2.0_DP/(j-1))*y(i)*phi_ij(i,j-1)-sqrt((j-1-1.0_DP)/(j-1))*phi_ij(i,j-2) &
                * exp(-tmp(j-3)+.0_DP)
                tmp(j-2) = int8(log(phi_ij(i,j)))
                ! print *,tmp(j-2)
                phi_ij(i,j) = exp(-tmp(j-2)+ .0_DP)*phi_ij(i,j)
            end if
        end do
        do j = 1,n_order
            if(j-1<0) then
                return
            else if(j-1 <2) then
                phi_ij(i,j) = phi_ij(i,j) *exp(-0.5_DP*y(i)**2)
            else if(j-1>=2) then
                phi_ij(i,j) = phi_ij(i,j) *exp(-0.5_DP*y(i)**2+sum(tmp(1:j-2)))
            end if
        end do
    end do
end subroutine

subroutine phi_grid_array2(y,ny,n_order,phi_ij)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer, intent(in) :: n_order,ny
    real(kind = DP), intent(in) :: y(ny)
    complex(kind = DP),intent(out) :: phi_ij(ny,n_order)
    real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
    integer :: i,j

    do i = 1,ny
        ! print *,i
        do j = 1,n_order
            if(j-1<0) then
                return
            else if(j-1 ==0) then
                phi_ij(i,j) = PI**(-0.25_DP)*exp(-0.5_DP*y(i)**2)
            else if(j-1 ==1) then
                phi_ij(i,j) = PI**(-0.25_DP)*sqrt(2.0_DP)*y(i)*exp(-0.5_DP*y(i)**2)
            else if(j-1>=2) then
                phi_ij(i,j) = sqrt(2.0_DP/(j-1))*y(i)*phi_ij(i,j-1)-sqrt((j-1-1.0_DP)/(j-1))*phi_ij(i,j-2)
            end if
        end do
    end do
end subroutine


! function ddy(y,p,o)
!     implicit none
!     interface
!         Recursive function phi(n,y) result(out)
!             implicit none
!             integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
!             real(kind = DP) :: y
!             real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
!             integer :: n
!             complex(kind = DP) :: out
!         end function
!     end interface
!     integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
!     real(kind = DP) :: y,p
!     integer :: o
!     real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
!     complex(kind = DP) :: ddy
!     if(o ==0) then
!         ddy = -y*PI**(-0.25_DP)*exp(-0.5_DP*y**2)
!     else if(o ==1) then
!         ddy = sqrt(2.0_DP)*PI**(-0.25_DP)*exp(-0.5_DP*y**2)*(1.0_DP-y**2)
!     else if(o>=2) then
!         ddy = -y*phi(o,y) + sqrt(2.0_DP*o)*phi(o-1,y)
!     end if
! end function

! function dddy(y,p,o)
!     implicit none
!     interface
!         Recursive function phi(n,y) result(out)
!             implicit none
!             integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
!             real(kind = DP) :: y
!             real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
!             integer :: n
!             complex(kind = DP) :: out
!         end function
!     end interface
!     integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
!     real(kind = DP) :: y,p
!     integer :: o
!     complex(kind = DP) :: dddy
!     dddy = (-(2*o+1)+y**2/p**2) *phi(o,y/p)
! end function

subroutine phi_matrix_collo(y_collo,n_order,p,phi_ij)
    implicit none
    interface
        Recursive function phi(n,y) result(out)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y
            real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
            integer :: n
            complex(kind = DP) :: out
        end function
    end interface
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer, intent(in) :: n_order
    real(kind = DP), intent(in) :: y_collo(n_order),p
    complex(kind = DP),intent(out) :: phi_ij(n_order,n_order)
    integer :: i,j
    do i = 1,n_order
        ! print *,i
        do j = 1,n_order
            phi_ij(i,j) = phi(j-1,y_collo(i)/p)
        end do
    end do
end subroutine

subroutine phi_matrix_grid(y,ny,n_order,p,phi_ik)
    implicit none
    interface
        Recursive function phi(n,y) result(out)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y
            real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
            integer :: n
            complex(kind = DP) :: out
        end function
    end interface
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer, intent(in) :: n_order,ny
    real(kind = DP), intent(in) :: y(ny),p
    complex(kind = DP),intent(out) :: phi_ik(ny,n_order)
    integer :: i,k
    do i = 1,ny
        ! print *,i
        do k = 1,n_order
            phi_ik(i,k) = phi(k-1,y(i)/p)
        end do
    end do
end subroutine


program main
    implicit none
    interface
        Recursive function phi(n,y) result(out)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y
            real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
            integer :: n
            complex(kind = DP) :: out
        end function

        function ddy(y,p,o)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y,p
            integer :: o
            real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
            complex(kind = DP) :: ddy
        end function
    end interface
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real :: PI=4*DATAN(1.D0)
    real(kind = 8) :: out(5,4)
    integer :: i,j
    complex(kind = DP) :: test
    ! do j = 1,40
    !     do i = 1,40
    !         out(j,i) = j-i
    !     end do
    ! end do
    test = phi(1000,1.0_DP)
    print *, test
end program