C     ******************************************************************
C     ******************************************************************

      program latticema

      implicit none

C     PARAMETER
      integer mmax
      double precision pi
      parameter ( mmax = 2000 )
      parameter ( pi   = 3.1415926535898d0 )

C     LOCAL SCALARS
      integer i,ignum,m
      double precision a,ang,b,ddummy,r,seed

C     LOCAL ARRAYS
      double precision u(mmax),v(mmax),s(mmax)

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

c     write(*,*) 'Enter a and b (a >= b): '
c     read(*,*) a,b

c     write(*,*) 'Enter r: '
c     read(*,*) r

c     write(*,*) 'Enter m: '
c     read(*,*) m

      a = 30.0d0
      b = 15.0d0
      r = 1.0d0
      m = 375

      write(*,*) 'Enter type of random initial guess (from 1 to 4): '
      read(*,*) ignum

      seed = 123456.0d0
      ddummy = drand( seed )

      if ( ignum .eq. 1 ) then
          do i = 1,m
              ang = 2.0d0 * pi * drand( seed )
              u(i) = - a + 2.0d0 * a * drand( seed )
              v(i) = - b + 2.0d0 * b * drand( seed )
              s(i) = drand( seed )
          end do

      else if ( ignum .eq. 2 ) then
          do i = 1,m
              ang = 2.0d0 * pi * drand( seed )
              u(i) = - 0.5d0 * a + a * drand( seed )
              v(i) = - 0.5d0 * b + b * drand( seed )
              s(i) = drand( seed )
          end do

      else if ( ignum .eq. 3 ) then
          do i = 1,m
              ang = 2.0d0 * pi * drand( seed )
              u(i) = a * cos( ang )
              v(i) = b * sin( ang )
              s(i) = drand( seed )
          end do

      else if ( ignum .eq. 4 ) then
          do i = 1,m
              ang = 2.0d0 * pi * drand( seed )
              u(i) = a * cos( ang )
              v(i) = b * sin( ang )
              s(i) = 0.5d0 * drand( seed )
          end do
      end if
 
      call drawsol(a,b,r,m,u,v,s)

      stop

      end

C     ******************************************************
C     ******************************************************

      subroutine drawsol(a,b,r,m,u,v,s)

      implicit none 

C     SCALAR ARGUMENTS
      integer m
      double precision a,b,r

C     ARRAY ARGUMENTS
      double precision u(m),v(m),s(m)

C     LOCAL SCALARS
      integer i
      double precision scale,x,y

      open(unit=10,file='iniguess.mp')

C     SCALING
      scale = 10.0d0 / max( 2.0d0 * a, 2.0d0 * b )
      write(10,10) scale

C     ELLIPSE
      write(10,20) 2.0d0*a,2.0d0*b,0.0d0,0.0d0

C     CIRCULAR ITEMS
      do i = 1,m
          x = ( 1.0d0 + ( s(i) - 1.0d0 ) * ( b**2 / a**2 ) ) * u(i)
          y = s(i) * v(i)

          write(10,20) 2.0d0*r,2.0d0*r,x,y
      end do

      write(10,30)

      close(10)

C     NON-EXECUTABLE STATEMENTS

  10  format('beginfig(-1);'/,'u = ',f20.10,'cm;')
  20  format('draw fullcircle',
     +       /,'     xscaled  ',f20.10,'u',
     +       /,'     yscaled  ',f20.10,'u',
     +       /,'     shifted (',f20.10,'u,',f20.10,'u);')
  30  format('endfig;',/,'end;') 

      end

C     ******************************************************************
C     ******************************************************************

      double precision function drand(ix)

C     This is the random number generator of Schrage:
C
C     L. Schrage, A more portable Fortran random number generator, ACM
C     Transactions on Mathematical Software 5 (1979), 132-138.

      double precision ix

      double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k
      data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/

      xhi= ix/b16
      xhi= xhi - dmod(xhi,1.d0)
      xalo= (ix-xhi*b16)*a
      leftlo= xalo/b16
      leftlo= leftlo - dmod(leftlo,1.d0)
      fhi= xhi*a + leftlo
      k= fhi/b15
      k= k - dmod(k,1.d0)
      ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k
      if (ix.lt.0) ix= ix + p
      drand= ix*4.656612875d-10

      return

      end
