      program prog10
C
C     ２次元の粒子衝突を Monte-Carlo 法でシミュレーションするプログラム
C     function rnd は前回のものを使って(コピー・ペースト)下さい。 
C
      implicit real*8(a-h,o-z)
      integer*4 iseed
      real*8 hist(-50:50),p(2,10000)
C
C 準備 ----------------------------------------------------------------*
C パラメータ設定
      iseed = 123456789
      pi = 4*atan(1.0d0)
      ncol = 100000
      n = 10000
      dx = 0.02d0
C ファイルを開く
      open(16,file='prog10.dat')
C
C ヒストグラムの初期化
      do j = -50,50
         hist(j) = 0.0d0
      end do
C 運動量の初期値設定
      px = 0.0d0
      py = 0.0d0
      do i = 1,n
         p(1,i) = rnd(iseed)
         p(2,i) = rnd(iseed)
         px = px + p(1,i)
         py = py + p(2,i)
      end do
      px = px / n
      py = py / n
C 運動量の平均を 0 にする。
      do i = 1, n
         p(1,i) = p(1,i) - px
         p(2,i) = p(2,i) - py
      end do
C
C 実際の衝突計算: ２次元の弾性散乱を ncol 回繰り返す。-----------------*
C
      do k = 1, ncol
C 衝突する粒子対をランダムに選ぶ。
         i1 =int(rnd(iseed)*n)+1
         i2 =int(rnd(iseed)*n)+1
C
C ２次元の弾性散乱
         px = (p(1,i1) + p(1,i2))/2
         py = (p(2,i1) + p(2,i2))/2
         pr = sqrt((p(1,i1)-p(1,i2))**2+(p(2,i1)-p(2,i2))**2)/2
C 散乱の方向は乱数で選ぶ
         theta = 2*pi*rnd(iseed)
         p(1,i1) = px + pr*cos(theta)
         p(2,i1) = py + pr*sin(theta)
         p(1,i2) = px - pr*cos(theta)
         p(2,i2) = py - pr*sin(theta)
      end do
C
C 分布を求めて出力 ----------------------------------------------------*
      energy = 0.0d0
      do i = 1, n
C Px を dx ごとの短冊に割り当てる。
         k = nint(p(1,i)/dx)
          if(k.ge.-50.and.k.le.50) then
             hist(k) = hist(k) + 1.0d0
          endif
         energy = energy + (p(1,i)**2 + p(2,i)**2)/2
      end do
      T = energy / n
C 結果出力
      do j = -50, 50
         x = dx*j
         write(16,100) x, hist(j)/n/dx, exp(-x**2/2/T)/sqrt(2*pi*T)
      end do
C フォーマット文
 100  format(3f10.4)
      end
