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

      subroutine genpro(nx,ny,xstep,ystep,procit,propol,nvmipp,nvmapp,
     +rmin,rmax,inform)

C     This subroutine generates a location problem (see the Figure).
C     First, a regular grid with nx horizontal points and ny vertical 
C     points in the positive orthant is considered. The points of the 
C     grid start at the origin with an horizontal distance of xstep and 
C     a vertical distance of ystept. This grid will be the space at 
C     which the cities (represented by polygons) will be distributed. 
C     Before building the cities, an area (rectangle) of preservation 
C     where almost nothing can be done, is defined. This area of 
C     preservation will receive, after the construction of the cities, 
C     an hydraulic plant of energy generation (to supply the energy to 
C     the cities). Then, in the rest of the space, the cities are built. 
C     At each point of the grid (out of the central region) a city 
C     (represented by a polygon) will be built with probability procit. 
C     The definition of the polygon uses variables nvmipp, nvmapp, rmin 
C     and rmax in a way described in the genpol (generate polygon) 
C     subroutine. To transmit the energy from the plant to the cities, a 
C     tower inside each city and a tower inside the central region must 
C     be built. The objective of the problem is to determine the 
C     location of this towers in order to minimize the sum of the 
C     distances from each city tower to the central one. 
C
C     On Entry:
C
C     nx    integer,
C           number of horizontal points in the grid,
C
C     ny    integer,
C           number of vertical points in the grid,
C
C     xstep double precision,
C           horizontal distance between points of the grid,
C
C     ystep double precision,
C           vertical distance between points of the grid,
C
C     procit  double precision, 
C           probability of defining a city at point of the grid
C           (0 <= procit <= 1),
C
C     propol double precision, 
C           probability of the city been defined by a polygon (0 <= 
C           propol <= 1), if the city is not defined by a polygon then 
C           it is defined by a circle,
C
C     nvmipp integer,
C           parameter for the polygon generation (described in genpol 
C           subroutine),
C
C     nvmapp integer,
C           parameter for the polygon generation (described in genpol 
C           subroutine),
C
C     rmin  double precision,
C           parameter for the polygon generation (described in genpol 
C           subroutine),
C
C     rmax  double precision,
C           parameter for the polygon generation (described in genpol 
C           subroutine).
C
C     On output:
C
C     As described in the genpol subroutine, the output is saved in the 
C     polyg common block.


C     PARAMETERS
      integer ncmax
      parameter ( ncmax  =                630000 )
      integer npmax
      parameter ( npmax  =               1000000 )
      integer nvsmax
      parameter ( nvsmax =            13 * npmax )

C     SCALAR ARGUMENTS
      integer inform,nvmapp,nvmipp,nx,ny
      double precision procit,propol,rmax,rmin,xstep,ystep

C     COMMON SCALARS
      integer nc,np,totnvs
      double precision xdisp,xscal,ydisp,yscal

C     COMMON ARRAYS
      integer nvs(npmax)
      double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax),
     +        pol(nvsmax),radii(ncmax),vert(2*nvsmax)

C     LOCAL SCALARS
      integer i,j,k
      logical inscre,outcre,insrec
      double precision c,cx,cy,seed,vx,vy,lx,ly,ux,uy

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

C     EXTERNAL SUBROUTINES
      external genpol

C     COMMON BLOCKS
      common /circl/ ccent,radii,nc
      common /poly1/ nvs,np,totnvs
      common /poly2/ vert
      common /poly3/ edges
      common /poly4/ pcent
      common /poly5/ pol
      common /ellip/ xdisp,xscal,ydisp,yscal

      inform = 0

C     SEED FOR THE RANDOM GENERATION

      seed = 760013.0d0

C     DEFINE CENTRAL REGION

      xdisp = 0.20d0 * ( nx - 1 ) * xstep
      ydisp = 0.50d0 * ( ny - 1 ) * ystep

      xscal = 0.15d0 * ( nx - 1 ) * xstep
      yscal = 0.50d0 * ( ny - 1 ) * ystep

C     DEFINE CENTRAL POINT POLYGON

      nvs(1)    =       4

      lx        =   xdisp + 0.50d0 * xscal 
      ux        =   xdisp + 1.50d0 * xscal
      ly        =   ydisp - 0.75d0 * yscal
      uy        =   ydisp + 0.75d0 * yscal

      vert(1)   =      lx
      vert(2)   =      ly
      vert(3)   =      lx
      vert(4)   =      uy
      vert(5)   =      ux
      vert(6)   =      uy
      vert(7)   =      ux
      vert(8)   =      ly

      edges(1)  = - 1.0d0
      edges(2)  =   0.0d0
      edges(3)  =      lx
      edges(4)  =   0.0d0
      edges(5)  =   1.0d0
      edges(6)  = -    uy
      edges(7)  =   1.0d0
      edges(8)  =   0.0d0
      edges(9)  = -    ux
      edges(10) =   0.0d0
      edges(11) = - 1.0d0
      edges(12) =      ly

      np        =       1
      totnvs    =       4

      pcent(1)  = ( lx + ux ) / 2.0d0
      pcent(2)  = ( ly + uy ) / 2.0d0

      do i = 1,4
          pol(i) = 1
      end do

C     THESE THREE LINES BELOW SHOULD BE USED INSTEAD OF THE 27 LINES
C     ABOVE IF THE CENTRAL RECTANGLE WOULD NOT TO BE CONSIDERED.

C     nvs(1)    =       0
C     np        =       1
C     totnvs    =       0

      nc        =       0

C     DEFINE CITY-POLYGONS AND CITY-CIRCLES CENTERED AT THE GRID 
C     POINTS. THE CITIES MUST BE OUTSIDE THE CENTRAL REGION AND NOT
C     COMPLETELY INSIDE THE CASSINIAN CURVE

      do i = 0,nx - 1
          do j = 0,ny - 1

              cx = i * xstep
              cy = j * ystep

              if ( drand(seed) .le. propol ) then

                ! GENERATE A CITY-POLYGON

                  call genpol(cx,cy,nvmipp,nvmapp,rmin,rmax,seed,inform)

                  if ( inform .ne. 0 ) then
                      return
                  end if

                ! TEST WHETHER THE POLYGON IS IN THE INTERIOR, IT IS 
                ! OUTSIDE OR IT CUTS CENTRAL REGION

                ! 1) CITIES IN THE INTERIOR ARE FORBIDDEN TO AVOID 
                ! INFEASIBLE PROBLEMS

                ! 2) CITIES THAT TOUCH THE CENTRAL REGION ARE DESIRED
                ! WITH THE AIM OF MAKING THE CENTRAL-REGION CONSTRAINT
                ! ACTIVE AT THE SOLUTION AS MANY TIMES AS POSSIBLE

                ! 3) CITIES THAT DO NOT TOUCH THE CENTRAL REGION ARE
                ! ARE INTRODUCED INTO THE PROBLEM WITH PROBABILITY PROCIT

                ! IT IS NOT EASY TO DETERMINE THE PREVIOUS RELATIONS (INSIDE,
                ! OUTSIDE, HALF AND A HALF) BETWEEN THE POLYGON AND THE 
                ! NON-C0NVEX CENTRAL REGION. SO WE WILL CONSIDER THAT:

                ! A POLYGON IS INSIDE THE CENTRAL REGION IF ALL ITS VERTICES
                ! ARE INSIDE OR IN THE BORDER OF THE CENTRAL REGION.

                ! A POLYGON IS OUTSIDE THE CENTRAL REGION IF ALL ITS VERTICES
                ! ARE OUTSIDE OR IN THE BORDER OF THE CENTRAL REGION.

                ! A POLYGON CUTS THE CENTRAL REGION IF IT HAS AT LEAST A 
                ! VERTEX INSIDE OR IN THE BORDER OF THE CENTRAL REGION AND
                ! AT LEAST A VERTEX OUTSIDE OR IN THE BORDER OF THE CENTRAL
                ! REGION. SO, TO BE HALF AND HALF IS EQUIVALENT TO NOT TO
                ! INSIDE, NOR OUTISDE.
                 
                  inscre = .true.
                  outcre = .true.

                  do k = 1,nvs(np+1)

                      vx = ( vert(2*(totnvs+k)-1) - xdisp ) / xscal
                      vy = ( vert(2*(totnvs+k))   - ydisp ) / yscal 

                      c = vx ** 2 + vy ** 2 - 1.0d0

                      if ( c .gt. 0.0d0 ) then
                          inscre = .false.
                      end if

                      if ( c .lt. 0.0d0 ) then
                          outcre = .false.
                      end if

                  end do

                  if ( cx .ge. lx - rmax .and. cx .le. ux + rmax .and.
     +                 cy .ge. ly - rmax .and. cy .le. uy + rmax ) then
                      insrec = .true.
                  else
                      insrec = .false.
                  end if

                ! ADD THE CITY-POLYGON TO THE PROBLEM
 
                  if ( .not. insrec .and.
     +               ( ( .not. inscre .and. .not. outcre ) .or.
     +                 ( outcre .and. drand(seed) .le. procit ) ) ) then

                      np = np + 1
                      totnvs = totnvs + nvs(np)

                  end if

              else

                ! GENERATE A CITY-CIRCLE

                  call gencir(cx,cy,rmin,rmax,seed,inform)

                  if ( inform .ne. 0 ) then
                      return
                  end if

                ! VERIFY WHERE IT IS

                  inscre = .true.
                  outcre = .true.

                  do k = 1,4

                      if ( k .eq. 1 ) then
                          vx = ( cx - radii(nc + 1) - xdisp ) / xscal
                          vy = ( cy                 - ydisp ) / yscal 
                      else if ( k .eq. 2 ) then
                          vx = ( cx + radii(nc + 1) - xdisp ) / xscal
                          vy = ( cy                 - ydisp ) / yscal 
                      else if ( k .eq. 3 ) then
                          vx = ( cx                 - xdisp ) / xscal
                          vy = ( cy - radii(nc + 1) - ydisp ) / yscal 
                      else if ( k .eq. 4 ) then
                          vx = ( cx                 - xdisp ) / xscal
                          vy = ( cy + radii(nc + 1) - ydisp ) / yscal 
                      end if

                      c = vx ** 2 + vy ** 2 - 1.0d0

                      if ( c .gt. 0.0d0 ) then
                          inscre = .false.
                      end if

                      if ( c .lt. 0.0d0 ) then
                          outcre = .false.
                      end if

                 end do

                 if ( cx .ge. lx - rmax .and. cx .le. ux + rmax .and.
     +                cy .ge. ly - rmax .and. cy .le. uy + rmax ) then
                     insrec = .true.
                 else
                     insrec = .false.
                 end if

                ! ADD THE CITY-CIRCLE TO THE PROBLEM

                  if ( .not. insrec .and.
     +               ( ( .not. inscre .and. .not. outcre ) .or.
     +                 ( outcre .and. drand(seed) .le. procit ) ) ) then

                      nc = nc + 1

                  end if

              end if

          end do
      end do

      end

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

      subroutine gencir(cx,cy,rmin,rmax,seed,inform)

C     PARAMETERS
      integer ncmax
      parameter ( ncmax  =                630000 )

C     SCALAR ARGUMENTS
      integer inform
      double precision cx,cy,rmax,rmin,seed

C     COMMON SCALARS
      integer nc

C     COMMON ARRAYS
      double precision ccent(2*ncmax),radii(ncmax)

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

C     COMMON BLOCKS
      common /circl/ ccent,radii,nc

      inform = 0

C     VERIFY SPACE AVAILABILITY FOR CIRCLES

      if ( nc .eq. ncmax ) then
          write (*,fmt=*) 'THE MAXIMUM NUMBER OF CIRCLES WAS ACHIEVED.'
          write (*,fmt=*) 'INCREASE NCMAX OR REDUCE THE NUMBER OF GRID '
          write (*,fmt=*) 'POINTS (NX*NY) OR THE PROBABILITY OF HAVING '
          write (*,fmt=*) 'A CITY-CIRCLE AT A GRID POINT'
          write (*,fmt=*) '(PROCIT*(1-PROPOL)).'
          inform = - 1
          return
      end if

C     SAVE THE CENTER 

      ccent(2 * nc + 1) = cx
      ccent(2 * nc + 2) = cy

C     GENERATE THE RADIUS

      radii(nc + 1) = rmin + ( rmax - rmin ) * drand(seed)

      end

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

      subroutine genpol(cx,cy,nvmipp,nvmapp,rmin,rmax,seed,inform)

C     This subroutine generates a polygon in R^2 with its 
C     vertices in a sphere centered at point (cx,cy). The
C     number of vertices is randomly generated 
C     satisfying nvmipp <= number of vertices  <= 
C     nvmapp. The ratio of the sphere is also randomly 
C     generated satisfying rmin <= ratio < rmax. The 
C     generated polygon is stored in the common block "polyg".
C 
C     On Entry:
C
C     cx    double precision,
C           first coordinate of the center of the sphere,
C
C     cy    double precision,
C           second coordinate of the center of the sphere,
C
C     nvmipp integer,
C           minimum number of vertices,
C
C     nvmapp integer,
C           maximum number of vertices,
C
C     rmin double precision,
C           minimum ratio of the sphere,
C
C     rmax double precision,
C           maximum ratio of the sphere,
C
C     seed double precision,
C           seed for the random generation.
C
C     On Output:
C
C     inform integer
C            termination parameter
C
C            0 = everything was ok
C
C          < 0 = there was not enough memory space.
C
C     The generated polygon is stored in the polyg common block
C     described below.
C
C     Common block polyg:
C
C     common /polyg/nvs,vert,edges,np,totnvs
C 
C     This structure represents, at any time, np polygons.
C     Position i of array nvs indicates the number of vertices
C     of polygon i. Arrays vert and edges store the 
C     vertices and edges of the polygons.
C
C     For example, if nvs(1) = 3 it indicates that the first
C     polygon has 3 vertices (edges). Then, if the vertices
C     are (x1,y1), (x2,y2) and (x3,y3), we have that vert(1) 
C     = x1, vert(2) = y1, vert(3) = x2, vert(4) = y2, vert(5) 
C     = x3, and vert(6) = y3. And, if the edges (written as
C     ax + by + c = 0) are a1 x + b1 y + c1 = 0, a2 x + b2 y 
C     + c2 = 0, and a3 x + b3 y + c3 = 0 then edges(1) = a1,
C     edges(2) = b1, edges(3) = c1, edges(4) = a2, edges(5) = 
C     b2, edges(6) = c2, edges(7) = a3, edges(8) = b3 and
C     edges(9) = c3.
C
C     totnvs indicates the total number of vertices  
C     of the set of polygons. This information is used when
C     a new polygon is created to know the first free position
C     of arrays vert and edges at which the vertices and edges 
C     of the new polygon will be saved.
C
C     Two additional details: 
C
C     1) For each polygon, the vertices are ordered clockwise 
C     and edge i corresponds to the edge between vertices 
C     i and i+1 (0 if i=n).
C
C     2) For each edge of the form ax + bx + c = 0, constants
C     a, b and c are chosen in such a way that 
C     (|a| = 1 or |b| = 1) and (a cx + b cy + c <= 0).

C     PARAMETERS
      integer npmax
      parameter ( npmax  =               1000000 )
      integer nvsmax
      parameter ( nvsmax =            13 * npmax )
      double precision pi
      parameter ( pi     =     3.1415926535898d0 )

C     SCALAR ARGUMENTS
      integer inform,nvmapp,nvmipp
      double precision cx,cy,rmax,rmin,seed

C     COMMON SCALARS
      integer np,totnvs

C     COMMON ARRAYS
      integer nvs(npmax)
      double precision edges(3*nvsmax),pcent(2*npmax),pol(nvsmax),
     +        vert(2*nvsmax)

C     LOCAL SCALARS
      integer i
      double precision dist,lseed,mindist,r

C     LOCAL ARRAYS
      double precision angl(nvsmax)

C     EXTERNAL FUNCTIONS
      double precision drand
      external drand

C     EXTERNAL SUBROUTINES
      external class,constr

C     INTRINSIC FUNCTIONS
      intrinsic cos,int,sin

C     COMMON BLOCKS
      common /poly1/ nvs,np,totnvs
      common /poly2/ vert
      common /poly3/ edges
      common /poly4/ pcent
      common /poly5/ pol

      inform = 0

C     VERIFY SPACE AVAILABILITY FOR POLYGONS

      if ( np .eq. npmax ) then
          write (*,fmt=*) 'THE MAXIMUM NUMBER OF POLYGONS WAS ACHIEVED.'
          write (*,fmt=*) 'INCREASE NPMAX OR REDUCE THE NUMBER OF GRID '
          write (*,fmt=*) 'POINTS (NX*NY) OR THE PROBABILITY OF HAVING '
          write (*,fmt=*) 'A CITY-POLYGON AT A GRID POINT'
          write (*,fmt=*) '(PROCIT*PROPOL).'
          inform = - 1
          return
      end if

C     GENERATE THE NUMBER OF VERTICES 

      lseed = 157318.0d0 + seed

      nvs(np+1) = nvmipp + int((nvmapp-nvmipp+1)*drand(lseed))

C     VERIFY SPACE AVAILABILITY FOR VERTICES AND SIDES

      if ( totnvs + nvs(np+1) .gt. nvsmax ) then
          write (*,fmt=*) 'THE MAXIMUM NUMBER OF POLYGONS VERTICES WAS'
          write (*,fmt=*) 'ACHIEVED. INCREASE NVSMAX OR REDUCE THE'
          write (*,fmt=*) 'AVERAGE NUMBER OF VERTICES (NVMAPP-NVMIPP)/2'
          write (*,fmt=*) 'OR THE NUMBER OF POLYGONS. TO REDUCE THE'
          write (*,fmt=*) 'NUMBER OF POLYGONS, REDUCE THE NUMBER OF'
          write (*,fmt=*) 'GRID POINTS (NX*NY) OR THE PROBABILITY OF'
          write (*,fmt=*) 'HAVING A CITY-POLYGONS AT A GRID POINT'
          write (*,fmt=*) '(PROCIT*PROPOL).'
          inform = - 1
          return
      end if

C     SAVE THE "CENTER" TO BE USED AS INITIAL POINT

      pcent(2*np+1)  = cx
      pcent(2*np+2)  = cy

C     IDENTIFY THE CONSTRAINTS WITH THE POLYGON

      do i = totnvs + 1,totnvs + nvs(np+1)
          pol(i) = np + 1
      end do

C     GENERATE THE RADIUS OF THE SPHERE

      r = rmin + ( rmax - rmin ) * drand( seed )

C     GENERATE ALL ANGLES SATISFYING 0 <= ANGLE_I < 2*PI

 10   continue

      do i = 1,nvs(np + 1)
          angl(i) = 2 * pi * drand( lseed )
      end do

C     CLASSIFY THE ANGLES IN DECREASING ORDER

      call class(nvs(np + 1),angl)

C     CONSTRUCT THE VERTICES

      do i = 1,nvs(np + 1)
          vert(2 * ( totnvs + i ) - 1) = cx + r * cos( angl(i) )
          vert(2 * ( totnvs + i )    ) = cy + r * sin( angl(i) )
      end do

C     FOR NUMERICAL STABILITY IN THE CONSTRAINT CALCULATION,
C     AVOID TOO SIMILAR ANGLES (EQUIVALENT TO TOO NEAR POINTS)

      mindist = 1.0D+99

      do i = totnvs + 1,totnvs + nvs(np+1) - 1

          dist = ( vert(2*i-1) - vert(2*i+1) ) ** 2 +
     +           ( vert(2*i)   - vert(2*i+2) ) ** 2

          if ( dist .lt. mindist ) then
              mindist = dist
          end if

      end do

      i = totnvs + nvs(np+1)

      dist = ( vert(2*i-1) - vert(2*totnvs+1) ) ** 2 +
     +       ( vert(2*i)   - vert(2*totnvs+2) ) ** 2

      if ( dist .lt. mindist ) then
          mindist = dist
      end if

      if ( mindist .lt. 1.0D-8 ) then
C         write (*,fmt=*) 'DISCARDING ANGLES.'
          go to 10
      end if

C     CONSTRUCT THE EDGES

      do i = totnvs + 1,totnvs + nvs(np+1) - 2
          call constr(vert(2*i-1),vert(2*i),vert(2*i+1),vert(2*i+2),
     +                vert(2*i+3),vert(2*i+4),edges(3*i-2),edges(3*i-1),
     +                edges(3*i),inform)
          if ( inform .ne. 0 ) then
              return
          end if
      end do

      i = totnvs + nvs(np+1) - 1

      call constr(vert(2*i-1),vert(2*i),vert(2*i+1),vert(2*i+2),
     +            vert(2*totnvs+1),vert(2*totnvs+2),edges(3*i-2),
     +            edges(3*i-1),edges(3*i),inform)
      if ( inform .ne. 0 ) then
          return
      end if

      i = totnvs + nvs(np+1)

      call constr(vert(2*i-1),vert(2*i),vert(2*totnvs+1),
     +            vert(2*totnvs+2),vert(2*totnvs+3),vert(2*totnvs+4),
     +            edges(3*i-2),edges(3*i-1),edges(3*i),inform)
      if ( inform .ne. 0 ) then
          return
      end if

      end

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

      subroutine constr(x1,y1,x2,y2,x3,y3,a,b,c,inform)

C     This subroutine computes the real constants a, b and c of 
C     the straight line ax + by + c = 0 in R^2 defined by the 
C     points (x1,y1) and (x2,y2); such that the point (x3,y3) 
C     satisfies the constraint ax + by + c <= 0.
C
C     On Entry:
C
C     x1    double precision,
C           first coordinate of point (x1,y1),
C
C     y1    double precision,
C           second coordinate of point (x1,y1),
C
C     x2    double precision,
C           first coordinate of point (x2,y2),
C
C     y2    double precision,
C           second coordinate of point (x2,y2),
C
C     x3    double precision,
C           first coordinate of point (x3,y3),
C
C     y3    double precision,
C           second coordinate of point (x3,y3).
C
C     On Return
C
C     a,b,c double precision
C           the desired constants.

C     SCALAR ARGUMENTS
      integer inform
      double precision a,b,c,x1,x2,x3,y1,y2,y3

      if (x1.eq.x2 .and. y1.eq.y2) then
          write (*,fmt=*)
     +      'ERROR IN FUNCTION CONSTRAINT: X1=X2 AND Y1=Y2'
          inform = - 1
          return
      end if

      if (y1.ne.y2) then
          a = 1.0d0
          b = - (x2-x1)/ (y2-y1)
          c = - (x1+b*y1)

      else
          a = 0.0d0
          b = 1.0d0
          c = -y1
      end if

      if (a*x3+b*y3+c.gt.0.0d0) then
          a = -a
          b = -b
          c = -c
      end if

      end

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

      subroutine class(n,x)

C     This subroutine classifies the elements of a vector in
C     decreasing order, i.e., on output: x(1) >= x(2) >= 
C     ... >= x(n).
C
C     On Entry:
C
C     n     integer,
C           number of elements of the vector to be classified,
C
C     x     double precision x(n),
C           vector to be classified.
C
C     On Return
C
C     x     double precision x(n),
C           classified vector.

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      double precision aux,xmax
      integer i,j,pos

      do i = 1,n

          xmax = x(i)
          pos = i

          do j = i + 1,n
              if (x(j).gt.xmax) then
                  xmax = x(j)
                  pos = j
              end if
          end do

          if (pos.ne.i) then
              aux = x(i)
              x(i) = x(pos)
              x(pos) = aux
          end if

      end do

      end
