      subroutine hetrfm2_complex(psi_in,psi_out,iop)
!******************************************
!     2-D transform in Serie of Hermite FUNCTIONS !
! Input: a,grw: output from heinit.f
!     iop=1: spectral to physics
!     iop=0: physics to spectral
!      x: sequence to be transformed
! Output: y: transformed sequence
!*****************************************      
	use constant, only : nx,ny,ax,ay,gxw,gyw
	implicit none
	integer(4) :: iop
	complex(8),intent(in) :: psi_in(0:nx,0:ny)
	complex(8),intent(out) :: psi_out(0:nx,0:ny)

	integer(4) :: i,j
	complex(8) :: psi_ast(0:nx,0:ny)

	do j = 0,ny
		call hetrfm1_complex(nx+1,psi_in(:,j),psi_ast(:,j),ax,gxw,iop)
	enddo

	do i = 0,nx
		call hetrfm1_complex(ny+1,psi_ast(i,:),psi_out(i,:),ay,gyw,iop)
	enddo

	end subroutine hetrfm2_complex

	subroutine hetrfm1_complex(n,x,y,a,grw,iop)
!******************************************
!     1-D transform in Serie of Hermite FUNCTIONS !
! Input: a,grw: output from heinit.f
!     iop=1: spectral to physics
!     iop=0: physics to spectral
!      x: sequence to be transformed
! Output: y: transformed sequence
!
	implicit none
	integer(4),intent(in) :: n,iop
	complex(8),intent(in) :: x(n)
	complex(8),intent(out) :: y(n)
	real(8),intent(in) :: a(n,0:n-1),grw(n)

	integer(4) :: j
	real(8) :: xr(n),xi(n)
	real(8) :: xr1(n),xi1(n)

	do j = 1,n
		xr(j) = dreal(x(j))
		xi(j) = dimag(x(j))
	enddo

	call hetrfm1(n,xr,xr1,a,grw,iop)
	call hetrfm1(n,xi,xi1,a,grw,iop)

	do j = 1,n
		y(j) = dcmplx(xr1(j),xi1(j))
	enddo

	end subroutine hetrfm1_complex

      subroutine heinit(n,x,a,grw)
!******************************************
!* Initialization for Scaled Hemite functions U_i(x)=alpha_i H_i e^(-x*x/2)
!* such that \int U_i u_j dx=delta_ij
!* 
!* Input: x(j), j=1,n: Hermit-Gauss points
!* Output: a(j,i)=U_i(x_j):=alpha_i H_i(x(j))*exp(-x(j)*x(j)/2)
!*         grw(i)= scaled Hermite-Gauss weight
!*
!*  It works for n up to at least 1400 ! It will eventually overflow for n very
!*  large...
!*  To make it work for even larger n, one should do it as described in
!*  my SIAM laguerre paper, see ~/laguerre/lib/lainit.f !
!*****************************************
      implicit double precision (a-h,o-z)
      dimension x(n),a(n,0:n-1),grw(n)
      spi=sqrt(acos(-1.d0))

      do j=1,n
         t0=1.d0/sqrt(spi)
         a(j,0)=t0*exp(-.5d0*x(j)*x(j)/2)
         a(j,1)=2*x(j)*a(j,0)/sqrt(2.d0)
         
         do i=2,n-1
             a(j,i)=2*x(j)*a(j,i-1)/sqrt(2.d0*i)-2*(i-1)*a(j,i-2)/sqrt(4.d0*i*(i-1))
          enddo
          do i=0,n-1
             a(j,i)=a(j,i)*exp(-.5d0*x(j)*x(j)/2)
          enddo
      enddo

!**** Compute the scaled Gauss-Radau weight

      do j=1,n
         grw(j)=1.d0/(n*a(j,n-1)**2)
      enddo
      return
      end

      subroutine hetrfm1(n,x,y,a,grw,iop)
!******************************************
!*     1-D transform in Serie of Hermite FUNCTIONS !
!* Input: a,grw: output from heinit.f
!*     iop=1: spectral to physics
!*     iop=0: physics to spectral
!*      x: sequence to be transformed
!* Output: y: transformed sequence
!*****************************************      
	implicit none
	integer(4) :: n,iop
      real(8):: x(n),y(n),a(n,0:n-1),grw(n)

	real(8) :: xt(n)
	integer(4) :: j

      if (iop.eq.0) then
         do j=1,n
            xt(j)=x(j)*grw(j)
         enddo
         call dgemv('t',n,n,1.d0,a,n,xt,1,0.d0,y,1)
      else
         call dgemv('n',n,n,1.d0,a,n,x,1,0.d0,y,1)
      endif

      return
      end
