PROGRAM scatter
*     微分断面積と全断面積の計算
      IMPLICIT NONE
      REAL*8 N, R, b, dN(0:360), delta_b, dt, dtheta, vx, vx0,
     +       vy, x0
      INTEGER IR_, nbin
*     LOGICAL TBkeyinput
      DATA dN/361*0.0D0/
      CALL GWopen(IR_, 0)
      CALL initial(x0,vx0,R,delta_b,dt,nbin,dtheta)
      N = 0                                 ! 入射粒子の数
      b = 0
*     b = R としていても，b が R をわずかに越える可能性があることに注意
      DO WHILE(b .LE. R + delta_b/2)
         CALL trajectory(vx,vy,vx0,x0,b,dt)
         CALL detector(dN,vx,vy,b,dtheta,N)
         b = b + delta_b
      END DO
*     WRITE(*,*) 'Hit any key for the differential cross section.'
*     DO WHILE(.NOT.TBkeyinput())
*        CALL GWsleep(IR_, 10)
*     END DO
      CALL output(dN,N,nbin,dtheta,R)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(x0,vx0,R,delta_b,dt,nbin,dtheta)
      IMPLICIT NONE
      REAL*8 E, x0, vx0, R, delta_b, dt, dtheta
      REAL Rs, xwin, ywin
      INTEGER IR_, nbin
      WRITE(*,'(A,$)') 'kinetic energy of incident particles = '
      READ(*,*) E
      WRITE(*,'(A,$)') 'incremental change in impact parameter = '
      READ(*,*) delta_b
      R = 1.0D0                             ! ビームの半径
      vx0 = sqrt(2*E)                       ! 質量 = 1
      dt = 0.01D0                           ! 分割時間
      nbin = 9                              ! 検出器の数
      dtheta = 180.0D0/nbin                 ! 角度
*     目標物から測ったビームの初期位置
      x0 = -3
      CALL compute_aspect_ratio(REAL(abs(x0)),xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      CALL GWsetpen(IR_, 16, -1, -1, -1)
      Rs = 1.0                              ! 散乱領域の半径
      CALL GWellipse(IR_, -Rs, -Rs, Rs, Rs)
      CALL GWsetpxl(IR_, 0.0, 0.0, 16)
      END
C
      SUBROUTINE trajectory(vx,vy,vx0,x0,b,dt)
      IMPLICIT NONE
      REAL*8 x, y, vx, vy, vx0, x0, b, dt
      INTEGER IR_
      LOGICAL Ltmp1_
      y = b
      x = x0
      vx = vx0
      vy = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL EulerRichardson(x,y,vx,vy,dt) ! 分割時間の1ステップ
         CALL GWsetpxl(IR_, REAL(x), REAL(y), 13)
         Ltmp1_ = x*x + y*y .LE. x0*x0 + b*b
      END DO
      END
C
      SUBROUTINE detector(dN,vx,vy,b,dtheta,N)
      IMPLICIT NONE
      REAL*8 dN(0:360), vx, vy, b, dtheta, N, arctan
      INTEGER theta
      theta = int(arctan(vx,vy)/dtheta)     ! 散乱角
      dN(theta) = dN(theta) + b
      N = N + b
      END
C
      SUBROUTINE EulerRichardson(x,y,vx,vy,dt)
*     配列を使うとプログラムの行数を減らすことができる
      IMPLICIT NONE
      REAL*8 fx, fxmid, fy, fymid, vxmid, vymid, xmid, ymid, x, y, vx,
     +       vy, dt
      CALL force(x,y,fx,fy)
*     プログラム中では質量は1である
      vxmid = vx + 0.5D0*fx*dt
      vymid = vy + 0.5D0*fy*dt
      xmid = x + 0.5D0*vx*dt
      ymid = y + 0.5D0*vy*dt
      CALL force(xmid,ymid,fxmid,fymid)
      vx = vx + fxmid*dt
      vy = vy + fymid*dt
      x = x + vxmid*dt
      y = y + vymid*dt
      END
C
      SUBROUTINE force(x,y,fx,fy)
      IMPLICIT NONE
      REAL*8 f, r, r2, x, y, fx, fy
      r2 = x*x + y*y
      r = sqrt(r2)
*     簡単な水素原子のモデルの場合の相互作用
*     力の作用する範囲を 1 に設定
      IF(r2 .LE. 1) THEN
         f = 1/r2 - r
      ELSE
         f = 0
      END IF
      fx = f*x/r
      fy = f*y/r
      END
C
      SUBROUTINE output(dN,N,nbin,dtheta,R)
      IMPLICIT NONE
      REAL*8 d_Omega, diffcross, dtheta_rad, pi, target_density, theta,
     +       totalcross, dN(0:360), N, dtheta, R
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER i, nbin
*     INTEGER IR_
      DATA totalcross/0.0D0/
*     CALL GWclear(IR_, -1)
      dtheta_rad = dtheta*pi/180            ! ラジアン
      target_density = 1/(pi*R*R)
      WRITE(*,*)
      WRITE(*,*) '       theta         dN/N   sigma(theta)'
      WRITE(*,*)
      DO i = 0, nbin
         theta = (i + 0.5D0)*dtheta_rad
         IF(dN(i) .GT. 0) THEN
            d_Omega = 2*pi*sin(theta)*dtheta_rad
            diffcross = dN(i)/(d_Omega*N*target_density)
            totalcross = totalcross + diffcross*d_Omega
            WRITE(*,'(3F13.6)') (i+1)*dtheta, dN(i)/N, diffcross
         END IF
      END DO
      WRITE(*,*)
      WRITE(*,*) 'total cross section =', totalcross
      END
C
      FUNCTION arctan(x,y)
*     -pi/2 ラジアンから pi/2 ラジアンを 0 度から 180 度に変換
      IMPLICIT NONE
      REAL*8 arctan, pi, theta, x, y
      PARAMETER(pi = 3.141592653589793D0)
      IF(x .GT. 0) THEN
         theta = atan(abs(y/x))
      ELSE IF(y .GT. 0) THEN
         theta = atan(y/x) + pi
      ELSE
         theta = pi - atan(y/x)
      END IF
      arctan = theta*180/pi
      END