implicit none
       double precision :: theta_a, theta_b, theta_c, theta_d
       double precision :: tb, tc, td, alpha, freq, omega, pi, t, theta
       double precision :: wg, hg, r0, rad2deg, deg2rad, area(10000), a0
       integer :: nt, i
c
       pi = 2. * asin(1.)
       rad2deg = 180. / pi
       deg2rad = pi / 180.
c
       write(*,1,advance='no')
  1    format(/,5x,'Enter the chopper frequency :: ')
       read(*,*) freq
       omega = 2. * pi * freq
       write(*,2,advance='no')
  2    format(/,5x,'Enter the width, height and R0 :: ')
       read(*,*) wg, hg, r0
c
       theta_a = -1. * atan(0.5 * wg / (r0 - 0.5 * hg))
       theta_b = -1. * atan(0.5 * wg / (r0 + 0.5 * hg))
       theta_c = atan(0.5 * wg / (r0 + 0.5 * hg))
       theta_d = atan(0.5 * wg / (r0 - 0.5 * hg))
       a0 = wg * hg
c
       tb = (theta_b - theta_a) / omega
       tc = (theta_c - theta_a) / omega
       td = (theta_d - theta_a) / omega
       write(*,*) (1.e6*tb), (1.e6*tc), (1.e6*td)
c
       write(*,3) (rad2deg * (theta_d - theta_a))
  3    format(/,5x,'The angle (theta_d - theta_a) = ',f6.2,' degs')
       write(*,4,advance='no')
  4    format(/,5x,'Enter the alpha angle :: ')
       read(*,*) alpha
       alpha = deg2rad * alpha
c
       if(alpha.lt.(theta_d - theta_a)) then
         write(*,5)
  5      format(/,5x,'Error :: alpha < (theta_d - theta_a)')
         stop
       end if
       nt = nint(1000000. * (abs(theta_a) + 0.5 * alpha) / omega)
       do i = 1, (nt+1)
         t = real(i-1)/1000000.
         theta = omega * t + theta_a
         if(theta.lt.theta_b) then
           area(i) = 0.5 * (0.5 * wg + (r0 - 0.5 * hg) * tan(theta)) *
     &                  ((-0.5 * wg / tan(theta)) - (r0 - 0.5 * hg))
         end if
         if(theta.ge.theta_b.and.theta.lt.(0.0)) then
           area(i) = ( 0.5 * wg + r0 * tan(theta)) * hg
         end if
         if(theta.ge.(0.0).and.theta.lt.theta_c) then
           area(i) = ( 0.5 * wg + r0 * tan(theta)) * hg
         end if
         if(theta.ge.theta_c.and.theta.lt.theta_d) then
           area(i) = a0 - 0.5 * (0.5 * wg - (r0 - 0.5 *hg) *
     &          tan(theta)) * (0.5 * wg / tan(theta) - (r0 - 0.5 *hg))
         end if
         if(theta.ge.theta_d) then
           area(i) = a0
         end if
         area(i) = area(i) / a0
       end do
       open(unit=42,file='t1b_open.dat',status='unknown')
       do i = 1, (nt+1)
         write(42,6) (i-1), area(i)
  6      format(5x,i5,5x,f6.4)
       end do
       do i = 1, nt
         write(42,6) (nt + i), area(nt + 1 - i)
       end do
       close(unit=42)
       stop
       end