program cube53
c
c     Mike Collins -- cinclodes@yahoo.com
c
c     This FORTRAN code simulates on the screen a 5x5x5 version of 
c     Rubik's Cube that contains a 3x3x3 cube. The inner cube is not
c     affected when the outer layer of a face is twisted. A face of 
c     the inner cube turns along with a face of the outer cube when
c     both layers are twisted. 
c
c     The output file cube53.ps contains the display. After compiling 
c     the code, run it and enter zeros when prompted. The code will 
c     quit after you enter the second zero. This initializes the 
c     display file, which you can now open with your Postscript 
c     displayer. Keep the display open and run the code again. This 
c     time, you can enter a random integer to mix up the code or zero 
c     again if you want to start with the puzzle solved. Then you can 
c     make a series of moves by entering one integer per line. Enter
c     negatives to twist in the opposite directions. You can reopen 
c     the display after each move.  
c
      real x(150,4),y(150,4),x3(54,4),y3(54,4)
      integer iclr(150),iclr3(54)
    1 open(unit=1,status='unknown',file='cube53.ps')
      call start()
c
      x0=6.0
      y0=7.5
      x03=6.0
      y03=4.2
      cube=0.20
      pi=4.0*atan(1.0)
      theta=45.0
      cost=cos(theta*pi/180.0)
      sint=sin(theta*pi/180.0)
c
      call setup(x,y,x0,y0,cube,cost,sint,iclr)
      call setup3(x3,y3,x03,y03,cube,cost,sint,iclr3)
c
      do 2 m=1,150
      n=iclr(m)
      call box(x,y,m,n)
    2 continue
      call frame(x,y)
c
      do 3 m=1,54
      n=iclr3(m)
      call box3(x3,y3,m,n)
    3 continue
      call frame3(x3,y3)
c
      write(1,*)' showpage'
      close(1)
c
c     Mix up the cube.
c
      write(*,*)' '
      write(*,*)'   Randomly mix up the puzzle. Enter an integer seed'
      write(*,*)'   for the random number generator. Enter 0 if you'
      write(*,*)'   dont want to mix up the cube.'
      write(*,*)' '
      read(*,*)iseed
c
      if(iseed.ne.0)then
      do 4 i=1,1000
      xran=1.0+11.999*ran0(iseed)
      n=ifix(xran)
      call opera(iclr,n,1)
      call opera3(iclr3,n,1)
    4 continue
c
      open(unit=1,status='unknown',file='cube53.ps')
      call start()
c
      do 5 m=1,150
      n=iclr(m)
      call box(x,y,m,n)
    5 continue
      call frame(x,y)
c
      do 6 m=1,54
      n=iclr3(m)
      call box3(x3,y3,m,n)
    6 continue
      call frame3(x3,y3)
      write(1,*)
      write(1,*)' showpage'
      close(1)
      end if
c
c     Solve the cube.
c
      write(*,*)' '
      write(*,*)' Enter an operator number one at a time.'
      write(*,*)' Definition of the operators is given in the display.'
      write(*,*)' Use negative numbers for inverse operators.'
      write(*,*)' Enter 99 to reset the cube. Enter 0 to quit.'
      write(*,*)' '
c
    7 read(*,*)i
      if(i.eq.0)go to 10
      k=1
      if(i.lt.0)then
      k=3
      i=-i
      end if
c
      call opera(iclr,i,k)
      call opera3(iclr3,i,k)
c
      open(unit=1,status='unknown',file='cube53.ps')
      call start()
c
      do 8 m=1,150
      n=iclr(m)
      call box(x,y,m,n)
    8 continue
      call frame(x,y)
      do 9 m=1,54
      n=iclr3(m)
      call box3(x3,y3,m,n)
    9 continue
      call frame3(x3,y3)
      write(1,*)
      write(1,*)' showpage'
      close(1)
      if(i.eq.99)go to 1
      go to 7
c
   10 stop
      end
c
      subroutine opera(iclr,i,kin)
      integer iclr(150)
      k=kin
      if(k.lt.0)k=3
c
      do 1 j=1,k
c
      if((i.eq.1).or.(i.eq.2).or.(i.eq.13))then
      call cycle(iclr,46,26,30,50)
      call cycle(iclr,41,27,35,49)
      call cycle(iclr,36,28,40,48)
      call cycle(iclr,31,29,45,47)
      call cycle(iclr,1,51,96,146)
      call cycle(iclr,6,56,91,141)
      call cycle(iclr,11,61,86,136)
      call cycle(iclr,16,66,81,131)
      call cycle(iclr,21,71,76,126)
      call cycle(iclr,42,32,34,44)
      call cycle(iclr,37,33,39,43)
      end if
c
      if((i.eq.2).or.(i.eq.13))then
      call cycle(iclr,2,52,97,147)
      call cycle(iclr,7,57,92,142)
      call cycle(iclr,12,62,87,137)
      call cycle(iclr,17,67,82,132)
      call cycle(iclr,22,72,77,127)
      end if
c
      if((i.eq.3).or.(i.eq.4))then
      call cycle(iclr,121,125,105,101)
      call cycle(iclr,122,120,104,106)
      call cycle(iclr,123,115,103,111)
      call cycle(iclr,124,110,102,116)
      call cycle(iclr,75,25,130,80)
      call cycle(iclr,70,20,135,85)
      call cycle(iclr,65,15,140,90)
      call cycle(iclr,60,10,145,95)
      call cycle(iclr,55,5,150,100)
      call cycle(iclr,117,119,109,107)
      call cycle(iclr,118,114,108,112)
      end if
c
      if(i.eq.4)then
      call cycle(iclr,74,24,129,79)
      call cycle(iclr,69,19,134,84)
      call cycle(iclr,64,14,139,89)
      call cycle(iclr,59,9,144,94)
      call cycle(iclr,54,4,149,99)
      end if
c
      if((i.eq.5).or.(i.eq.6).or.(i.eq.14))then
      call cycle(iclr,71,51,55,75)
      call cycle(iclr,72,56,54,70)
      call cycle(iclr,73,61,53,65)
      call cycle(iclr,74,66,52,60)
      call cycle(iclr,46,21,125,100)
      call cycle(iclr,47,22,124,99)
      call cycle(iclr,48,23,123,98)
      call cycle(iclr,49,24,122,97)
      call cycle(iclr,50,25,121,96)
      call cycle(iclr,67,57,59,69)
      call cycle(iclr,68,62,58,64)
      end if
c
      if((i.eq.6).or.(i.eq.14))then
      call cycle(iclr,41,16,120,95)
      call cycle(iclr,42,17,119,94)
      call cycle(iclr,43,18,118,93)
      call cycle(iclr,44,19,117,92)
      call cycle(iclr,45,20,116,91)
      end if
c
      if((i.eq.7).or.(i.eq.8))then
      call cycle(iclr,130,126,146,150)
      call cycle(iclr,129,131,147,145)
      call cycle(iclr,128,136,148,140)
      call cycle(iclr,127,141,149,135)
      call cycle(iclr,134,132,142,144)
      call cycle(iclr,133,137,143,139)
      call cycle(iclr,5,30,76,101)
      call cycle(iclr,4,29,77,102)
      call cycle(iclr,3,28,78,103)
      call cycle(iclr,2,27,79,104)
      call cycle(iclr,1,26,80,105)
      end if
c
      if(i.eq.8)then
      call cycle(iclr,10,35,81,106)
      call cycle(iclr,9,34,82,107)
      call cycle(iclr,8,33,83,108)
      call cycle(iclr,7,32,84,109)
      call cycle(iclr,6,31,85,110)
      end if
c
      if((i.eq.9).or.(i.eq.10).or.(i.eq.15))then
      call cycle(iclr,21,1,5,25)
      call cycle(iclr,16,2,10,24)
      call cycle(iclr,11,3,15,23)
      call cycle(iclr,6,4,20,22)
      call cycle(iclr,17,7,9,19)
      call cycle(iclr,18,12,8,14)
      call cycle(iclr,51,30,130,125)
      call cycle(iclr,52,35,129,120)
      call cycle(iclr,53,40,128,115)
      call cycle(iclr,54,45,127,110)
      call cycle(iclr,55,50,126,105)
      end if
c
      if((i.eq.10).or.(i.eq.15))then
      call cycle(iclr,56,29,135,124)
      call cycle(iclr,57,34,134,119)
      call cycle(iclr,58,39,133,114)
      call cycle(iclr,59,44,132,109)
      call cycle(iclr,60,49,131,104)
      end if
c
      if((i.eq.11).or.(i.eq.12))then
      call cycle(iclr,96,100,80,76)
      call cycle(iclr,97,95,79,81)
      call cycle(iclr,98,90,78,86)
      call cycle(iclr,99,85,77,91)
      call cycle(iclr,92,94,84,82)
      call cycle(iclr,93,89,83,87)
      call cycle(iclr,71,121,150,26)
      call cycle(iclr,72,116,149,31)
      call cycle(iclr,73,111,148,36)
      call cycle(iclr,74,106,147,41)
      call cycle(iclr,75,101,146,46)
      end if
c
      if(i.eq.12)then
      call cycle(iclr,66,122,145,27)
      call cycle(iclr,67,117,144,32)
      call cycle(iclr,68,112,143,37)
      call cycle(iclr,69,107,142,42)
      call cycle(iclr,70,102,141,47)
      end if
c
      if(i.eq.13)then
      call cycle(iclr,101,105,125,121)
      call cycle(iclr,106,104,120,122)
      call cycle(iclr,111,103,115,123)
      call cycle(iclr,116,102,110,124)
      call cycle(iclr,80,130,25,75)
      call cycle(iclr,85,135,20,70)
      call cycle(iclr,90,140,15,65)
      call cycle(iclr,95,145,10,60)
      call cycle(iclr,100,150,5,55)
      call cycle(iclr,107,109,119,117)
      call cycle(iclr,112,108,114,118)
      call cycle(iclr,79,129,24,74)
      call cycle(iclr,84,134,19,69)
      call cycle(iclr,89,139,14,64)
      call cycle(iclr,94,144,9,59)
      call cycle(iclr,99,149,4,54)
      call cycle(iclr,3,53,98,148)
      call cycle(iclr,8,58,93,143)
      call cycle(iclr,13,63,88,138)
      call cycle(iclr,18,68,83,133)
      call cycle(iclr,23,73,78,128)
      end if
c
      if(i.eq.14)then
      call cycle(iclr,150,146,126,130)
      call cycle(iclr,145,147,131,129)
      call cycle(iclr,140,148,136,128)
      call cycle(iclr,135,149,141,127)
      call cycle(iclr,144,142,132,134)
      call cycle(iclr,139,143,137,133)
      call cycle(iclr,101,76,30,5)
      call cycle(iclr,102,77,29,4)
      call cycle(iclr,103,78,28,3)
      call cycle(iclr,104,79,27,2)
      call cycle(iclr,105,80,26,1)
      call cycle(iclr,106,81,35,10)
      call cycle(iclr,107,82,34,9)
      call cycle(iclr,108,83,33,8)
      call cycle(iclr,109,84,32,7)
      call cycle(iclr,110,85,31,6)
      call cycle(iclr,11,115,90,36)
      call cycle(iclr,12,114,89,37)
      call cycle(iclr,13,113,88,38)
      call cycle(iclr,14,112,87,39)
      call cycle(iclr,15,111,86,40)
      end if
c
      if(i.eq.15)then
      call cycle(iclr,76,80,100,96)
      call cycle(iclr,81,79,95,97)
      call cycle(iclr,86,78,90,98)
      call cycle(iclr,91,77,85,99)
      call cycle(iclr,82,84,94,92)
      call cycle(iclr,87,83,89,93)
      call cycle(iclr,26,150,121,71)
      call cycle(iclr,31,149,116,72)
      call cycle(iclr,36,148,111,73)
      call cycle(iclr,41,147,106,74)
      call cycle(iclr,46,146,101,75)
      call cycle(iclr,27,145,122,66)
      call cycle(iclr,32,144,117,67)
      call cycle(iclr,37,143,112,68)
      call cycle(iclr,42,142,107,69)
      call cycle(iclr,47,141,102,70)
      call cycle(iclr,65,48,136,103)
      call cycle(iclr,64,43,137,108)
      call cycle(iclr,63,38,138,113)
      call cycle(iclr,62,33,139,118)
      call cycle(iclr,61,28,140,123)
      end if
c
    1 continue
c
      return
      end
c
      subroutine cycle(iclr,i1,i2,i3,i4)
      integer iclr(150)
c
      temp=iclr(i4)
      iclr(i4)=iclr(i3)
      iclr(i3)=iclr(i2)
      iclr(i2)=iclr(i1)
      iclr(i1)=temp
c
      return
      end
c
      subroutine box(x,y,m,n)
      real x(150,4),y(150,4)
c
      write(1,*)x(m,1),y(m,1),' moveto'
      write(1,*)x(m,2),y(m,2),' lineto'
      write(1,*)x(m,3),y(m,3),' lineto'
      write(1,*)x(m,4),y(m,4),' lineto'
      write(1,*)' closepath'
      call color(n)
      write(1,*)' fill'
c      
      return
      end
c
      subroutine frame(x,y)
      real x(150,4),y(150,4)
c
      do 1 m=1,150
      write(1,*)x(m,1),y(m,1),' moveto'
      write(1,*)x(m,2),y(m,2),' lineto'
      write(1,*)x(m,3),y(m,3),' lineto'
      write(1,*)x(m,4),y(m,4),' lineto'
      write(1,*)' closepath'
      call color(7)
      write(1,*)' stroke'
    1 continue
c      
      return
      end
c
      subroutine color(n)
c
      if(n.eq.1)write(1,*)' 0 0 1 setrgbcolor'
      if(n.eq.2)write(1,*)' 1 0.5 0 setrgbcolor'
      if(n.eq.3)write(1,*)' 1 1 0 setrgbcolor'
      if(n.eq.4)write(1,*)' 1 1 1 setrgbcolor'
      if(n.eq.5)write(1,*)' 1 0 0 setrgbcolor'
      if(n.eq.6)write(1,*)' 0 1 0 setrgbcolor'
      if(n.eq.7)write(1,*)' 0 0 0 setrgbcolor'
c
      return
      end
c
      subroutine setup(x,y,x0,y0,cube,cost,sint,iclr)
      real x(150,4),y(150,4)
      integer iclr(150)
c
      x(1,1)=x0-2.5*cube
      y(1,1)=y0-2.5*cube
      x(1,2)=x(1,1)+cube
      y(1,2)=y(1,1)
      x(1,3)=x(1,1)+cube
      y(1,3)=y(1,1)+cube
      x(1,4)=x(1,1)
      y(1,4)=y(1,1)+cube
c
      dx=cube
      dy=0.0
      delx=0.0
      dely=cube
      call shft(x,y,1,dx,dy,delx,dely)
c
      x(26,1)=x(1,1)-5.0*cost*cube
      y(26,1)=y(1,1)+5.0*sint*cube
      x(26,2)=x(26,1)+cost*cube
      y(26,2)=y(26,1)-sint*cube
      x(26,3)=x(26,2)
      y(26,3)=y(26,2)+cube
      x(26,4)=x(26,1)
      y(26,4)=y(26,1)+cube
c
      dx=cost*cube
      dy=-sint*cube
      delx=0.0
      dely=cube
      call shft(x,y,26,dx,dy,delx,dely)
c
      x(51,1)=x(21,4)
      y(51,1)=y(21,4)
      x(51,2)=x(51,1)+cube
      y(51,2)=y(51,1)
      x(51,3)=x(51,2)-cost*cube
      y(51,3)=y(51,2)+sint*cube
      x(51,4)=x(51,3)-cube
      y(51,4)=y(51,3)
c
      dx=cube
      dy=0.0
      delx=-cost*cube
      dely=sint*cube
      call shft(x,y,51,dx,dy,delx,dely)
c
      do 2 i=1,25
      do 1 j=1,4
      x(i+75,j)=x(i,j)-cost*8.0/3.0
      y(i+75,j)=y(i,j)+sint*8.0/3.0
    1 continue
    2 continue
c
      do 4 i=26,50
      do 3 j=1,4
      x(i+75,j)=x(i,j)+2.0
      y(i+75,j)=y(i,j)
    3 continue
    4 continue
c
      do 6 i=51,75
      do 5 j=1,4
      x(i+75,j)=x(i,j)
      y(i+75,j)=y(i,j)-2.0
    5 continue
    6 continue
c
      do 8 i=1,150
      do 7 j=1,4
      x(i,j)=x(i,j)*72.0
      y(i,j)=y(i,j)*72.0+50.0
    7 continue
    8 continue
c
      do 10 n=1,6
      i1=1+(n-1)*25
      do 9 i=i1,i1+24
      iclr(i)=n
    9 continue
   10 continue
c
      return
      end
c
      subroutine shft(x,y,n0,dx,dy,delx,dely)
      real x(150,4),y(150,4)
c
      n=n0
      do 2 j=1,5
      do 1 i=1,5
      x(n,1)=x(n0,1)+float(i-1)*dx+float(j-1)*delx
      y(n,1)=y(n0,1)+float(i-1)*dy+float(j-1)*dely
      x(n,2)=x(n0,2)+float(i-1)*dx+float(j-1)*delx
      y(n,2)=y(n0,2)+float(i-1)*dy+float(j-1)*dely
      x(n,3)=x(n0,3)+float(i-1)*dx+float(j-1)*delx
      y(n,3)=y(n0,3)+float(i-1)*dy+float(j-1)*dely
      x(n,4)=x(n0,4)+float(i-1)*dx+float(j-1)*delx
      y(n,4)=y(n0,4)+float(i-1)*dy+float(j-1)*dely
      n=n+1
    1 continue
    2 continue
c
      return
      end
c
      subroutine start()
c
      iz=100
      write(1,*)'%!'
      write(1,*)' /Times-Roman findfont'
      write(1,*)' 14 scalefont setfont'
      write(1,*)72, 560+iz,' moveto'
      write(1,*)' (Rotate a face:) show'
      write(1,*)72, 540+iz,' moveto'
      write(1,*)' ((outer layer, both layers) ) show'
      write(1,*)72, 520+iz,' moveto'
      write(1,*)' (  1, 2 = left) show'
      write(1,*)72, 500+iz,' moveto'
      write(1,*)' (  3, 4 = right) show'
      write(1,*)72, 480+iz,' moveto'
      write(1,*)' (  5, 6 = top) show'
      write(1,*)72, 460+iz,' moveto'
      write(1,*)' (  7, 8 = bottom) show'
      write(1,*)72, 440+iz,' moveto'
      write(1,*)' (  9, 10 = front) show'
      write(1,*)72, 420+iz,' moveto'
      write(1,*)' (  11, 12 = back) show'
      write(1,*)72, 380+iz,' moveto'
      write(1,*)' (Rotate the entire cube:) show'
      write(1,*)72, 360+iz,' moveto'
      write(1,*)' (  13 = left-right axis) show'
      write(1,*)72, 340+iz,' moveto'
      write(1,*)' (  14 = top-bottom axis) show'
      write(1,*)72, 320+iz,' moveto'
      write(1,*)' (  15 = front-back axis) show'
      write(1,*)72, 280+iz,' moveto'
      write(1,*)' (All rotations are counterclockwise.) show'
      write(1,*)300, 455,' moveto'
      write(1,*)' (Above: Exterior of the 5x5x5 cube.) show'
      write(1,*)300, 435,' moveto'
      write(1,*)' (Below: The 3x3x3 cube nested inside.) show'
c
      return
      end
c
      subroutine opera3(iclr,i,kin)
      integer iclr(54)
      k=kin
      if(k.lt.0)k=3
c
      do 1 j=1,k
      if(i.eq.2)then
      call cycle(iclr,52,1,19,34)
      call cycle(iclr,49,4,22,31)
      call cycle(iclr,46,7,25,28)
      call cycle(iclr,12,18,16,10)
      call cycle(iclr,11,15,17,13)
      end if
c
      if(i.eq.4)then
      call cycle(iclr,43,45,39,37)
      call cycle(iclr,44,42,38,40)
      call cycle(iclr,21,3,54,36)
      call cycle(iclr,24,6,51,33)
      call cycle(iclr,27,9,48,30)
      end if
c
      if(i.eq.6)then
      call cycle(iclr,19,21,27,25)
      call cycle(iclr,22,20,24,26)
      call cycle(iclr,9,43,34,18)
      call cycle(iclr,8,44,35,17)
      call cycle(iclr,7,45,36,16)
      end if
c
      if(i.eq.8)then
      call cycle(iclr,52,54,48,46)
      call cycle(iclr,53,51,47,49)
      call cycle(iclr,1,10,30,39)
      call cycle(iclr,2,11,29,38)
      call cycle(iclr,3,12,28,37)
      end if
c
      if(i.eq.10)then
      call cycle(iclr,7,1,3,9)
      call cycle(iclr,8,4,2,6)
      call cycle(iclr,18,46,39,21)
      call cycle(iclr,15,47,42,20)
      call cycle(iclr,12,48,45,19)
      end if
c
      if(i.eq.12)then
      call cycle(iclr,34,36,30,28)
      call cycle(iclr,35,33,29,31)
      call cycle(iclr,25,43,54,10)
      call cycle(iclr,26,40,53,13)
      call cycle(iclr,27,37,52,16)
      end if
c
      if(i.eq.13)then
      call cycle(iclr,52,1,19,34)
      call cycle(iclr,49,4,22,31)
      call cycle(iclr,46,7,25,28)
      call cycle(iclr,12,18,16,10)
      call cycle(iclr,11,15,17,13)
      call cycle(iclr,37,39,45,43)
      call cycle(iclr,40,38,42,44)
      call cycle(iclr,36,54,3,21)
      call cycle(iclr,33,51,6,24)
      call cycle(iclr,30,48,9,27)
      call cycle(iclr,26,29,47,8)
      call cycle(iclr,23,32,50,5)
      call cycle(iclr,20,35,53,2)
      end if
c
      if(i.eq.14)then
      call cycle(iclr,19,21,27,25)
      call cycle(iclr,22,20,24,26)
      call cycle(iclr,9,43,34,18)
      call cycle(iclr,8,44,35,17)
      call cycle(iclr,7,45,36,16)
      call cycle(iclr,46,48,54,52)
      call cycle(iclr,49,47,51,53)
      call cycle(iclr,39,30,10,1)
      call cycle(iclr,38,29,11,2)
      call cycle(iclr,37,28,12,3)
      call cycle(iclr,6,40,31,15)
      call cycle(iclr,5,41,32,14)
      call cycle(iclr,4,42,33,13)
      end if
c
      if(i.eq.15)then
      call cycle(iclr,7,1,3,9)
      call cycle(iclr,8,4,2,6)
      call cycle(iclr,18,46,39,21)
      call cycle(iclr,15,47,42,20)
      call cycle(iclr,12,48,45,19)
      call cycle(iclr,28,30,36,34)
      call cycle(iclr,31,29,33,35)
      call cycle(iclr,10,54,43,25)
      call cycle(iclr,13,53,40,26)
      call cycle(iclr,16,52,37,27)
      call cycle(iclr,24,17,49,38)
      call cycle(iclr,23,14,50,41)
      call cycle(iclr,22,11,51,44)
      end if
c
    1 continue
c
      return
      end
c
      subroutine box3(x,y,m,n)
      real x(54,4),y(54,4)
c
      write(1,*)x(m,1),y(m,1),' moveto'
      write(1,*)x(m,2),y(m,2),' lineto'
      write(1,*)x(m,3),y(m,3),' lineto'
      write(1,*)x(m,4),y(m,4),' lineto'
      write(1,*)' closepath'
      call color(n)
      write(1,*)' fill'
c      
      return
      end
c
      subroutine frame3(x,y)
      real x(54,4),y(54,4)
c
      do 1 m=1,54
      write(1,*)x(m,1),y(m,1),' moveto'
      write(1,*)x(m,2),y(m,2),' lineto'
      write(1,*)x(m,3),y(m,3),' lineto'
      write(1,*)x(m,4),y(m,4),' lineto'
      write(1,*)' closepath'
      call color(7)
      write(1,*)' stroke'
    1 continue
c      
      return
      end
c
      subroutine setup3(x,y,x0,y0,cube,cost,sint,iclr)
      real x(54,4),y(54,4)
      integer iclr(54)
c
      x(1,1)=x0-1.5*cube
      y(1,1)=y0-1.5*cube
      x(1,2)=x(1,1)+cube
      y(1,2)=y(1,1)
      x(1,3)=x(1,1)+cube
      y(1,3)=y(1,1)+cube
      x(1,4)=x(1,1)
      y(1,4)=y(1,1)+cube
c
      do 1 j=1,4
      x(2,j)=x(1,j)+cube
      y(2,j)=y(1,j)
      x(3,j)=x(1,j)+2.0*cube
      y(3,j)=y(1,j)
      x(4,j)=x(1,j)
      y(4,j)=y(1,j)+cube
      x(5,j)=x(1,j)+cube
      y(5,j)=y(1,j)+cube
      x(6,j)=x(1,j)+2.0*cube
      y(6,j)=y(1,j)+cube
      x(7,j)=x(1,j)
      y(7,j)=y(1,j)+2.0*cube
      x(8,j)=x(1,j)+cube
      y(8,j)=y(1,j)+2.0*cube
      x(9,j)=x(1,j)+2.0*cube
      y(9,j)=y(1,j)+2.0*cube
    1 continue
c
      x(12,1)=x(1,4)
      y(12,1)=y(1,4)
      x(12,2)=x(12,1)
      y(12,2)=y(12,1)-cube
      x(12,3)=x(12,2)-cost*cube
      y(12,3)=y(12,2)+sint*cube
      x(12,4)=x(12,3)
      y(12,4)=y(12,3)+cube
c
      do 2 j=1,4
      x(11,j)=x(12,j)-cost*cube
      y(11,j)=y(12,j)+sint*cube
      x(10,j)=x(11,j)-cost*cube
      y(10,j)=y(11,j)+sint*cube
      x(13,j)=x(10,j)
      y(13,j)=y(10,j)+cube
      x(14,j)=x(11,j)
      y(14,j)=y(11,j)+cube
      x(15,j)=x(12,j)
      y(15,j)=y(12,j)+cube
      x(16,j)=x(13,j)
      y(16,j)=y(13,j)+cube
      x(17,j)=x(14,j)
      y(17,j)=y(14,j)+cube
      x(18,j)=x(15,j)
      y(18,j)=y(15,j)+cube
    2 continue
c
      x(19,1)=x(7,4)
      y(19,1)=y(7,4)
      x(19,2)=x(19,1)+cube
      y(19,2)=y(19,1)
      x(19,3)=x(19,2)-cost*cube
      y(19,3)=y(19,2)+sint*cube
      x(19,4)=x(19,3)-cube
      y(19,4)=y(19,3)
c
      do 3 j=1,4
      x(20,j)=x(19,j)+cube
      y(20,j)=y(19,j)
      x(21,j)=x(20,j)+cube
      y(21,j)=y(20,j)
      x(22,j)=x(19,j)-cost*cube
      y(22,j)=y(19,j)+sint*cube
      x(23,j)=x(20,j)-cost*cube
      y(23,j)=y(20,j)+sint*cube
      x(24,j)=x(21,j)-cost*cube
      y(24,j)=y(21,j)+sint*cube
      x(25,j)=x(22,j)-cost*cube
      y(25,j)=y(22,j)+sint*cube
      x(26,j)=x(23,j)-cost*cube
      y(26,j)=y(23,j)+sint*cube
      x(27,j)=x(24,j)-cost*cube
      y(27,j)=y(24,j)+sint*cube
    3 continue
c
      do 5 i=1,9
      do 4 j=1,4
      x(i+27,j)=x(i,j)-cost*5.0/3.0
      y(i+27,j)=y(i,j)+sint*5.0/3.0
    4 continue
    5 continue
c
      do 7 i=10,18
      do 6 j=1,4
      x(i+27,j)=x(i,j)+1.25
      y(i+27,j)=y(i,j)
    6 continue
    7 continue
c
      do 9 i=19,27
      do 8 j=1,4
      x(i+27,j)=x(i,j)
      y(i+27,j)=y(i,j)-1.25
    8 continue
    9 continue
c
      do 11 i=1,54
      do 10 j=1,4
      x(i,j)=x(i,j)*72.0
      y(i,j)=y(i,j)*72.0
   10 continue
   11 continue
c
      do 13 n=1,6
      i1=1+(n-1)*9
      do 12 i=i1,i1+8
      iclr(i)=n
   12 continue
   13 continue
c
      return
      end
c
      function ran0(iseed)
      if(iseed.eq.0)iseed=178544878
      m=2147483647
      n=16807
      k=iseed/127773
      iseed=n*(iseed-k*127773)-k*2836
      if(iseed.lt.0)iseed=iseed+m
      ran0=float(iseed)/float(m)
      return
      end