c***************************************************************************
c      mcubsub.f                                                           *
c                                                                          *
c          cube sub program                                                *    
c          main program   :   mcubvrm.f                                    *
c***************************************************************************
c-------------------------
       subroutine initvrml
c-------------------------
       write(10,10)
   10  format("#VRML V2.0 utf8")
       return
       end
c
c--------------------------------------
       subroutine pointlight(px,py,pz)
c--------------------------------------
       write(10,20)
   20  format("PointLight {")
       write(10,21)
   21  format("ambientIntensity 0.5")
       write(10,22)
   22  format("color 1.0 0.0 0.0")
       write(10,23)
   23  format("intensity 1.0")
       write(10,24) px,py,pz
   24  format("location ",f7.2,x,f7.2,x,f7.2)
       write(10,25)
   25  format("radius 50")
       write(10,29)
   29  format("}")
       return
       end
c-----------------------------------------
       subroutine directlight(ix,iy,iz)
c-----------------------------------------
       write(10,20)
   20  format("DirectionalLight {")
       write(10,21)
   21  format("ambientIntensity 0.5")
       write(10,22)
   22  format("color 1.0 0.0 0.0")
       write(10,23) ix,iy,iz
   23  format("direction ",i4,x,i4,x,i4)
       write(10,24)
   24  format("intensity 1.0")
       write(10,29)
   29  format("}")
       return
       end
c
c---------------------------------
       subroutine color1p(nx,ico)
c---------------------------------
       do 100 i=1,nx
c
       r=0.0
       g=0.0
       b=0.0
c
       if(ico.eq.1) then
         g=1.0
       end if
       if(ico.eq.2) then 
         r=1.0
         b=1.0
       end if 
       if(ico.eq.3) then 
         r=1.0
       end if
       if(ico.eq.4) then 
         r=1.0
         b=1.0
       end if
       if(ico.eq.5) then 
         b=1.0
       end if
       if(ico.eq.6) then
         r=1.0 
         g=1.0
       end if
       if(ico.eq.7) then 
         g=1.0
         b=1.0
       end if
       if(ico.eq.8) then
         r=1.0
         g=1.0
       end if
c
   10  continue
       if(i.ne.nx) write(10,20) r,g,b
       if(i.eq.nx) write(10,30) r,g,b
   20  format(f3.1,1x,f3.1,1x,f3.1,",")
   30  format(f3.1,1x,f3.1,1x,f3.1)
  100  continue
       return
       end
c
c-----------------------------
       subroutine coordp(nx,p)
c-----------------------------
       dimension p(1)
c
       do 100 i=1,nx
       i1=1+3*(i-1)
       i2=i1+1
       i3=i2+1
       x=p(i1)
       y=p(i2)
       z=p(i3)
       if(i.ne.nx) write(10,20) x,y,z
       if(i.eq.nx) write(10,30) x,y,z
   20  format(f7.3,1x,f7.3,1x,f7.3,",")
   30  format(f7.3,1x,f7.3,1x,f7.3)
  100  continue
       return
       end
c-----------------------
       subroutine facebe
c-----------------------
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("geometry IndexedFaceSet {")
       write(10,23)
   23  format("color Color {")
       write(10,24)
   24  format("color [")
       return
       end
c------------------------ 
       subroutine coordbe
c------------------------
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("coord Coordinate {")
       write(10,24)
   24  format("point [")
       return
       end
c-------------------------------
       subroutine indexp(nx,ib)
c-------------------------------
       dimension ib(1)
c
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex FALSE")
c
c----( colorIndexed )
       ibb=0
       do 30 i=1,nx
       if (ib(i).eq.-1) then
          ibb=ibb+1
       end if
   30  continue
       write(10,31)
   31  format("colorIndex [ ")
       write(10,32) (0,i=1,ibb+1)
   32  format(20i3)
       write(10,33)
   33  format("]")
c
c----( coordIndex )
       write(10,41)
   41  format("coordIndex [")
       write(10,42) (ib(i),i=1,nx)
   42  format(10i5)
c
       return
       end
c-----------------------
       subroutine faceen
c-----------------------
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("solid FALSE")
       write(10,23)
       write(10,23)
   23  format("}")
       return
       end
c
c----------------------------
       subroutine tranen
c       * transparency 0.5 *
c----------------------------
       write(10,20)
   20  format("]")
       write(10,21)
   21  format("solid FALSE",/"}")
       write(10,22)
   22  format("appearance Appearance {")
       write(10,23)
   23  format("material Material {")
       write(10,24)
   24  format("ambientIntensity 0.4")
       write(10,25)
   25  format("diffuseColor 0.4 0.4 0.0")
       write(10,26)
   26  format("transparency 0.5")
       write(10,27)
   27  format("emissiveColor 0.4 0.4 0.0")
c      write(10,28)
c  28  format("specularColor 0.1 0.1 0.1")
c      write(10,29)
c  29  format("shininess 0.3")
       write(10,30)
   30  format("}}}")
       return
       end
c
c----------------------------------------
       subroutine csub2(ic,ib,aa,p,pt,ct)
c----------------------------------------
       parameter(nol=12,nop=3*nol)
       dimension pp(nop),p(8),nf(nol,6),ii(10),jj(10)
       dimension ic(1),ib(1),aa(1),pt(1),ct(1)
c
       ico=ic(1)
       ua=aa(1)
       x1=aa(2)
       x2=aa(3)
       y1=aa(4)
       y2=aa(5)
       z1=aa(6)
       z2=aa(7)
c
       do 10 j1=1,6
       do 10 i1=1,nol
       nf(i1,j1)=0
   10  continue
c
       do 12 i1=1,nop
       pp(i1)=-1.0
   12  continue
c
c      u1=(ua-p(1))*(ua-p(2))
       u1=(p(1)-ua)*(p(2)-ua)
       if(u1.lt.0.0) then
       nf(1,1)=1
       nf(1,2)=1
c      x=(x1*(p(2)-ua)+x2*(ua-p(1)))/(p(2)-p(1))
       x=(x1-x2)/(p(1)-p(2))*(ua-p(1))+x1
       y=y1
       z=z1
       i1=1
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(2))*(ua-p(3))
       u1=(p(2)-ua)*(p(3)-ua)
       if(u1.lt.0.0) then
       nf(2,1)=1
       nf(2,3)=1
       x=x2
c      y=(y1*(p(3)-ua)+y2*(ua-p(2)))/(p(3)-p(2))
       y=(y1-y2)/(p(2)-p(3))*(ua-p(2))+y1
       z=z1
       i1=1+3*1
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(3))*(ua-p(4))
       u1=(p(3)-ua)*(p(4)-ua)
       if(u1.lt.0.0) then
       nf(3,1)=1
       nf(3,4)=1
c      x=(x2*(p(4)-ua)+x1*(ua-p(3)))/(p(4)-p(3))
       x=(x1-x2)/(p(4)-p(3))*(ua-p(4))+x1
       y=y2
       z=z1
       i1=1+3*2
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(4))*(ua-p(1))
       u1=(p(4)-ua)*(p(1)-ua)
       if(u1.lt.0.0) then
       nf(4,1)=1
       nf(4,5)=1
       x=x1
c      y=(y2*(p(1)-ua)+y1*(ua-p(4)))/(p(1)-p(4))
       y=(y1-y2)/(p(1)-p(4))*(ua-p(1))+y1
       z=z1
       i1=1+3*3
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c
c      u1=(ua-p(5))*(ua-p(6))
       u1=(p(5)-ua)*(p(6)-ua)
       if(u1.lt.0.0) then
       nf(5,2)=1
       nf(5,6)=1
c      x=(x1*(p(6)-ua)+x2*(ua-p(5)))/(p(6)-p(5))
       x=(x1-x2)/(p(5)-p(6))*(ua-p(5))+x1
       y=y1
       z=z2
       i1=1+3*4
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(6))*(ua-p(7))
       u1=(p(6)-ua)*(p(7)-ua)
       if(u1.lt.0.0) then
       nf(6,3)=1
       nf(6,6)=1
       x=x2
c      y=(y1*(p(7)-ua)+y2*(ua-p(6)))/(p(7)-p(6))
       y=(y1-y2)/(p(6)-p(7))*(ua-p(6))+y1
       z=z2
       i1=1+3*5
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(7))*(ua-p(8))
       u1=(p(7)-ua)*(p(8)-ua)
       if(u1.lt.0.0) then
       nf(7,4)=1
       nf(7,6)=1
c      x=(x2*(p(8)-ua)+x1*(ua-p(7)))/(p(8)-p(7))
       x=(x1-x2)/(p(8)-p(7))*(ua-p(8))+x1
       y=y2
       z=z2
       i1=1+3*6
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(8))*(ua-p(5))
       u1=(p(8)-ua)*(p(5)-ua)
       if(u1.lt.0.0) then
       nf(8,5)=1
       nf(8,6)=1
       x=x1
c      y=(y2*(p(5)-ua)+y1*(ua-p(8)))/(p(5)-p(8))
       y=(y1-y2)/(p(5)-p(8))*(ua-p(5))+y1
       z=z2
       i1=1+3*7
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c
c      u1=(ua-p(1))*(ua-p(5))
       u1=(p(1)-ua)*(p(5)-ua)
       if(u1.lt.0.0) then
       nf(9,2)=1
       nf(9,5)=1
       x=x1
       y=y1
c      z=(z1*(p(5)-ua)+z2*(ua-p(1)))/(p(5)-p(1))
       z=(z1-z2)/(p(1)-p(5))*(ua-p(1))+z1
       i1=1+3*8
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(2))*(ua-p(6))
       u1=(p(2)-ua)*(p(6)-ua)
       if(u1.lt.0.0) then
       nf(10,2)=1
       nf(10,3)=1
       x=x2
       y=y1
c      z=(z1*(p(6)-ua)+z2*(ua-p(2)))/(p(6)-p(2))
       z=(z1-z2)/(p(2)-p(6))*(ua-p(2))+z1
       i1=1+3*9
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(3))*(ua-p(7))
       u1=(p(3)-ua)*(p(7)-ua)
       if(u1.lt.0.0) then
       nf(11,3)=1
       nf(11,4)=1
       x=x2
       y=y2
c      z=(z1*(p(7)-ua)+z2*(ua-p(3)))/(p(7)-p(3))
       z=(z1-z2)/(p(3)-p(7))*(ua-p(3))+z1
       i1=1+3*10
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c      u1=(ua-p(4))*(ua-p(8))
       u1=(p(4)-ua)*(p(8)-ua)
       if(u1.lt.0.0) then
       nf(12,4)=1
       nf(12,5)=1
       x=x1
       y=y2
c      z=(z1*(p(8)-ua)+z2*(ua-p(4)))/(p(8)-p(4))
       z=(z1-z2)/(p(4)-p(8))*(ua-p(4))+z1
       i1=1+3*11
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c
       call csub1(ic,ib,p,pp,nf,pt,ct)
c
       return
       end
c
c-------------------------------------------       
       subroutine csub1(ic,ib,p,pp,nf,pt,ct)
c-------------------------------------------
       parameter(nol=12,nop=3*nol)
       dimension p(8),pp(nop),nf(nol,6),ii(10),jj(10),ic(1)
       dimension ib(1),pa(nop),pt(1),ct(1)
c
       ico=ic(1)
       iin=0
       do 10 i=1,4
       if(nf(i,1).eq.1) then
       iin=iin+1
       ii(iin)=i
       jj(iin)=1
       end if
   10  continue
c
       do 12 i=5,8
       if(nf(i,6).eq.1) then
       iin=iin+1
       ii(iin)=i
       jj(iin)=6
       end if
   12  continue
c
       if(iin.eq.0) go to 100
       if(iin.eq.2) go to 200
       if(iin.eq.4) go to 400
       go to 990
c
c-----
  100  continue
       do 110 i=9,12
       j=i-7
       if(nf(i,j).eq.1) then
       iin=iin+1
       ii(iin)=i
       jj(iin)=j
       end if
  110  continue
       go to 1000
c
c-----
  200  continue
       i1=ii(1)
       i2=ii(2)
       j1=jj(1)
       j2=jj(2)
c
       i3=0
       j3=0
       do 210 j=2,5
       if(nf(i2,j).eq.1) then
       j3=j
       end if
  210  continue
       do 212 i=9,12
       if(i.eq.i2) go to 212
       if(nf(i,j3).eq.1) then
       iin=iin+1
       i3=i
       ii(iin)=i3
       jj(iin)=j3
       end if
  212  continue
c
       i4=0
       do 220 j=2,5
       if(j.eq.j3) go to 220
       if(nf(i3,j).eq.1) then
       j4=j
       end if
  220  continue
       do 222 i=9,12
       if(i.eq.i2) go to 222
       if(i.eq.i3) go to 222
       if(nf(i,j4).eq.1) then
       iin=iin+1
       i4=i
       ii(iin)=i4
       jj(iin)=j4
       end if
  222  continue
       if(i4.eq.0) go to 1000
c
       i5=0
       j5=0
       do 230 j=2,5
       if(j.eq.j3) go to 230
       if(j.eq.j4) go to 230
       if(nf(i4,j).eq.1) then
       j5=j
       end if
  230  continue
       do 232 i=9,12
       if(i.eq.i2) go to 232
       if(i.eq.i3) go to 232
       if(i.eq.i4) go to 232
       if(nf(i,j5).eq.1) then
       iin=iin+1
       i5=i
       ii(iin)=i5
       jj(iin)=j5
       end if
  232  continue
       go to 1000
c
c-----
  400  continue
       iin=2
       i1=ii(1)
       i2=ii(2)
       j1=jj(1)
       j2=jj(2)
c
       do 410 j=2,5
       if(nf(i2,j).eq.1) then
       j3=j
       end if
  410  continue
       do 412 i=5,12
       if(i.eq.i2) go to 412
       if(nf(i,j3).eq.1) then
       iin=iin+1
       i3=i
       ii(iin)=i3
       jj(iin)=j3
       end if
  412  continue
c
       i4=0
       j4=0
       do 420 j=2,6
       if(j.eq.j3) go to 420
       if(nf(i3,j).eq.1) then
       j4=j
       end if
  420  continue
       do 422 i=5,12
       if(i.eq.i2) go to 422
       if(i.eq.i3) go to 422
       if(nf(i,j4).eq.1) then
       iin=iin+1
       i4=i
       ii(iin)=i4
       jj(iin)=j4
       end if
  422  continue
       if(i4.eq.0) go to 1000
c
       i5=0
       j5=0
       do 430 j=2,6
       if(j.eq.j3) go to 430
       if(j.eq.j4) go to 430
       if(nf(i4,j).eq.1) then
       j5=j
       end if
  430  continue
       do 432 i=5,12
       if(i.eq.i2) go to 432
       if(i.eq.i3) go to 432
       if(i.eq.i4) go to 432
       if(nf(i,j5).eq.1) then
       iin=iin+1
       i5=i
       ii(iin)=i5
       jj(iin)=j5
       end if
  432  continue
       if(i5.eq.0) go to 1000
c
       i6=0
       j6=0
       do 440 j=2,6
       if(j.eq.j3) go to 440
       if(j.eq.j4) go to 440
       if(j.eq.j5) go to 440
       if(nf(i5,j).eq.1) then
       j6=j
       end if
  440  continue
       do 442 i=5,12
       if(i.eq.i2) go to 442
       if(i.eq.i3) go to 442
       if(i.eq.i4) go to 442
       if(i.eq.i5) go to 442
       if(nf(i,j6).eq.1) then
       iin=iin+1
       i6=i
       ii(iin)=i6
       jj(iin)=j6
       end if
  442  continue
c
c-----
 1000  continue
       do 1010 j=1,iin
       ii1=ii(j)
       jj1=jj(j)
       i1=1+3*(ii1-1)
       i2=i1+1
       i3=i2+1
       j1=1+3*(j-1)
       j2=j1+1
       j3=j2+1
c
       ic(5)=ic(5)+1
       ic(6)=ic(6)+1
       k=ic(5)
       ka=ic(6)
       k1=1+3*(k-1)
       k2=k1+1
       k3=k2+1
       ct(k)=p(j)
       ib(ka)=k-1
c
       x=pp(i1)
       y=pp(i2)
       z=pp(i3)
       pt(k1)=x
       pt(k2)=y
       pt(k3)=z
 1010  continue
c
       if(iin.eq.0) go to 990
       ic(6)=ic(6)+1
       ka=ic(6)
       ib(ka)=-1
c
c
  990  continue
       return
       end
c--------------------------------------------
c      cube subroutime program
       subroutine csuba(ic,aa,u,mxl,myl,mzl)
c--------------------------------------------
       parameter(ipt=99999,ipp=ipt*3)
       parameter(nol=12,nop=3*nol)
       dimension pp(nop),p(8),nf(nol,6)
       dimension ic(1),aa(1),u(mxl,myl,mzl)
       dimension ib(ipt),ct(ipt),pt(ipp)
c
       ico=ic(1)
       nx=ic(2)
       ny=ic(3)
       nz=ic(4)
       ua=aa(1)
       x1=aa(2)
       x2=aa(3)
       y1=aa(4)
       y2=aa(5)
       z1=aa(6)
       z2=aa(7)
       xl=aa(8)
       yl=aa(9)
       zl=aa(10)
       x0=aa(11)
       y0=aa(12)
       z0=aa(13)
c
       n1=nx
       n2=n1*ny
       n3=n2*nz
       nxs=nx-1
       nys=ny-1
       nzs=nz-1
       dx=xl/float(nxs)
       dy=yl/float(nys)
       dz=zl/float(nzs)
c###
       x0=-1*dx*(nx/2)/2
       y0=-1*dy*ny/2
       z0=-1*dz*nz/2
c###
       do 200 k=1,nzs
       z1=z0+dz*float(k-1)
       z2=z1+dz
       ic(5)=0
       ic(6)=0
c
       do 100 j=1,nys
       y1=y0+dy*float(j-1)
       y2=y1+dy
       do 100 i=1,nxs
       x1=x0+dx*float(i-1)
       x2=x1+dx
c
       p(1)=u(i,j,k)
       p(2)=u(i+1,j,k)
       p(3)=u(i+1,j+1,k)
       p(4)=u(i,j+1,k)
       p(5)=u(i,j,k+1)
       p(6)=u(i+1,j,k+1)
       p(7)=u(i+1,j+1,k+1)
       p(8)=u(i,j+1,k+1)
c
       aa(1)=ua
       aa(2)=x1
       aa(3)=x2
       aa(4)=y1
       aa(5)=y2
       aa(6)=z1
       aa(7)=z2
c
       call csub2(ic,ib,aa,p,pt,ct)
c
  100  continue
c
       iin1=ic(5)
       iin2=ic(6)
c
       if(iin1.eq.0) go to 200
c
       call facebe
       call color1p(1,ico)
       call coordbe
       call coordp(iin1,pt)
       iin21=iin2-1
       call indexp(iin21,ib)
       call tranen
c
  200  continue
       return
       end
c
c      ------------------------
       subroutine viewp(iviewp)
c      ------------------------
c      for VRML plug-in diplay of personal computer
c
       write(10,10)
   10  format("Viewpoint {")
       write(10,20)
   20  format("orientation 0 0 -1 0")
       write(10,30) iviewp
   30  format("position 0 0 ",i4)
       write(10,40)
   40  format("}")
       return
       end
c
c      -------------------------
       subroutine backgr1(r,g,b)
c      -------------------------
       write(10,10)
   10  format("Background {", "groundColor [")
       write(10,12) r,g,b
   12  format(f3.1,2(1x,f3.1),",")
       write(10,12) r,g,b
       write(10,14) r,g,b
   14  format(f3.1,2(1x,f3.1))
       write(10,16)
   16  format("]")
       write(10,18)
   18  format("groundAngle [1.05, 1.57]")
       write(10,20)
   20  format("skyColor [")
       write(10,12) r,g,b
       write(10,12) r,g,b
       write(10,14) r,g,b
       write(10,16)
       write(10,22)
   22  format("skyAngle [1.05, 1.57]")
       write(10,24)
   24  format("}")
       return
       end
