program cube3
c
c     Mike Collins -- cinclodes@yahoo.com
c
c     This FORTRAN code simulates on the screen Rubik's Cube. There
c     are about 4.33x10^19 configurations of this classic puzzle.
c
c     The output file cube3.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(54,4),y(54,4)
      integer iclr(54)
    1 open(unit=1,status='unknown',file='cube3.ps')
      call start()
c
      x0=5.0
      y0=6.5
      cube=0.5
      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)
c
      do 2 m=1,54
      n=iclr(m)
      call box(x,y,m,n)
    2 continue
      call frame(x,y)
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 3 i=1,1000
      xran=1.0+5.999*ran0(iseed)
      n=ifix(xran)
      call opera(iclr,n,1)
    3 continue
c
      open(unit=1,status='unknown',file='cube3.ps')
      call start()
c
      do 4 m=1,54
      n=iclr(m)
      call box(x,y,m,n)
    4 continue
      call frame(x,y)
      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
    5 read(*,*)i
      if(i.eq.0)go to 7
      k=1
      if(i.lt.0)then
      k=3
      i=-i
      end if
c
      call opera(iclr,i,k)
c
      open(unit=1,status='unknown',file='cube3.ps')
      call start()
c
      do 6 m=1,54
      n=iclr(m)
      call box(x,y,m,n)
    6 continue
      call frame(x,y)
      write(1,*)
      write(1,*)' showpage'
      close(1)
      if(i.eq.99)go to 1
      go to 5
c
    7 stop
      end
c
      subroutine opera(iclr,i,kin)
      integer iclr(54)
      k=kin
      if(k.lt.0)k=3
c
      do 1 j=1,k
      if(i.eq.1)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.2)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.3)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.4)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.5)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.6)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.7)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.8)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.9)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 cycle(iclr,i1,i2,i3,i4)
      integer iclr(54)
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(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 frame(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 color(n)
c
      if(n.eq.1)write(1,*)' 0 0 1 setrgbcolor'
      if(n.eq.5)write(1,*)' 1 0 0 setrgbcolor'
      if(n.eq.3)write(1,*)' 1 1 0 setrgbcolor'
      if(n.eq.6)write(1,*)' 0 1 0 setrgbcolor'
      if(n.eq.2)write(1,*)' 1 0.5 0 setrgbcolor'
      if(n.eq.4)write(1,*)' 1 1 1 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(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*4.0
      y(i+27,j)=y(i,j)+sint*4.0
    4 continue
    5 continue
c
      do 7 i=10,18
      do 6 j=1,4
      x(i+27,j)=x(i,j)+3.0
      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)-3.0
    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
      subroutine start()
c
      write(1,*)'%!'
      write(1,*)' /Times-Roman findfont'
      write(1,*)' 14 scalefont setfont'
      write(1,*)' 72 500 moveto'
      write(1,*)' (Rotate a face:) show'
      write(1,*)' 72 480 moveto'
      write(1,*)' (  1 = left) show'
      write(1,*)' 72 460 moveto'
      write(1,*)' (  2 = right) show'
      write(1,*)' 72 440 moveto'
      write(1,*)' (  3 = top) show'
      write(1,*)' 72 420 moveto'
      write(1,*)' (  4 = bottom) show'
      write(1,*)' 72 400 moveto'
      write(1,*)' (  5 = front) show'
      write(1,*)' 72 380 moveto'
      write(1,*)' (  6 = back) show'
      write(1,*)' 72 350 moveto'
      write(1,*)' (Rotate the entire cube:) show'
      write(1,*)' 72 330 moveto'
      write(1,*)' (  7 = left-right axis) show'
      write(1,*)' 72 310 moveto'
      write(1,*)' (  8 = top-bottom axis) show'
      write(1,*)' 72 290 moveto'
      write(1,*)' (  9 = front-back axis) show'
      write(1,*)' 72 260 moveto'
      write(1,*)' (All rotations are counterclockwise.) show'
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