      program boxm
c Box-Mueller algoritme for aa generere Gaussisk-fordelte tall
      parameter(nmax=250000,nh=100,nm=-nh,hmax=8.,nc=500)
      dimension nhist(nm:nh),corr(nc),xold(nc)
      hinv=1./hmax
      nbin=nh-nm+1
      nc2=nc/2
      do i=nm,nh
      nhist(i)=0
      enddo
      rinv=0.5/(2.**31-1.)
      ibm=1955
      do i=1,1000
      ibm=ibm*16807
      enddo
      pi2=8*atan(1.)
      do i=1,nmax
      ibm=ibm*16807
      y1=rinv*float(ibm)+0.5
      ibm=ibm*16807
      y2=rinv*float(ibm)+0.5
      x1=sqrt(-2.*log(y1))*cos(pi2*y2)
      x2=sqrt(-2.*log(y1))*sin(pi2*y2)
      nh1=(x1*hinv)*nbin
      nhist(nh1)=nhist(nh1)+1
      nh2=(x2*hinv)*nbin
      nhist(nh2)=nhist(nh2)+1
      do j=2,nc
      xold(j-1)=xold(j)
      enddo
      xold(nc)=x1
      if(i.gt.nc2) then
      do j=1,nc
      corr(j)=corr(j)+xold(j)*xold(nc)
      enddo
      endif
      do j=2,nc
      xold(j-1)=xold(j)
      enddo
      xold(nc)=x2
      if(i.gt.nc2) then
      do j=1,nc
      corr(j)=corr(j)+xold(j)*xold(nc)
      enddo
      endif
      enddo
      nhist(0)=nhist(0)/2
      open(unit=1,file='boxm1.dat',status='unknown')
      open(unit=2,file='boxm2.dat',status='unknown')
      do i=nm,nh
      write(1,*) float(i)*hmax/nbin,nhist(i)
      enddo
      close(1)
      do j=nc,1,-1
      write(2,*) nc-j,corr(j)/corr(nc)
      enddo
      close(2)
      end

