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


function ddy(y,o,phi_ij,phi_ij01)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: y
    integer :: o
    real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
    complex(kind = DP) :: ddy, phi_ij,phi_ij01
    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_ij + sqrt(2.0_DP*o)*phi_ij01
    end if
end function

function dddy(y,o,phi_ij)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: y
    integer :: o
    complex(kind = DP) :: dddy,phi_ij
    dddy = (-(2*o+1)+y**2) *phi_ij
end function

function Bx(y,V_A,Lm,Bx0)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: y,Bx,V_A
    real(kind = DP),intent(in) :: Lm,Bx0
    Bx = y/Lm* V_A* exp(-y**2/(2*Lm**2)) * exp(0.5_DP) +Bx0
end function

function dBxdy(y,V_A,Lm)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: y,dBxdy,V_A
    real(kind = DP),intent(in) :: Lm
    dBxdy = exp(0.5_DP) * V_A / Lm * exp(-y**2/(2*Lm**2)) * (1.0_DP -y**2/(Lm**2))
end function

function ddBxdy(y,V_A,Lm)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: y,ddBxdy,V_A
    real(kind = DP),intent(in) :: Lm
    ddBxdy = exp(0.5_DP) * V_A *y / Lm**3 * exp(-y**2/(2*Lm**2)) * (y**2/(Lm**2) -3.0_DP)
end function


subroutine assign(n_order,n_p,n_eqt,n_var,p,y_collo,Bx_mean,By0,B_0,V_A,mu,k_x,Lm,tao_r,tao_d,H,f,phi_const)
    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,o,phi_ij,phi_ij01)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y
            integer :: o
            real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
            complex(kind = DP) :: ddy, phi_ij,phi_ij01
        end function

        function dddy(y,o,phi_ij)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y
            integer :: o
            complex(kind = DP) :: dddy,phi_ij
        end function
        
        function Bx(y,V_A,Lm,Bx0)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y,Bx,V_A
            real(kind = DP),intent(in) :: Lm,Bx0
        end function

        function dBxdy(y,V_A,Lm)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y,dBxdy,V_A
            real(kind = DP),intent(in) :: Lm
        end function

        function ddBxdy(y,V_A,Lm)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: y,ddBxdy,V_A
            real(kind = DP),intent(in) :: Lm
        end function
    end interface
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP),parameter :: PI=4.D0*DATAN(1.D0)
    integer(kind = 4), intent(in) ::n_order,n_p,n_eqt,n_var
    complex(kind = DP), intent(in) :: phi_const(n_order,n_order)
    real(kind = DP), intent(in) :: p,By0,B_0,V_A,mu,k_x,Lm,tao_r,tao_d,y_collo(n_p),Bx_mean

    ! OUTPUT varibles: 
    complex(kind = DP), intent(out) :: H(n_p*n_eqt,n_order*n_var),f(n_eqt*n_p)

    integer :: i,j
    real(kind = DP) :: Bx0
    complex(kind = DP):: phi_ij,ddy_value,j1=(0.0_DP,1.0_DP),dBx0dy

    do i = 1,n_p
        ! print *, i
        do j = 1,n_order
            phi_ij = phi_const(i,j)
            Bx0 = Bx(y_collo(i),V_A,Lm,Bx_mean)
            dBx0dy = dBxdy(y_collo(i),V_A,Lm)
            if(j-1>=2) then
                ddy_value = ddy(y_collo(i),j-1,phi_ij,phi_const(i,j-1))
            else 
                ddy_value = ddy(y_collo(i),j-1,phi_ij,(1.0_DP,0.0_DP))
            end if
            H(i+n_p*0,j+n_order*0) = 1/tao_d*phi_ij
            H(i+n_p*1,j+n_order*0) = y_collo(i) * phi_ij
            H(i+n_p*2,j+n_order*0) = j1 * k_x * phi_ij
            H(i+n_p*3,j+n_order*0) = (B_0 + Bx0*j1*k_x)* phi_ij
            H(i+n_p*4,j+n_order*0) = 0
            
            H(i+n_p*0,j+n_order*1) = -H(i+n_p*1,j+n_order*0)
            H(i+n_p*1,j+n_order*1) = 1/tao_d*phi_ij
            H(i+n_p*2,j+n_order*1) = ddy_value
            H(i+n_p*3,j+n_order*1) = -(dBx0dy * phi_ij)
            H(i+n_p*4,j+n_order*1) = (j1*k_x*Bx0+B_0)* phi_ij

            H(i+n_p*0,j+n_order*2) = j1*k_x*phi_ij
            H(i+n_p*1,j+n_order*2) = ddy_value
            H(i+n_p*2,j+n_order*2) = 1/tao_r*phi_ij

            H(i+n_p*0,j+n_order*3) = (B_0 -j1*Bx0*k_x)*phi_ij
            H(i+n_p*1,j+n_order*3) = 0
            H(i+n_p*3,j+n_order*3) = mu*(-k_x**2*phi_ij+dddy(y_collo(i),j-1,phi_ij))

            H(i+n_p*0,j+n_order*4) = -dBx0dy * phi_ij
            H(i+n_p*1,j+n_order*4) = (-j1*k_x*Bx0+B_0)*phi_ij
            H(i+n_p*4,j+n_order*4) = H(i+n_p*3,j+n_order*3)
        end do
    end do
   
    do i = 1,n_p
        ! f(2*n_p+i) = phi_const(i,1) / tao_r
        f(2*n_p+i) = exp(-0.5_DP*(y_collo(i))**2) / tao_r
        ! if(abs(y_collo(i)) <= 2.544496_DP) then
        !     f(2*n_p+i) = cos(y_collo(i)/1.6198765_DP) / tao_r
        ! end if
        ! f(3*n_p+i) = -mu*ddBxdy(y_collo(i),V_A,Lm)
    end do
end subroutine

subroutine assign_3d(n_var,k_x,n_order,ny,nx,u,v,eta,bx,by,a,x,phi_const)
    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(kind = 4), intent(in) ::n_order,n_var,ny,nx
    real(kind = 8), intent(in) :: k_x,x(nx)
    complex(kind = 8), intent(in) :: a(n_order*n_var),phi_const(ny,n_order)
    ! OUTPUT varibles: 
    complex(kind = 8), intent(out) :: u(ny,nx),v(ny,nx),eta(ny,nx),bx(ny,nx),by(ny,nx)

    integer :: i,j,k
    complex(kind = 8) :: phi_ik

    do i =1,ny
        ! print *, i
        do j = 1,nx
            do k = 1,n_order
                phi_ik = phi_const(i,k)
                u(i,j) =  phi_ik*a(k) + u(i,j)
                v(i,j) =  phi_ik*a(n_order+k) + v(i,j)
                eta(i,j) =  phi_ik*a(n_order*2+k) +eta(i,j)
                bx(i,j) =  phi_ik*a(n_order*3+k) +bx(i,j)
                by(i,j) =  phi_ik*a(n_order*4+k) +by(i,j)
            end do
            u(i,j) =exp((0,1)*k_x*x(j)) * u(i,j)
            v(i,j) =exp((0,1)*k_x*x(j)) * v(i,j)
            eta(i,j) =exp((0,1)*k_x*x(j)) * eta(i,j)
            bx(i,j) =exp((0,1)*k_x*x(j)) * bx(i,j)
            by(i,j) = exp((0,1)*k_x*x(j)) *by(i,j)
        end do
    end do
end subroutine