      program drawsolma

      implicit none 

C     PARAMETERS
      integer mmax
      parameter ( mmax = 100 )

C     LOCAL SCALARS
      integer m
      double precision a,b,r

C     LOCAL ARRAYS
      double precision u(mmax),v(mmax),s(mmax)
      character * 30 filename 

      r = 1.0d0

      a = 2.0
      b = 1.0

      m = 2

      u(1) =  1.8648489729625890D+00
      u(2) = -1.4833079891008683D+00
      v(1) = -9.7868617995693452D-01
      v(2) = -7.7862651730129495D-01
      s(1) =  3.1206195674682435D-01
      s(2) =  8.2530732067849288D-01

      filename = "solution-a2-b1-m2.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 2.0
      b = 1.0

      m = 3

      u(1) =  1.8648489729625890D+00
      u(2) = -1.4833079891008683D+00
      u(3) = -1.9573723599138690D+00
      v(1) = -7.7862651730129495D-01
      v(2) = -3.7587608650635129D-01
      v(3) =  6.5061464135698577D-01
      s(1) =  9.4013937554954130D-01
      s(2) =  9.2248569513574175D-01
      s(3) =  2.1707896478123603D-01

      filename = "solution-a2-b1-m3.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 2.0
      b = 2.0

      m = 2

      u(1) =  1.9995171083079881D+00
      u(2) = -1.9995171077622456D+00
      v(1) =  4.3945307825742269D-02
      v(2) = -4.3945312500000000D-02
      s(1) =  5.0000001295669760D-01
      s(2) =  5.0000001299647479D-01

      filename = "solution-a2-b2-m2.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 2.0
      b = 2.0

      m = 3

      u(1) =  1.5228262083996769D+00
      u(2) = -2.0000000000000000D+00
      u(3) =  6.2619358443252182D-01
      v(1) = -1.5225347986836257D+00
      v(2) = -6.0696487753092587D-01
      v(3) =  2.0000000000000000D+00
      s(1) =  5.5653552809152507D-01
      s(2) =  5.3868141744650122D-01
      s(3) =  5.4018600489019997D-01

      filename = "solution-a2-b2-m3.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 3.0
      b = 2.0

      m = 2

      u(1) =  3.0000000000000000D+00
      u(2) = -3.0000000000000000D+00
      v(1) =  1.2836276281188849D-02
      v(2) = -1.2836294893363378D-02
      s(1) =  2.5003596543631057D-01
      s(2) =  2.5003596517676613D-01

      filename = "solution-a3-b2-m2.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 3.0
      b = 2.0

      m = 3

      u(1) =  2.9615662753969927D+00
      u(2) = -2.9794107357377499D+00
      u(3) =  3.8669170460519015D-01
      v(1) =  3.1916010135820849D-01
      v(2) =  2.3406923109497071D-01
      v(3) = -1.9833015131119722D+00
      s(1) =  2.6166544363498240D-01
      s(2) =  2.5634752364044311D-01
      s(3) =  4.9767166335181307D-01

      filename = "solution-a3-b2-m3.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 3.0
      b = 2.0

      m = 4

      u(1) =  3.5087895845462430D-01
      u(2) = -2.9977231315477422D+00
      u(3) = -3.5087895805073788D-01
      u(4) =  2.9977231314015400D+00
      v(1) = -1.9862558082961832D+00
      v(2) = -7.8351518595436181D-02
      v(3) =  1.9862558083278892D+00
      v(4) =  7.8351521067719657D-02
      s(1) =  4.9808366858181347D-01
      s(2) =  2.5072651428056730D-01
      s(3) =  4.9808366858623937D-01
      s(4) =  2.5072651432608611D-01

      filename = "solution-a3-b2-m4.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 4.0
      b = 2.0

      m = 2

      u(1) =  4.0000000000000000D+00
      u(2) = -4.0000000000000000D+00
      v(1) = -7.2061593439745738D-03
      v(2) =  7.2065119889902635D-03
      s(1) =  2.6355493674202789D-05
      s(2) =  2.6357675365711593D-05

      filename = "solution-a4-b2-m2.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 4.0
      b = 2.0

      m = 3

      u(1) =  3.9819011058168394D+00
      u(2) = -3.9962805728848982D+00
      u(3) =  1.3932752197159106D+00
      v(1) = -1.9004177347770981D-01
      v(2) = -8.7037461378811284D-02
      v(3) =  1.8747148923655830D+00
      s(1) =  1.3274379509982786D-02
      s(2) =  2.8498585630901741D-03
      s(3) =  4.7555888111584671D-01

      filename = "solution-a4-b2-m3.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 4.0
      b = 2.0

      m = 4

      u(1) =  3.9927534108054132D+00
      u(2) = -3.9927534559862314D+00
      u(3) =  1.3926235577070740D+00
      u(4) = -1.3926235411027119D+00
      v(1) = -1.2041903073767272D-01
      v(2) =  1.2041866656466509D-01
      v(3) =  1.8748678173549151D+00
      v(4) = -1.8748677887792791D+00
      s(1) =  5.3968833950756032D-03
      s(2) =  5.3968514842667572D-03
      s(3) =  4.7559373720362452D-01
      s(4) =  4.7559372552400414D-01

      filename = "solution-a4-b2-m4.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 2.0
      b = 1.0

      m = 4

      u(1) =  1.8648489729625890D+00
      u(2) = -1.4833079891008683D+00
      u(3) = -1.9573723599138690D+00
      u(4) = -1.5572530346025899D+00
      v(1) = -3.7587608650635129D-01
      v(2) =  6.5061464135698577D-01
      v(3) =  8.8027875109908260D-01
      v(4) =  8.4497139027148349D-01
      s(1) =  2.1707896478123603D-01
      s(2) =  4.4616127079147250D-01
      s(3) =  6.3247858805618473D-01
      s(4) =  6.7630021395144943D-02

      filename = "solution-a2-b1-m4.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 2.0
      b = 2.0

      m = 4

      u(1) =  1.6033282853841628D+00
      u(2) = -1.6006845901487263D+00
      u(3) =  1.6216492809849328D+00
      u(4) = -1.6219755334418720D+00
      v(1) = -1.6399633688024249D+00
      v(2) = -1.6412373851180029D+00
      v(3) =  1.6123646582589093D+00
      v(4) =  1.6165917200495119D+00
      s(1) =  6.1528077072452791D-01
      s(2) =  6.1147713035879303D-01
      s(3) =  6.0479443279070888D-01
      s(4) =  6.0631911997949084D-01

      filename = "solution-a2-b2-m4.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 4.0
      b = 2.0

      m = 5

      u(1) =  3.9855119351434056D+00
      u(2) = -3.9861007367856480D+00
      u(3) =  1.3930622116160876D+00
      u(4) =  5.8244362080668732D-01
      u(5) = -1.3930252252422035D+00
      v(1) =  1.7013731437284824D-01
      v(2) =  1.6666205531055994D-01
      v(3) = -1.8747856854962486D+00
      v(4) =  1.9786844041710712D+00
      v(5) = -1.8747912032869447D+00
      s(1) =  1.0684684571775264D-02
      s(2) =  1.0259861073582229D-02
      s(3) =  4.7557704618377367D-01
      s(4) =  4.9597668217944180D-01
      s(5) =  4.7557794664064218D-01

      filename = "solution-a4-b2-m5.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 3.0
      b = 2.0

      m = 5

      u(1) =  3.0000000000000000D+00
      u(2) = -2.1543743603892018D+00
      u(3) =  4.0729402875119608D-01
      u(4) =  9.9524221705778948D-01
      u(5) = -1.5148527479691261D+00
      v(1) = -9.3369384230012842D-01
      v(2) = -1.6723349942934891D+00
      v(3) = -1.9769455692047639D+00
      v(4) =  1.9451628090881612D+00
      v(5) =  1.9765303830381287D+00
      s(1) =  4.1739958836016722D-01
      s(2) =  5.0000000000000000D-01
      s(3) =  4.9598905398370380D-01
      s(4) =  5.0000000000000000D-01
      s(5) =  5.5316406332684365D-01

      filename = "solution-a3-b2-m5.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 2.0
      b = 2.0

      m = 5

      u(1) = -1.7420718438715945D+00
      u(2) = -1.9808628202554213D+00
      u(3) =  2.0000000000000000D+00
      u(4) =  2.0674445785123039D-02
      u(5) =  1.7554279082945392D+00
      v(1) =  1.3732015828332151D+00
      v(2) = -9.5070545493874781D-01
      v(3) = -9.6554170710190024D-01
      v(4) = -2.0000000000000000D+00
      v(5) =  1.3272728416289283D+00
      s(1) =  5.6014452557609506D-01
      s(2) =  1.0000000000000000D+00
      s(3) =  1.0000000000000000D+00
      s(4) =  4.8320827657172277D-01
      s(5) =  5.7476078145381582D-01

      filename = "solution-a2-b2-m5.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      a = 2.0
      b = 1.0

      m = 5

      u(1) =  1.8648489729625890D+00
      u(2) = -1.4833079891008683D+00
      u(3) = -1.9573723599138690D+00
      u(4) = -1.5572530346025899D+00
      u(5) = -7.5175217301270258D-01
      v(1) =  6.5061464135698577D-01
      v(2) =  8.8027875109908260D-01
      v(3) =  8.4497139027148349D-01
      v(4) = -5.6584207043752799D-01
      v(5) = -1.0767745841705501D-01
      s(1) =  6.3247858805618473D-01
      s(2) =  6.7630021395144943D-02
      s(3) =  6.5776964816414241D-01
      s(4) =  1.3447727827319586D-01
      s(5) =  1.5961605689559272D-01

      filename = "solution-a2-b1-m5.mp"
      call drawsol(a,b,r,m,u,v,s,filename)

      stop

      end

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

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

      implicit none 

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

C     ARRAY ARGUMENTS
      double precision u(m),v(m),s(m)
      character * 30 filename

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

      open(unit=10,file=trim(filename))

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

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;')
  11  format('pickup pencircle scaled 0.01;')
  20  format('draw fullcircle',
     +       /,'     xscaled  ',f20.10,'u',
     +       /,'     yscaled  ',f20.10,'u',
     +       /,'     shifted (',f20.10,'u,',f20.10,'u);')
  30  format('endfig;',/,'end;') 

      end
