PROGRAM hof c harmonic oscillator wavefunctions IMPLICIT none integer*4 n,nw,np,ii,jj,k,kk,i,j real*8 pi,akfac,x,dx,fik,xmax,alpha real*8,allocatable::p(:,:),ank(:,:),w(:),x0(:),am0(:) pi=4.0d0*atan(1.0d0) write(6,*)'n, np:' read(5,*)n,np nw=1 allocate(p(n,n),ank(nw,n),w(nw),x0(nw),am0(nw)) am0(1)=1.0d0 w(1)=1.0d0 x0(1)=0.0d0 write(6,*)'Initializing HO basis, kmax = ',n do 103 ii=1,n do 103 jj=1,n 103 p(jj,ii)=0.0d0 p(1,1)=1.0d0 do 105 ii=1,nw 105 ank(ii,1)=dsqrt(dsqrt(am0(ii)*w(ii)/pi)) k=1 102 k=k+1 akfac=1.0d0 do 101 kk=1,k-1 101 akfac=akfac*dble(kk) do 1 ii=1,nw 1 ank(ii,k)=dsqrt(dsqrt(am0(ii)*w(ii)/pi)/2.0d0**(k-1)/akfac) p(k,1 )=-p(k-1,2) do 104 j=2,k-2 104 p(k,j )=(2.0d0*p(k-1,j-1)-dble(j)*p(k-1,j+1)) if(k-2.gt.0)p(k,k-1)=2.0d0*p(k-1,k-2) p(k,k )=2.0d0*p(k-1,k-1) write(6,*)k,n if(k.lt.n)goto 102 write(6,*)'initiated' c alpha=am0(1)*w(1) xmax=3.0*dsqrt(dble(n)/alpha) dx=2.0d0*xmax/dble(np) x=-xmax-dx open(9,file='t.prn') do 2 i=1,np x=x+dx 2 write(9,900)x,fik(1,0,n,x,n,1,p,ank,w,x0,am0) 900 format(g12.4,' ',g12.4) close(9) end c ------------------------------------------------------------ function fik(ivar,id,k,x,n,nvar0,p,ank,w,x0,am0) implicit none integer nvar0,ivar,n real*8 fik,x,hk,hkp,hkpp,y, 1ff,dw,p(n,n),ank(nvar0,n),w(nvar0),x0(nvar0),am0(nvar0) integer*4 k,id,kk dw=dsqrt(w(ivar)*am0(ivar)) y=dw*(x-x0(ivar)) ff=ank(ivar,k)*exp(-y**2/2.0d0) hk=0.0d0 hkp=0.0d0 hkpp=0.0d0 do 1 kk=1,k 1 hk=hk+p(k,kk)*y**(kk-1) if(id.gt.0)then do 2 kk=2,k 2 hkp=hkp+dble(kk-1)*p(k,kk)*y**(kk-2) if(id.gt.1)then do 3 kk=3,k 3 hkpp=hkpp+dble((kk-1)*(kk-2))*p(k,kk)*y**(kk-3) endif endif if(id.eq.0)then fik=ff*hk else if(id.eq.1)then fik=ff*(hkp-y*hk)*dw else if(id.eq.2)then fik=ff*(hk*(y**2-1.0d0)-2.0d0*y*hkp+hkpp)*dw**2 else write(6,*)'higher ders not implemented' stop endif endif endif return end