program trisectrix of maclaurin
      implicit none
      integer pgopen, i
      real x1(3600), x2(3600), r, theta, delta_theta, a, pi

      if (pgopen('/xwin').le.0) stop

      write(*,*) 'Please input a non-zero real numbers "a" (ex. 1.0):'
      read(*,*) a

      pi = 3.14159268
c     delta_theta = 2 * 2*pi/3600
      delta_theta = 3 * 2*pi/3600
      
      do i=1, 3600
        theta = (i-1)*delta_theta
        r = a/cos(theta/3)
        x1(i)=r*cos(theta)
        x2(i)=r*sin(theta)
      end do

      call pgenv(-4.0,4.0,-4.0,4.0,1,0)
      call pglab('x1','x2','Trisectrix of Maclaurin')
      call pgline(3600,x1,x2)
      call pgclos

      end