PROGRAM interfere
      IMPLICIT NONE
      REAL*8 L, a, lambda, ymax
      INTEGER IR_, nbar, nslit
      CALL GWopen(IR_, 0)
      CALL initial(nslit,a,L,lambda,ymax,nbar)
      CALL pattern(nslit,a,L,lambda,ymax,nbar)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(nslit,a,L,lambda,ymax,nbar)
      IMPLICIT NONE
      REAL*8 Imax, Lt, dy, a, L, lambda, ymax
      INTEGER IR_, nslit, nbar, itick, ntick
      WRITE(*,'(A,$)') 'number of slits = '
      READ(*,*) nslit
      WRITE(*,'(A,$)') 'distance between slits (mm) = '
      READ(*,*) a
      WRITE(*,'(A,$)') 'distance to the screen (mm) = '
      READ(*,*) L
      WRITE(*,'(A,$)') 'wavelength (angstroms)= '
      READ(*,*) lambda
      WRITE(*,'(A,$)') 'max. position on photographic plate (mm) = '
      READ(*,*) ymax
      WRITE(*,'(A,$)') '# times to average intensity = '
      READ(*,*) nbar
      lambda = lambda*1d-7                  ! オングストロームを mm に変換
      Imax = (nslit/L)**2                   ! 強度の最大値
      CALL GWindow(IR_, REAL(-1.1D0*ymax), REAL(-0.1D0*Imax)
     +, REAL(1.1D0*ymax), REAL(1.1D0*Imax))
      CALL GWline(IR_, REAL(-ymax), 0.0, REAL(ymax), 0.0) ! 横軸を描く
      CALL GWline(IR_, 0.0, REAL(Imax), 0.0, 0.0) ! 強度に対する軸を描く
      WRITE(*,*) 'intensity versus vertical distance on screen'
      dy = lambda*L/a                       ! L >> aならば，極大値間の距離
      ntick = int(ymax/dy)
      Lt = 0.01D0*Imax                      ! 縦軸の刻み
      DO itick = -ntick, ntick
         CALL GWline(IR_, REAL(itick*dy), 0.0, REAL(itick*dy), REAL(Lt))
      END DO
      END
C
      SUBROUTINE pattern(nslit,a,L,lambda,ymax,nbar)
      IMPLICIT NONE
      REAL*8 L2, amplitude, delta, dy, k, phase,
     +       pi, r, r2, y, yslit, a, L, lambda, ymax, intensity
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, nslit, nbar, ipoint, islit, npoint, t
      CALL GWsetpen(IR_, -1, -1, -1, -1);
      k = 2*pi/lambda                       ! 波数
      npoint = 200                          ! y 軸の正の領域にある点の数
      dy = ymax/npoint                      ! プロットする点の間の距離
      delta = a/(nslit - 1)                 ! スリット間の距離
      phase = 2*pi/nbar
      L2 = L*L
      DO ipoint = -npoint, npoint
         y = ipoint*dy                      ! スクリーンの位置
         intensity = 0                      ! 光の強度      
         DO t = 1, nbar
            amplitude = 0
            DO islit = 1, nslit
               yslit =  -0.5D0*a + (islit - 1)*delta ! スリットの位置
               yslit = yslit - y            ! スリットとスクリーンの位置の縦の距離
               r2 = L2 + yslit*yslit
               r = sqrt(r2)                 ! スリットのスクリーンからの距離
               amplitude = (1/r)*cos(k*r - phase*t) + amplitude
            END DO
            intensity = intensity + amplitude*amplitude
         END DO
         intensity = intensity/nbar
         CALL GWline2(IR_, REAL(y), REAL(intensity))
      END DO
      END