c************************************************************
c           subroutine file     << sub2a.f >>               *
c                                                           *
c      update ; 01/09/12     /home/g3/ogino/VRML/jst2k1     *
c************************************************************
c
c      ---------------------
c      subroutine initvrml
c      ---------------------
c      initialization of VRML
       subroutine initvrml
       write(10,10)
   10  format("#VRML V2.0 utf8")
       return
       end
c      ---------------------
c      subroutine point3d
c      ---------------------
       subroutine point3d(nx,ny,nz,x0,y0,z0,xl,yl,zl,ico,u)
c
c      nx: number in x-direction
c      ny: number in y-direction
c      nz: number in z-direction
c      x0: x-value at the left and bottom point
c      y0: y-value at the left and bottom point
c      z0: z-value at the left and bottom point
c      xl: length in x-direction
c      yl: length in y-direction
c      zl: length in z-direction
c      ico: color code
c      u(nx*ny): input data [0.0,1.0]
c
       dimension u(1)
c
       hx=xl/float(nx-1)
       hy=yl/float(ny-1)
       hz=zl/float(nz-1)
c
       call colorb
       call colorpb(ico,nx,ny,nz,u)
       call colore
c
       call pointb
       call pointp(nx,ny,nz,x0,y0,z0,hx,hy,hz)
       call pointe
c
       return
       end
c      ---------------------
c      subroutine colorb
c      ---------------------
       subroutine colorb
c
       write(10,*) " Shape {"
       write(10,*) " geometry PointSet {"
       write(10,*) " color Color {"
       write(10,*) " color ["
c
       return
       end
c      ---------------------
c      subroutine colore
c      ---------------------
       subroutine colore
       write(10,*) " ]"
       write(10,*) " }"
c
       return
       end
c      ---------------------
c      subroutine pointb
c      ---------------------
       subroutine pointb
       write(10,*) " coord Coordinate {"
       write(10,*) " point ["
       return
       end
       subroutine pointe
       write(10,*) " ]"
       write(10,*) " }"
       write(10,*) " }"
       write(10,*) " }"
c
       return
       end
c      ---------------------
c      subroutine colorpb
c      ---------------------
       subroutine colorpb(ico,nx,ny,nz,u)
       dimension u(1)
c
       n1=nx
       n2=n1*ny
       n3=n2*nz
       last=nx*ny*nz
c
       do 10 k=1,nz
       do 10 j=1,ny
       do 10 i=1,nx
       x=hx*float(i-1)
       y=hy*float(j-1)
       z=hz*float(k-1)
       ii=i+n1*(j-1)+n2*(k-1)
c
       x1=u(ii)
c
       a0=1.20
       a0=1.00
       if(ico.eq.2) go to 32
       if(ico.eq.3) go to 33
       if(ico.eq.4) go to 34
       if(ico.eq.5) go to 35
       if(ico.eq.6) go to 36
       if(ico.eq.7) go to 37
       r=a0-(a0-0.5)*x1
       g=a0-(a0-0.5)*x1
       b=a0-(a0-0.5)*x1
       go to 40
   32  continue
       r=a0-a0*x1
       g=1.0
       b=a0-a0*x1
       go to 40
   33  continue
       r=1.0
       g=a0-a0*x1
       b=a0-a0*x1
       go to 40
   34  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       g=a0+3.0-y1
       b=a0+1.0-y1
       if(y1.ge.a0+2.0) b=y1-(a0+3.0)
       go to 40
   35  continue
       a0=2.0
       y1=x1*(a0+4.0)
       b=a0-y1
       if(y1.ge.a0) b=y1-(a0+1.0)
       r=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   36  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       b=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   37  continue
       x1=10.0*x1
       b=x1
       if(x1.ge.3.0) b=4.0-x1
       if(x1.ge.6.0) b=0.25*(x1-6.0)
       g=x1-2.0
       if(x1.ge.5.0) g=6.0-x1
       if(x1.ge.6.0) g=0.25*(x1-6.0)
       r=x1
       if(x1.ge.1.0) r=2.0-x1
       if(x1.ge.4.0) r=x1-4.0
       go to 40
c
   40  continue
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       if(ii.lt.last) then
       write(10,20) r,g,b
   20  format(1h ,f3.1,1x,f3.1,1x,f3.1,1h,)
       else
       write(10,22) r,g,b
   22  format(1h ,f3.1,1x,f3.1,1x,f3.1)
       end if
   10  continue
c
       return
       end
c      ---------------------
c      subroutine pointp
c      ---------------------
       subroutine pointp(nx,ny,nz,x0,y0,z0,hx,hy,hz)
c
       n1=nx
       n2=n1*ny
       n3=n2*nz
       last=nx*ny*nz
c
       do 10 k=1,nz
       z=z0+hz*float(k-1)
       do 10 j=1,ny
       y=y0+hy*float(j-1)
       do 10 i=1,nx
       x=x0+hx*float(i-1)
       ii=i+n1*(j-1)+n2*(k-1)
       if(ii.lt.last) then
       write(10,20) x,z,y
   20  format(1h ,f5.2,1x,f5.2,1x,f5.2,1h,)
       else
       write(10,22) x,z,y
   22  format(1h ,f5.2,1x,f5.2,1x,f5.2)
       end if
   10  continue
       return
       end
c      ---------------------
c      subroutine line3d
c      ---------------------
       subroutine line3d(lasl,r,g,b,po)
c
c      lasl: number of points on line
c      r: red color   [0.0,1.0]
c      g: green color [0.0,1.0]
c      b: blue color  [0.0,1.0]
c      po(3*lasl): x,y,z-values in position
c
       dimension po(1)
c
       call linebe(r,g,b)
c
       do 40 ii=1,lasl
       lin=ii
       li1=3*(lin-1)+1
       li2=3*(lin-1)+2
       li3=3*lin
       x=po(li1)
       y=po(li2)
       z=po(li3)
       call linep(x,y,z,ii,lasl)
   40  continue
c
       call lineen(lasl)
       return
       end
c      ---------------------
c      subroutine symbvb
c      ---------------------
       subroutine symbvb
       write(10,10)
   10  format("Viewpoint { position 10 4 8 }")
       return
       end
c      ---------------------
c      subroutine symblv
c      ---------------------
       subroutine symblv(x,y,z,h,chr,n)
c
c      x: x-position of start point of character
c      y: y-position of start point of character
c      z: z-position of start point of character
c      h: height of character
c      chr: character
c      n: number of character
c
       character chr*80,ica*80,ich(80)*1
       equivalence (ica,ich)
       ica=chr
       write(10,10)
c  10  format("Viewpoint { position 10 4 8 }"/"Transform {")
   10  format("Transform {")
       write(10,12) x,y,z
   12  format("translation ",3(1x,f5.1))
       write(10,14)
   14  format("children ["/"Shape {"/"geometry Text {")  
c      write(10,16) (ich(i),i=1,n)
c  16  format("string ",'"',80a1,'"')
       write(10,*) "string ",'"',(ich(i),i=1,n),'"'
       write(10,18)
   18  format("fontStyle FontStyle {")
       write(10,20) h
   20  format("size ",f5.1)      
       write(10,22)
   22  format('family "SERIF"'/'style "ROMAN"')
       write(10,24)
   24  format("}"/"}"/"}"/"]"/"}")
       return
       end
c      ---------------------
c       subroutine symblvc
c      ---------------------
       subroutine symblvc(x,y,z,h,r,g,b,chr,n)
c
c      x: x-position of start point of character
c      y: y-position of start point of character
c      z: z-position of start point of character
c      h: height of character
c      r: red   [0.0,1.0]
c      g: green [0.0,1.0]
c      b: blue  [0.0,1.0]
c      chr: character
c      n: number of character
c
       character chr*80,ica*80,ich(80)*1
       equivalence (ica,ich)
       ica=chr
       write(10,10)
c  10  format("Viewpoint { position 10 4 8 }"/"Transform {")
   10  format("Transform {")
       write(10,12) x,y,z
   12  format("translation ",3(1x,f5.1))
       write(10,14) r,g,b
   14  format("children ["/"Shape {"/"appearance Appearance {"
     1  /"material Material {"/"diffuseColor ",3(1x,f3.1)/
     2  "}"/"}"/"geometry Text {")
c      write(10,16) (ich(i),i=1,n)
c  16  format("string ",'"',80a1,'"')
       write(10,*) "string ",'"',(ich(i),i=1,n),'"'
       write(10,18)
   18  format("fontStyle FontStyle {")
       write(10,20) h
   20  format("size ",f5.1)
       write(10,22)
   22  format('family "SERIF"'/'style "ROMAN"')
       write(10,24)
   24  format("}"/"}"/"}"/"]"/"}")
       return
       end
c      ---------------------
c      subroutine image1
c      ---------------------
       subroutine image1(nx,ny,xb,yb,xl,yl,ico,icc,zcc,u,p)
c
c      nx: number in x-direction
c      nx: number in x-direction
c      x0: x-value at the left and bottom point
c      y0: y-value at the left and bottom point
c      xl: length in x-direction
c      yl: length in y-direction
c      ico: color code
c      icc: selection of plane, i
c         : 1 for xy plane, s2 for xz plane, 3 for yz plane
c      zcc: position of selected plane
c      u(nx*ny): input data
c      p: work array
c
       dimension p(1),u(1)
       ny0=ny-1
       nx2=nx*2
c
       do 100 ii=1,ny0
c
       call trans2(nx,ny,ii,u,p)
       call facebe
c
       call colorp(nx2,p,ico)
c
       call coordbe
c
c      call coordp(nx2,p)
       call coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p)
c
       call coorden(nx)
       call faceen
c
  100  continue
c
       return
       end
c      ---------------------
c      subroutine image2
c      ---------------------
       subroutine image2(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p)
c
c      nx: number in x-direction
c      nx: number in x-direction
c      x0: x-value at the left and bottom point
c      y0: y-value at the left and bottom point
c      xl: length in x-direction
c      yl: length in y-direction
c      ico: color code
c      icc: selection of plane, i
c         : 1 for xy plane, s2 for xz plane, 3 for yz plane
c      zcc: position of selected plane
c      v: (nx-1)*(ny-1) additional input data, unnecessary
c      u(nx*ny): input data
c      p: work array
c
       dimension p(1),u(1),v(1)
       ny0=ny-1
       nx2=nx*2
c
       do 100 ii=1,ny0
c
       call trans2(nx,ny,ii,u,p)
       call facebe
c
       call colorp(nx2,p,ico)
c
       call coordbe
c
c      call coordp(nx2,p)
       call coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p)
c
       call coorden(nx)
       call faceen
c
  100  continue
c
       return
       end
c      ---------------------
c      subroutine trans1
c      ---------------------
       subroutine trans1(nx,ny,v,u)
       dimension v(1),u(1)
c
       nx0=nx-1
       ny0=ny-1
       nxa=nx0-1
       nya=ny0-1
c
       do 10 j=1,nya
       do 10 i=1,nxa
       i1=i+1+nx*j
       j1=i+nx0*(j-1)
       u(i1)=0.25*(v(j1)+v(j1+1)+v(j1+nx0)+v(j1+nx))
   10  continue
c
       do 20 i=1,nxa
       i1=i+1
       i2=i1+nx*ny0
       j1=i+1
       j2=i+1+nx0*nya
       u(i1)=0.75*(v(j1-1)+v(j1))-0.25*(v(j1-1+nx0)+v(j1+nx0))
       u(i2)=0.75*(v(j2-1)+v(j2))-0.25*(v(j2-1-nx0)+v(j2-nx0))
   20  continue
c
       do 30 j=1,nya
       i1=1+nx*j
       i2=i1+nx0
       j1=1+nx0*j
       j2=j1+nxa
       u(i1)=0.75*(v(j1-nx0)+v(j1))-0.25*(v(j1-nx0+1)+v(j1+1))
       u(i2)=0.75*(v(j2-nx0)+v(j2))-0.25*(v(j2-nx0-1)+v(j2-1))
   30  continue
c
       i1=1
       i2=i1+nx0
       i3=1+nx*ny0
       i4=i3+nx0
       j1=1
       j2=j1+nxa
       j3=1+nx0*nya
       j4=j3+nxa
       u(i1)=0.5*(3.0*v(j1)-v(j1+nx))
       u(i2)=0.5*(3.0*v(j2)-v(j2+nxa))
       u(i3)=0.5*(3.0*v(j3)-v(j3-nxa))
       u(i4)=0.5*(3.0*v(j4)-v(j4-nx))
c
c      write(6,*) i1,i2,i3,i4
c      write(6,*) j1,j2,j3,j4
c      write(6,*) nx,nx0,nxa,ny,ny0,nya
c
       return
       end
c      ---------------------
c      subroutine trans2
c      ---------------------
       subroutine trans2(nx,ny,ii,u,p)
       dimension p(1),u(1)
       do 10 i=1,nx
       j1=i+nx*(ii-1)
       j2=j1+nx
       i1=1+2*(i-1)
       i2=i1+1
       p(i1)=u(j1)
       p(i2)=u(j2)
   10  continue
       return
       end
c      ---------------------
c      subroutine colorp
c      ---------------------
       subroutine colorp(nx,p,ico)
       dimension p(1)
       r=1.0
       g=1.0
       b=1.0
       do 100 i=1,nx
       x1=p(i)
c
       a0=1.20
       a0=1.00
       if(ico.eq.2) go to 32
       if(ico.eq.3) go to 33
       if(ico.eq.4) go to 34
       if(ico.eq.5) go to 35
       if(ico.eq.6) go to 36
       if(ico.eq.7) go to 37
       r=a0-(a0-0.5)*x1
       g=a0-(a0-0.5)*x1
       b=a0-(a0-0.5)*x1
       go to 40
   32  continue
       r=a0-a0*x1
       g=1.0
       b=a0-a0*x1
       go to 40
   33  continue
       r=1.0
       g=a0-a0*x1
       b=a0-a0*x1
       go to 40
   34  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       g=a0+3.0-y1
       b=a0+1.0-y1
       if(y1.ge.a0+2.0) b=y1-(a0+3.0)
       go to 40
   35  continue
       a0=2.0
       y1=x1*(a0+4.0)
       b=a0-y1
       if(y1.ge.a0) b=y1-(a0+1.0)
       r=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   36  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       b=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   37  continue
       x1=10.0*x1
       b=x1
       if(x1.ge.3.0) b=4.0-x1
       if(x1.ge.6.0) b=0.25*(x1-6.0)
       g=x1-2.0
       if(x1.ge.5.0) g=6.0-x1
       if(x1.ge.6.0) g=0.25*(x1-6.0)
       r=x1
       if(x1.ge.1.0) r=2.0-x1
       if(x1.ge.4.0) r=x1-4.0
       go to 40
c
   40  continue
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       if(i.ne.nx) write(10,20) r,g,b
       if(i.eq.nx) write(10,30) r,g,b
   20  format(f4.2,1x,f4.2,1x,f4.2,",")
   30  format(f4.2,1x,f4.2,1x,f4.2)
  100  continue
       return
       end
c      ---------------------
c      subroutine coordp
c      ---------------------
       subroutine coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p)
       dimension p(1)
       x=1.0
       y=1.0
       z=1.0
c      xb=0.0
c      yb=0.0
c      dx=1.0
c      dy=1.0
       nxy=nx*ny
       dx=xl/float(nx-1)
       dy=yl/float(ny-1)
c      do 100 j1=1,2
       do 100 i=1,nx
       do 100 j1=1,2
       j=j1-1+ii
       ixy=i+nx*(j-1)
c
       if (icc.eq.2) go to 52
       if (icc.eq.3) go to 53
c
       x=xb+dx*float(i-1)
       y=yb+dy*float(j-1)
       z=zcc
       go to 60
   52  continue
       y=xb+dx*float(i-1)
       z=yb+dy*float(j-1)
       x=zcc
       go to 60
   53  continue
       x=xb+dx*float(i-1)
       z=yb+dy*float(j-1)
       y=zcc
   60  continue
c
       if(i.ne.nx) write(10,20) x,y,z
       if(i.eq.nx) write(10,30) x,y,z
   20  format(f4.2,1x,f4.2,1x,f4.2,",")
   30  format(f4.2,1x,f4.2,1x,f4.2)
  100  continue
       return
       end
c      ---------------------
c      subroutine facebe
c      ---------------------
       subroutine facebe
       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      ---------------------
c      subroutine coordbe
c      ---------------------
       subroutine coordbe
       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      ---------------------
c      subroutine coorden
c      ---------------------
       subroutine coorden(nx)
       dimension ia(2000)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex TRUE")
       write(10,24)
   24  format("coordIndex [")
c
       do 10 i=1,nx
       i1=1+5*(i-1)
       i2=2+5*(i-1)
       i3=3+5*(i-1)
       i4=4+5*(i-1)
       i5=5+5*(i-1)
       ia(i1)=0+2*(i-1)
       ia(i2)=1+2*(i-1)
       ia(i3)=3+2*(i-1)
       ia(i4)=2+2*(i-1)
       ia(i5)=-1
   10  continue
c
       na2=5*(nx-1)-1
       write(10,102) (ia(i),i=1,na2)
  102  format(4(4i2,i3),2i2,3i3,220i3,2i3,2i4,i3,
     1        449(4i4,i3),2i4,2i5,i3,3000(4i5,i3))
c
       return
       end
c      ---------------------
c      subroutine faceen
c      ---------------------
       subroutine faceen
       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 linep
c      ---------------------
       subroutine linep(x,y,z,ii,lasl)
c
       if(ii.lt.lasl) then
       write(10,10) x,z,y
c  10  format(f5.1,1x,f5.1,1x,f5.1,1h,)
c  10  format(f5.2,1x,f5.2,1x,f5.2,1h,)
   10  format(f6.2,1x,f6.2,1x,f6.2,1h,)
       else
       write(10,20) x,z,y
c  20  format(f5.1,1x,f5.1,1x,f5.1)
c  20  format(f5.2,1x,f5.2,1x,f5.2)
   20  format(f6.2,1x,f6.2,1x,f6.2)
       end if
       return
       end
c      ---------------------
c      subroutine linep2
c      ---------------------
       subroutine linep2(x,y,z,ii,lasl)
c
       if(ii.lt.lasl) then
       write(10,10) x,z,y
   10  format(f5.1,1x,f5.1,1x,f5.1,1h,)
       else
       write(10,20) x,z,y
   20  format(f5.1,1x,f5.1,1x,f5.1)
       end if
       return
       end
c      ---------------------
c      subroutine linebe
c      ---------------------
       subroutine linebe(r,g,b)
c
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("geometry IndexedLineSet {")
       write(10,23)
   23  format("color Color {")
       write(10,20) r,g,b
   20  format("color [ ",f3.1,1x,f3.1,1x,f3.1," ]")
       write(10,25)
   25  format("}")
       write(10,26)
   26  format("coord Coordinate {")
       write(10,27)
   27  format("point [")
       return
       end
c      ---------------------
c      subroutine lineen
c      ---------------------
       subroutine lineen(lasl)
       integer ic(9000)
c
       do 10 i=1,9000
       ic(i)=i-1
   10  continue
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorIndex [ 0 ]")
       write(10,24)
   24  format("coordIndex [")
c
       write(10,20) (ic(i),i=1,lasl)
   20  format(10i2,90i3,900i4/(2000i5))
c
       write(10,25)
   25  format("]")
       write(10,26)
   26  format("colorPerVertex FALSE")
       write(10,27)
   27  format("}")
       write(10,28)
   28  format("}")
       return
       end
c      ---------------------
c      subroutine rect1
c      ---------------------
       subroutine rect1(icc,zcc,r,g,b,x1,y1,x2,y2)
       lasl=5
       call linebe(r,g,b)
c
       if (icc.eq.2) go to 52
       if (icc.eq.3) go to 53
c
       z=zcc
       call linep(x1,y1,z,1,lasl)
       call linep(x2,y1,z,2,lasl)
       call linep(x2,y2,z,3,lasl)
       call linep(x1,y2,z,4,lasl)
       call linep(x1,y1,z,5,lasl)
       go to 60
c
   52  continue
       x=zcc
       call linep(x,x1,y1,1,lasl)
       call linep(x,x2,y1,2,lasl)
       call linep(x,x2,y2,3,lasl)
       call linep(x,x1,y2,4,lasl)
       call linep(x,x1,y1,5,lasl)
       go to 60
c
   53  continue
       y=zcc
       call linep(x1,y,y1,1,lasl)
       call linep(x2,y,y1,2,lasl)
       call linep(x2,y,y2,3,lasl)
       call linep(x1,y,y2,4,lasl)
       call linep(x1,y,y1,5,lasl)
   60  continue
       call lineen(lasl)
c
       return
       end       
c      zcsub1.f subroutine program
c      ---------------------
c      subroutine zcoordp
c      ---------------------
       subroutine zcoordp(nx,p)
       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
c  20  format(f3.1,1x,f3.1,1x,f3.1,",")
c  30  format(f3.1,1x,f3.1,1x,f3.1)
   20  format(f5.3,1x,f5.3,1x,f5.3,",")
   30  format(f5.3,1x,f5.3,1x,f5.3)
  100  continue
       return
       end
c      ---------------------
c      subroutine zcoorden
c      ---------------------
       subroutine zcoorden(nx,ib)
       dimension ib(1)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex TRUE")
       write(10,24)
   24  format("coordIndex [")
c
c      do 10 i=1,nx
c      ib(i)=i-1
c  10  continue
c
       write(10,102) (ib(i),i=1,nx)
  102  format(10i5)
c
       return
       end
c
c      zsub12.f  cube subroutine program
c      ---------------------
c      subroutine csub2
c      ---------------------
       subroutine csub2(ic,ib,aa,p,pt,ct)
       parameter(nol=12,nop=3*nol)
c      dimension pp(nop),p(8),nf(nol,6),ii(6),jj(6)
       dimension pp(nop),p(8),nf(nol,6),ii(10),jj(10)
       dimension ic(1),ib(1),aa(1),pt(1),ct(1)
c
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
       u1=(ua-p(1))*(ua-p(2))
       if(u1.lt.0.0) then
       nf(1,1)=1
       nf(1,3)=1
       x=(x1*(p(2)-ua)+x2*(ua-p(1)))/(p(2)-p(1))
       y=y1
       z=z1
       i1=1
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(2))*(ua-p(3))
       if(u1.lt.0.0) then
       nf(2,1)=1
       nf(2,4)=1
       x=x2
       y=(y1*(p(3)-ua)+y2*(ua-p(2)))/(p(3)-p(2))
       z=z1
       i1=1+3*1
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
ccc
c      do 242 j1=1,12
c      i1=1+3*(j1-1)
c      i2=3*j1
c      write(6,22) j1,(pp(ia),ia=i1,i2)
c 242  continue
c
       u1=(ua-p(3))*(ua-p(4))
       if(u1.lt.0.0) then
       nf(3,1)=1
       nf(3,5)=1
       x=(x2*(p(4)-ua)+x1*(ua-p(3)))/(p(4)-p(3))
       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
       u1=(ua-p(4))*(ua-p(1))
       if(u1.lt.0.0) then
       nf(4,1)=1
       nf(4,6)=1
       x=x1
       y=(y2*(p(1)-ua)+y1*(ua-p(4)))/(p(1)-p(4))
       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
       u1=(ua-p(5))*(ua-p(6))
       if(u1.lt.0.0) then
       nf(5,2)=1
       nf(5,3)=1
       x=(x1*(p(6)-ua)+x2*(ua-p(5)))/(p(6)-p(5))
       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
       u1=(ua-p(6))*(ua-p(7))
       if(u1.lt.0.0) then
       nf(6,2)=1
       nf(6,4)=1
       x=x2
       y=(y1*(p(7)-ua)+y2*(ua-p(6)))/(p(7)-p(6))
       z=z2
       i1=1+3*5
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(7))*(ua-p(8))
       if(u1.lt.0.0) then
       nf(7,2)=1
       nf(7,5)=1
       x=(x2*(p(8)-ua)+x1*(ua-p(7)))/(p(8)-p(7))
       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
       u1=(ua-p(8))*(ua-p(5))
       if(u1.lt.0.0) then
       nf(8,2)=1
       nf(8,6)=1
       x=x1
       y=(y2*(p(5)-ua)+y1*(ua-p(8)))/(p(5)-p(8))
       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
       u1=(ua-p(1))*(ua-p(5))
       if(u1.lt.0.0) then
       nf(9,3)=1
       nf(9,6)=1
       x=x1
       y=y1
       z=(z1*(p(5)-ua)+z2*(ua-p(1)))/(p(5)-p(1))
       i1=1+3*8
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(2))*(ua-p(6))
       if(u1.lt.0.0) then
       nf(10,3)=1
       nf(10,4)=1
       x=x2
       y=y1
       z=(z1*(p(6)-ua)+z2*(ua-p(2)))/(p(6)-p(2))
       i1=1+3*9
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(3))*(ua-p(7))
       if(u1.lt.0.0) then
       nf(11,4)=1
       nf(11,5)=1
       x=x2
       y=y2
       z=(z1*(p(7)-ua)+z2*(ua-p(3)))/(p(7)-p(3))
       i1=1+3*10
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(4))*(ua-p(8))
       if(u1.lt.0.0) then
       nf(12,5)=1
       nf(12,6)=1
       x=x1
       y=y2
       z=(z1*(p(8)-ua)+z2*(ua-p(4)))/(p(8)-p(4))
       i1=1+3*11
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c
c      write(6,20) (p(i1),i1=1,8)
   20  format(1h ,4f10.4)
       do 24 j1=1,12
       i1=1+3*(j1-1)
       i2=3*j1
c      write(6,22) j1,(pp(ia),ia=i1,i2)
   22  format(1h ,i5,3f10.4)
   24  continue
       do 28 j1=1,12
c      write(6,26) j1,(nf(j1,ia),ia=1,6)
   26  format(1h ,i5,6i5)
   28  continue
c
c
       call csub1(ic,ib,p,pp,nf,pt,ct)
c
c 100  continue
       return
       end
c
c       
c
c      cube subroutine program
c      ---------------------
c      subroutine csub1
c      ---------------------
       subroutine csub1(ic,ib,p,pp,nf,pt,ct)
       parameter(nol=12,nop=3*nol)
c      dimension p(8),pp(nop),nf(nol,6),ii(6),jj(6),ic(1)
       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,2).eq.1) then
       iin=iin+1
       ii(iin)=i
       jj(iin)=2
       end if
   12  continue
cc
cc
       if(iin.eq.0) go to 100
       if(iin.eq.2.and.ii(2).le.4) go to 200
       if(iin.eq.2.and.ii(2).ge.5) go to 300
       if(iin.eq.4) go to 400
       go to 990
c
  100  continue
       do 110 i=9,12
       j=i-6
       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
       do 210 j=3,6
       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=3,6
       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
       do 230 j=3,6
       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
  300  continue
       i1=ii(1)
       i2=ii(2)
       j1=jj(1)
       j2=jj(2)
c
       do 310 j=3,6
       if(nf(i2,j).eq.1) then
       j3=j
       end if
  310  continue
       do 312 i=9,12
       if(i.eq.i2) go to 312
       if(nf(i,j3).eq.1) then
       iin=iin+1
       i3=i
       ii(iin)=i3
       jj(iin)=j3
       end if
  312  continue
c
       i4=0
       do 320 j=3,6
       if(j.eq.j3) go to 320
       if(nf(i3,j).eq.1) then
       j4=j
       end if
  320  continue
       do 322 i=9,12
       if(i.eq.i2) go to 322
       if(i.eq.i3) go to 322
       if(nf(i,j4).eq.1) then
       iin=iin+1
       i4=i
       ii(iin)=i4
       jj(iin)=j4
       end if
  322  continue
       if(i4.eq.0) go to 1000
c
       do 330 j=3,6
       if(j.eq.j3) go to 330
       if(j.eq.j4) go to 330
       if(nf(i4,j).eq.1) then
       j5=j
       end if
  330  continue
       do 332 i=9,12
       if(i.eq.i2) go to 332
       if(i.eq.i3) go to 332
       if(i.eq.i4) go to 332
       if(nf(i,j5).eq.1) then
       iin=iin+1
       i5=i
       ii(iin)=i5
       jj(iin)=j5
       end if
  332  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=3,6
       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
       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
       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
       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
c      write(6,1020) j,ii1,jj1,x,y,z
 1020  format(1h ,3i4,3f10.4)
 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      cube subroutime program
c      ---------------------
c      subroutine csuba
c      ---------------------
       subroutine csuba(ic,aa,u)
c      isosurface plot
c
c      ic(10):
c      ic(1)=ico,ic(2)=nx,ic(3)=ny,ic(4)=nz
c      aa(20):
c      aa(1)=ua,aa(2)=x1,aa(3)=x2,aa(4)=y1,aa(5)=y2
c      aa(6)=z1,aa(7)=z2,aa(8)=xl,aa(9)=yl,aa(10)=zl
c      aa(11)=x0,aa(12)=y0,aa(13)=z0
c      u(nx*ny*nz): input data
c
       parameter(ipt=9999,ipp=ipt*3)
       parameter(nol=12,nop=3*nol)
c      dimension pp(nop),p(8),nf(nol,6)
       dimension pp(nop),p(8),nf(nol,6)
       dimension ic(1),aa(1),u(1)
       dimension ib(ipt),ct(ipt),pt(ipp)
c
c
c      dimension ic(1),aa(1),u(1),pp(nop),nf(nol,6)
c      dimension p(8)
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
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
       i1=i+n1*(j-1)+n2*(k-1)
       i2=i1+1
       i3=i2+nx
       i4=i1+nx
       i5=i1+n2
       i6=i5+1
       i7=i6+nx
       i8=i5+nx
       p(1)=u(i1)
       p(2)=u(i2)
       p(3)=u(i3)
       p(4)=u(i4)
       p(5)=u(i5)
       p(6)=u(i6)
       p(7)=u(i7)
       p(8)=u(i8)
c
c
       aa(1)=ua
       aa(2)=x1
       aa(3)=x2
       aa(4)=y1
       aa(5)=y2
       aa(6)=z1
       aa(7)=z2
c
c
c      write(6,12) (ic(i1),i1=1,10)
c      write(6,14) (aa(i1),i1=1,20)
c      write(6,16) (p(i1),i1=1,8)
   12  format(1h ,10i5)
   14  format(1h ,5f10.5)
   16  format(1h ,4f10.5)
c
       call csub2(ic,ib,aa,p,pt,ct)
c
  100  continue
c
       iin1=ic(5)
       iin2=ic(6)
c      write(6,12) k,iin1,iin2
c
       if(iin1.eq.0) go to 200
c
       call facebe
c
       call colorp(iin1,ct,ico)
c
       call coordbe
c
       call zcoordp(iin1,pt)
c
       iin21=iin2-1
       call zcoorden(iin21,ib)
       call faceen

c
  200  continue
       return
       end
c
c
c       
c************************************************************
c           subroutine file     << msubrtn2.f >>            *
c                                                           * 
c      update ; 01/03/24     /home/g3/inayoshi/VRML/subrtn4 *
c************************************************************
c-------------------------------------
       subroutine triangm1(nx,ny,u,p)
c-------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       call facebe
c
       do 10 j=1,ny1
       call colorp2(nx,j,u)
   10  continue
c
       call coordbe
c
       do 20 j=1,ny1
       call coordp2(nx,j,p)
   20  continue
c
       call coorden5(nx,ny1)
c
       call ifaceen
c
       return
       end
c-------------------------------------
       subroutine triangm(nx,ny,u,p)
c-------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call triangb(nx,j,u,p)
   10  continue  
c
       return
       end
c-------------------------------------
       subroutine triangb(nx,ny,u,p)
c-------------------------------------
       dimension p(1),u(1)
c
       call facebe
       call colorp2(nx,ny,u)
       call coordbe
       call coordp2(nx,ny,p)
       call coorden1(nx)
       call ifaceen
c
       return
       end
c-------------------------------------------
       subroutine linemesh(p,n1,rd,r,g,b)
c-------------------------------------------
       dimension p(1)
       dimension p1(3*2)
       dimension q(3*12*2)
c
       n11=n1-2
       do 20 i=1,n11
       do 10 j=1,2
c
       i1=3*(i-1)+1+3*(j-1)
       i2=i1+1
       i3=i1+2
c---<< chusin ten set >>
       p1(1+(j-1)*3)=p(i1)
       p1(2+(j-1)*3)=p(i2)
       p1(3+(j-1)*3)=p(i3)
c
       xa=abs(p(i1)-p(i1+3))
       ya=abs(p(i2)-p(i2+3))
       za=abs(p(i3)-p(i3+3))
       amx=amax1(xa,ya,za)  
c
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       if(amx.eq.xa)  then
               k=1
        else
          if(amx.eq.ya)  then
                 k=2
           else
             if(amx.eq.za)  then
                   k=3
             end if
          end if
       end if
c
c---<< ensyu no point wo motomeru >>
c
       call ensyu(x,y,z,rd,q,j,k)
   10  continue
c
       call linemeshb(12,1,q,rd,r,g,b,k,p1)
   20  continue
c---<< cone >>
       rr=1.57
       if(k.eq.1)   then
            x1=0.0
            y1=0.0
            z1=-1.0
       end if
       if(k.eq.3)   then
            x1=-1.0
            y1=0.0
            z1=0.0
       end if
c
       rd2=rd*2
       h=rd*3
       x=p1(4)
       y=p1(5)
       z=p1(6)
       call cone1(rd2,h,r,g,b,x,y,z,x1,y1,z1,rr)
c
       return
       end
c--------------------------------------------------
       subroutine linemeshb(nx,ny,q,rd,r,g,b,k,p1)
c                 * 1 color *
c--------------------------------------------------
       dimension q(1)
c
       call facebe
c
c----( color )
       write(10,10) r,g,b
   10  format(f9.2,1x,f9.2,1x,f9.2)
c
       call coordbe1(nx)
       call coordp4(nx,ny,q,rd,k,p1)
       call coorden3(nx)
       call ifaceen
c
       return
       end
c---------------------------------------
       subroutine ensyu(x,y,z,rd,q,j,k)
c---------------------------------------
       dimension q(1)
c
       x1=rd/2*1.73205
       y1=rd/2
       j2=36*(j-1)
c--[ max=xa ]
       if(k.eq.1)    then
c (x)
       do 10 i=1,12
       q(1+(i-1)*3+j2)=x
   10  continue
c (y)
       q(2+j2)=y+rd
       q(5+j2)=y+x1
       q(8+j2)=y+y1
       q(11+j2)=y
       q(14+j2)=y-y1
       q(17+j2)=y-x1
       q(20+j2)=y-rd
       q(23+j2)=y-x1
       q(26+j2)=y-y1
       q(29+j2)=y
       q(32+j2)=y+y1
       q(35+j2)=y+x1
c (z)
       q(3+j2)=z
       q(6+j2)=z+y1
       q(9+j2)=z+x1
       q(12+j2)=z+rd
       q(15+j2)=z+x1
       q(18+j2)=z+y1
       q(21+j2)=z
       q(24+j2)=z-y1
       q(27+j2)=z-x1
       q(30+j2)=z-rd
       q(33+j2)=z-x1
       q(36+j2)=z-y1
       go to 99
       end if
c--[ max=ya ]
       if(k.eq.2)    then
c (x)
       q(1+j2)=x
       q(4+j2)=x+y1
       q(7+j2)=x+x1
       q(10+j2)=x+rd
       q(13+j2)=x+x1
       q(16+j2)=x+y1
       q(19+j2)=x
       q(22+j2)=x-y1
       q(25+j2)=x-x1
       q(28+j2)=x-rd
       q(31+j2)=x-x1
       q(34+j2)=x-y1
c (y)
       do 20 i=1,12
       q(2+(i-1)*3+j2)=y
   20  continue
c (z)
       q(2+j2)=z+rd
       q(5+j2)=z+x1
       q(8+j2)=z+y1
       q(11+j2)=z
       q(14+j2)=z-y1
       q(17+j2)=z-x1
       q(20+j2)=z-rd
       q(23+j2)=z-x1
       q(26+j2)=z-y1
       q(29+j2)=z
       q(32+j2)=z+y1
       q(35+j2)=z+x1
       go to 99
       end if
c--[ max=za ]
       if(k.eq.3)    then
c (x)
       q(1+j2)=x
       q(4+j2)=x+y1
       q(7+j2)=x+x1
       q(10+j2)=x+rd
       q(13+j2)=x+x1
       q(16+j2)=x+y1
       q(19+j2)=x
       q(22+j2)=x-y1
       q(25+j2)=x-x1
       q(28+j2)=x-rd
       q(31+j2)=x-x1
       q(34+j2)=x-y1
c (y)
       q(2+j2)=y+rd
       q(5+j2)=y+x1
       q(8+j2)=y+y1
       q(11+j2)=y
       q(14+j2)=y-y1
       q(17+j2)=y-x1
       q(20+j2)=y-rd
       q(23+j2)=y-x1
       q(26+j2)=y-y1
       q(29+j2)=y
       q(32+j2)=y+y1
       q(35+j2)=y+x1
c (z)
       do 30 i=1,12
       q(3*i+j2)=z
   30  continue
       end if
c
   99  continue
       return
       end
c--------------------------------------------------
       subroutine cylinderm(nx,ny,u,p,rd,br,bg,bb)
c                 * gradation color *
c--------------------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call cylinderb(nx,j,u,p,rd,br,bg,bb)
   10  continue
c
       return
       end
c--------------------------------------------------
       subroutine cylinderb(nx,ny,u,p,rd,br,bg,bb)
c                 * gradation color *
c--------------------------------------------------
       dimension p(1),u(1)
c
       call facebe
       call colorp3(nx,ny,u,br,bg,bb)
       call coordbe
       call coordp3(nx,ny,p,rd)
       call coorden2(nx)
       call ifaceen
c
       return
       end
c--------------------------------------------------
       subroutine cylinderm1(nx,ny,p,rd,r,g,b)
c                 * 1 color *
c--------------------------------------------------
       dimension p(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call cylinderb1(nx,j,p,rd,r,g,b)
   10  continue
c
       return
       end
c--------------------------------------------------
       subroutine cylinderb1(nx,ny,p,rd,r,g,b)
c                 * 1 color *
c--------------------------------------------------
       dimension p(1)
c
       call facebe
c----( color )
       write(10,10) r,g,b
   10  format(f8.2,1x,f8.2,1x,f8.2)
c
       call coordbe1(nx)
       call coordp3(nx,ny,p,rd)
       call coorden3(nx)
       call ifaceen
c
       return
       end
c----------------------------------------------------------------------
       subroutine arrow(nx,ny,u,p,rd1,rd2,h,cx,cy,cz,br,bg,bb,cr,cg,cb)
c----------------------------------------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call cylinderb(nx,j,u,p,rd1,br,bg,bb)
   10  continue
c
       x1=0.0
       y1=0.0
       z1=1.0
       rr=1.57
c
       call cone1(rd2,h,cr,cg,cb,cx,cy,cz,x1,y1,z1,rr)
       return
       end
c----------------------------------------------------------------------
       subroutine arrow2(nx,ny,u,p,rd1,rd2,h,cx,cy,cz,br,bg,bb,cr,cg,cb)
c----------------------------------------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call cylinderb(nx,j,u,p,rd1,br,bg,bb) 
   10  continue
c
       jj=ny1
       i=nx
       j=2
       j1=j+jj-1
       i1=3*(i-1)+1+3*nx*(j1-1)
       i2=i1+1
       i3=i1+2
       cx=p(i1)+0.8*rd2
       cy=p(i2)
       cz=p(i3)+rd1
c
       x1=0.0
       y1=0.0
       z1=-1.0
       rr=1.57
c
       call cone1(rd2,h,cr,cg,cb,cx,cy,cz,x1,y1,z1,rr)
       return
       end
c-----------------------------------------------------
       subroutine cone1(rd,h,r,g,b,x,y,z,x1,y1,z1,rr)
c-----------------------------------------------------
       write(10,20) x,y,z
   20  format("Transform {",/"translation ",f7.2,f7.2,f7.2)
       write(10,21) x1,y1,z1,rr
   21  format("rotation ",f7.2,f7.2,f7.2,f7.2)
       write(10,22)
   22  format("children [")
       write(10,23)
   23  format("Shape {")
       write(10,24)
   24  format("appearance Appearance {")
       write(10,25)
   25  format("material Material {")
       write(10,26) r,g,b
   26  format("diffuseColor",f5.2,f5.2,f5.2)
       write(10,27)
   27  format("}",/"}")
       write(10,28)
   28  format("geometry Cone {")
       write(10,29) rd
   29  format("bottomRadius  ",f5.2)
       write(10,30) h
   30  format("height  ",f5.2)
       write(10,31)
   31  format("}",/"}",/"]",/"}")
c
       return
       end
c-------------------------------
       subroutine coordbe1(nx)
c-------------------------------
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
c----( colorindexed )
       nx4=nx*4
       write(10,20)
   20  format("colorIndex [ ")
       write(10,30) (0,i=1,nx4)
   30  format(5000i3)
       write(10,40)
   40  format(" ]")
c
       write(10,23)
   23  format("coord Coordinate {")
       write(10,24)
   24  format("point [")
       return
       end
c--------------------------------------
       subroutine colorp1(nx,ny,u)
c--------------------------------------
       dimension u(1)
c
       do 100 i=1,nx
c
       i1=3*(i-1)+1+3*nx*(ny-1)
       i2=i1+1
       i3=i1+2
       r=u(i1)
       g=u(i2)
       b=u(i3)
c
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       if(i.eq.nx) then
       write(10,10) r,g,b
   10  format(f6.2,1x,f6.2,1x,f6.2)
       else
       write(10,20) r,g,b
   20  format(f6.2,1x,f6.2,1x,f6.2,",")
       end if
c
  100  continue
       return
       end
c--------------------------------------
       subroutine colorp2(nx,ny,u)
c--------------------------------------
       dimension u(1)
c
       do 100 i=1,nx
       do 100 j=1,2
       j1=j+ny-1
c
       i1=3*(i-1)+1+3*nx*(j1-1)
       i2=i1+1
       i3=i1+2
       r=u(i1)
       g=u(i2)
       b=u(i3)
c
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       if(i.eq.nx.and.j.eq.2) then
       write(10,10) r,g,b
   10  format(f6.2,1x,f6.2,1x,f6.2)
       else
       write(10,20) r,g,b
   20  format(f6.2,1x,f6.2,1x,f6.2,",")
       end if
c
  100  continue
       return
       end
c--------------------------------------------
       subroutine colorp3(nx,ny,u,br,bg,bb)
c                        * bottom tuki *
c--------------------------------------------
       dimension u(1)
c
       do 100 i=1,nx
       do 100 j=1,2
       j1=j+ny-1
c
       i1=3*(i-1)+1+3*nx*(j1-1)
       i2=i1+1
       i3=i1+2
       r=u(i1)
       g=u(i2)
       b=u(i3)
c
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       write(10,10) r,g,b
   10  format(f8.2,1x,f8.2,1x,f8.2,",")
c
  100  continue
c----( bottom chusin color )
       write(10,120) br,bg,bb
  120  format(f8.2,1x,f8.2,1x,f8.2)
c
       return
       end
c-------------------------------------
       subroutine coordp1(nx,ny,p)
c-------------------------------------
       dimension p(1)
c
       do 100 i=1,nx
c
       i1=3*(i-1)+1+3*nx*(ny-1)
       i2=i1+1
       i3=i1+2
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       if(i.eq.nx) then
       write(10,10) x,y,z
   10  format(f9.2,1x,f9.2,1x,f9.2)
       else
       write(10,20) x,y,z
   20  format(f9.2,1x,f9.2,1x,f9.2,",")
       end if
c
  100  continue
       return
       end
c-------------------------------------
       subroutine coordp2(nx,ny,p)
c-------------------------------------
       dimension p(1)
c
       do 100 i=1,nx
       do 100 j=1,2
       j1=j+ny-1
c
       i1=3*(i-1)+1+3*nx*(j1-1)
       i2=i1+1
       i3=i1+2
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       if(i.eq.nx.and.j.eq.2) then
       write(10,10) x,y,z
   10  format(f9.2,1x,f9.2,1x,f9.2)
       else
       write(10,20) x,y,z
   20  format(f9.2,1x,f9.2,1x,f9.2,",")
       end if
c
  100  continue
       return
       end
c-------------------------------------------
       subroutine coordp3(nx,ny,p,rd)
c                 * bottom tuki *
c-------------------------------------------
       dimension p(1)
c
       do 100 i=1,nx
       do 100 j=1,2
       j1=j+ny-1
c
       i1=3*(i-1)+1+3*nx*(j1-1)
       i2=i1+1
       i3=i1+2
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       write(10,20) x,y,z
   20  format(f9.2,1x,f9.2,1x,f9.2,",")
c
  100  continue
c
c----( bottom no chusin )
c
       i1=3*nx*2+1
       i2=i1+1
       i3=i1+2
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       write(10,110) x,y,z
  110  format(f9.2,1x,f9.2,1x,f9.2)
c
       return
       end
c-------------------------------------------
       subroutine coordp4(nx,ny,q,rd,k,p1)
c                 * bottom tuki *
c-------------------------------------------
       dimension q(1)
       dimension p1(1)
c
       do 100 i=1,nx
       do 100 j=1,2
c
       i1=3*(i-1)+1+3*nx*(j-1)
       i2=i1+1
       i3=i1+2
       x=q(i1)
       y=q(i2)
       z=q(i3)
c
       write(10,20) x,y,z
   20  format(f9.2,1x,f9.2,1x,f9.2,",")
c
  100  continue
c
c----( bottom no chusin )
       x=p1(1)
       y=p1(2)
       z=p1(3)
c
       write(10,30) x,y,z
   30  format(f9.2,1x,f9.2,1x,f9.2,",")
c
c----( cap no chusin )
       x=p1(4)
       y=p1(5)
       z=p1(6)
c
       write(10,40) x,y,z
   40  format(f9.2,1x,f9.2,1x,f9.2)
c
       return
       end
c-------------------------------
       subroutine coorden1(nx)
c-------------------------------
       dimension ia(9000)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex TRUE")
       write(10,24)
   24  format("coordIndex [")
c
       nx2=nx*2
       do 100 i=1,nx2
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=0+i-1
       ia(i2)=1+i-1
       ia(i3)=2+i-1
       ia(i4)=-1
  100  continue
c
       jj=4*(nx-1)*2-1
       write(10,102) (ia(i),i=1,jj)
  102  format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3,
     1        398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3,
     2        /500(3i5,i3),/500(3i5,i3),/500(3i5,i3),/500(3i5,i3))
       return
       end
c-----------------------------------
       subroutine coorden2(nx)
c                 * bottom tuki *
c-----------------------------------
       dimension ia(9000)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex TRUE")
       write(10,24)
   24  format("coordIndex [")
c
       nx2=nx*2
       do 100 i=1,nx2
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=0+i-1
       ia(i2)=1+i-1
       ia(i3)=2+i-1
       ia(i4)=-1
  100  continue
c
       jj=4*(nx-1)*2
       write(10,102) (ia(i),i=1,jj)
  102  format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3,
     1        398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3,
     2        /500(3i5,i3),/500(3i5,i3),/500(3i5,i3),/500(3i5,i3))
c
       write(10,103) ia(jj-2),ia(jj-1),ia(1)
  103  format(3i5," -1")
       write(10,104) ia(jj-1),ia(1),ia(2)
  104  format(3i5," -1")
c
c---- bottom
       nx2=nx*2
       do 200 i=1,nx
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=0+2*(i-1)
       ia(i2)=ia(i1)+2
       ia(i3)=nx2
       ia(i4)=-1
  200  continue
c
       jj=4*(nx-1)
       write(10,202) (ia(i),i=1,jj)
  202  format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3))
c
       write(10,203) ia(jj-2),ia(1),nx2
  203  format(3i4," -1")
c
c---- cap
       nx2=nx*2
       do 300 i=1,nx
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=1+2*(i-1)
       ia(i2)=ia(i1)+2
       ia(i3)=nx2
       ia(i4)=-1
  300  continue
c
       jj=4*(nx-1)
       write(10,302) (ia(i),i=1,jj)
  302  format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3))
c
       write(10,303) ia(jj-2),ia(1),nx2
  303  format(3i4," -1")
c
       return
       end
c-----------------------------------
       subroutine coorden3(nx)
c                 * bottom tuki *
c                 * 1 color     *
c-----------------------------------
       dimension ia(9000)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex FALSE")
       write(10,24)
   24  format("coordIndex [")
c
       nx2=nx*2
       do 100 i=1,nx2
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=0+i-1
       ia(i2)=1+i-1
       ia(i3)=2+i-1
       ia(i4)=-1
  100  continue
c
       jj=4*(nx-1)*2
       write(10,102) (ia(i),i=1,jj)
  102  format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3,
     1        398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3,
     2        /500(3i5,i3),/500(3i5,i3),/500(3i5,i3),/500(3i5,i3))
c
       write(10,103) ia(jj-2),ia(jj-1),ia(1)
  103  format(3i5," -1")
       write(10,104) ia(jj-1),ia(1),ia(2)
  104  format(3i5," -1")
c
c---- bottom
       nx2=nx*2
       do 200 i=1,nx
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=0+2*(i-1)
       ia(i2)=ia(i1)+2
       ia(i3)=nx2
       ia(i4)=-1
  200  continue
c
       jj=4*(nx-1)
       write(10,202) (ia(i),i=1,jj)
  202  format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3))
c
       write(10,203) ia(jj-2),ia(1),nx2
  203  format(3i4," -1")
c
c---- cap
       nx21=nx*2+1
       do 300 i=1,nx
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=1+2*(i-1)
       ia(i2)=ia(i1)+2
       ia(i3)=nx21
       ia(i4)=-1
  300  continue
c
       jj=4*(nx-1)
       write(10,302) (ia(i),i=1,jj)
  302  format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3))
c
       write(10,303) ia(jj-2),ia(1),nx21
  303  format(3i4," -1")
c
       return
       end
c----------------------------------
       subroutine coorden5(nx,ny1)
c----------------------------------
       dimension ia(9000)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex TRUE")
       write(10,24)
   24  format("coordIndex [")
c
cc      nx2=(nx-1)*2
       nx2=nx*2
       do 200 j=1,ny1
c
       do 100 i=1,nx2
       i1=1+4*(i-1)
       i2=2+4*(i-1)
       i3=3+4*(i-1)
       i4=4+4*(i-1)
       ia(i1)=(j-1)*2*nx+0+i-1
       ia(i2)=(j-1)*2*nx+1+i-1
       ia(i3)=(j-1)*2*nx+2+i-1
       ia(i4)=-1
  100  continue
c
       jj=4*(nx-1)*2-1
       write(10,102) (ia(i),i=1,jj)
c 102  format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3,
c    1        398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3,
  102  format(250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     1        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     2        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     3        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     4        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     5        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     6        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     7        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3),
     8        /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3))
c
  200  continue
       return
       end
c------------------------
       subroutine ifaceen
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("transparency 0.000000",/"}",/"}")
       write(10,25)
   25  format("}")
       return
       end
c--------------------------
       subroutine faceen2
c--------------------------
       write(10,20)
   20  format("]")
       write(10,21)
   21  format("solid FALSE",/"}")
       write(10,25)
   25  format("}")
       return
       end
c--------------------------------------
       subroutine pointlight(px,py,pz)
c--------------------------------------
       write(10,20)
   20  format("PointLight {")
       write(10,21) px,py,pz
   21  format("location ",f8.2,f8.2,f8.2)
       write(10,22)
   22  format("}")
       return
       end
c-----------------------------------------
       subroutine directlight(px,py,pz)
c-----------------------------------------
       write(10,20)
   20  format("DirectionalLight {")
       write(10,21) px,py,pz
   21  format("direction ",f8.2,f8.2,f8.2)
       write(10,22)
   22  format("}")
       return
       end
c-------------------------------------
       subroutine lineset1(nx,ny,u,p)
c-------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny
       call line1(nx,j,u,p)
   10  continue
c
       return
       end
c-------------------------------------
       subroutine line1(nx,ny,u,p)
c-------------------------------------
       dimension p(1),u(1)
c
       call ilinebe
       call colorp1(nx,ny,u)
       call coordbe
       call coordp1(nx,ny,p)
       call ilineen(nx)
c
       return
       end
c------------------------
       subroutine ilinebe
c------------------------
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("geometry IndexedLineSet {")
       write(10,23)
   23  format("color Color {")
       write(10,24)
   24  format("color [")
       return
       end
c-------------------------------
       subroutine ilineen(nx)
c-------------------------------
       dimension ia(9000)
c
       do 10 i=1,9000
       ia(i)=i-1
   10  continue
c
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,24)
   24  format("coordIndex [")
c
       write(10,25) (ia(i),i=1,nx)
   25  format(10i2,90i3,900i4/(2000i5))
c
       write(10,27)
   27  format("]")
       write(10,29)
       write(10,29)
   29  format("}")
       return
       end
c----------------------------------------
       subroutine imagetx(nx,ny,p,ch1,n)
c----------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call imagea(nx,j,p,ch1,n)
   10  continue
c
       return
       end
c---------------------------------------
       subroutine imagea(nx,ny,p,ch1,n)
c---------------------------------------
       dimension p(1),u(1)
c
       call imagebe(ch1,n)
       call coordp2(nx,ny,p)
       call coorden1(nx)
       call ifaceen
c
       return
       end
c-------------------------------
      subroutine imagebe(chr,n)
c-------------------------------
       character chr*15,ica*15,ich(15)*1
       equivalence (ica,ich)
       ica=chr
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("appearance Appearance {")
       write(10,23)
   23  format("texture ImageTexture {")
   24  write(10,*) "url",'"',(ich(i),i=1,n),'"'
       write(10,25)
   25  format("}",/"}")
       write(10,26)
   26  format("geometry IndexedFaceSet {")
       write(10,27)
   27  format("coord Coordinate {")
       write(10,28)
   28  format("point [")
       return
       end
c----------------------------------------
       subroutine movietx(nx,ny,p,ch1,n)
c----------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call movie1(nx,j,p,ch1,n)
   10  continue
c
       return
       end
c---------------------------------------
       subroutine movie1(nx,ny,p,ch1,n)
c---------------------------------------
       dimension p(1),u(1)
c
       call moviebe(ch1,n)
       call coordp2(nx,ny,p)
       call coorden1(nx)
       call ifaceen
c
       return
       end
c-------------------------------
      subroutine moviebe(chr,n)
c-------------------------------
       character chr*15,ica*15,ich(15)*1
       equivalence (ica,ich)
       ica=chr
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("appearance Appearance {")
       write(10,23)
   23  format("texture MovieTexture {")
   24  write(10,*) "url",'"',(ich(i),i=1,n),'"'
       write(10,30)
cc 30  format("loop TRUE",/"startTime 0")
   30  format("loop FALSE",/"startTime 0",/"stopTime 0")
       write(10,25)
   25  format("}",/"}")
       write(10,26)
   26  format("geometry IndexedFaceSet {")
       write(10,27)
   27  format("coord Coordinate {")
       write(10,28)
   28  format("point [")
       return
       end
c-----------------------------------------------------------
       subroutine defusem(nx,ny,u,p,x1,y1,z1,x2,y2,z2,r,n)
c----------------------------------------------------------
       dimension p(1),u(1)
       ny1=ny-1
c
       do 10 j=1,ny1
       call defusem1(nx,j,u,p,x1,y1,z1,x2,y2,z2,r,n)
   10  continue
c
       return
       end
c-----------------------------------------------------------
       subroutine defusem1(nx,ny,u,p,x1,y1,z1,x2,y2,z2,r,n)
c-----------------------------------------------------------
       dimension p(1),u(1)
c
       call defface
       call colorp2(nx,ny,u)
       call coordbe
       call coordp2(nx,ny,p)
       call coorden1(nx)
       call ifaceen
c
       do 10 i=1,n
       call useface(x1,y1,z1,x2,y2,z2,r,i)
   10  continue
c
       return
       end
c--------------------------
       subroutine defface
c--------------------------
       write(10,21)
   21  format("DEF Face 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 useface(x1,y1,z1,x2,y2,z2,r,i)
c------------------------------------------------
       write(10,21)
   21  format("Transform {")
       write(10,22) x1,y1,z1,r
   22  format("rotation ",f6.2,f6.2,f6.2,f6.2)
       write(10,23) x2,y2*i,z2
   23  format("translation ",f9.2,f9.2,f9.2)
       write(10,24)
   24  format("children USE Face",/"}")
       return
       end
c-----------------------------------------
       subroutine sphere2(p,n1,rd,r,g,b)
c-----------------------------------------
       dimension p(1)
c
       do 20 i=1,n1
c
       i1=3*(i-1)+1
       i2=i1+1
       i3=i1+2
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       call sphere(x,y,z,r,g,b,rd)
c
   20  continue
c
       return
       end
c----------------------------------------
       subroutine sphere(x,y,z,r,g,b,rd)
c----------------------------------------
       write(10,20)
   20  format("Transform {")
       write(10,21) x,y,z
   21  format("translation ",f9.2,f9.2,f9.2)
       write(10,22)
   22  format("children [")
       write(10,23)
   23  format("Shape {")
       write(10,24)
   24  format("appearance Appearance {")
       write(10,25)
   25  format("material Material {")
       write(10,26)
   26  format("transparency 0.000000")
       write(10,27) r,g,b
   27  format("diffuseColor ",f5.2,f5.2,f5.2)
       write(10,28)
   28  format("}",/"}")
       write(10,29)
   29  format("geometry Sphere {")
       write(10,30) rd
   30  format("radius ",f6.2)
       write(10,31)
   31  format("}",/"}",/"]",/"}")
c
       return
       end
c-------------------------------------------------------
       subroutine cylinder1(p,n1,rx,ry,rz,dd,rd,h,r,g,b)
c-------------------------------------------------------
       dimension p(1)
c
       do 20 i=1,n1
c
       i1=3*(i-1)+1
       i2=i1+1
       i3=i1+2
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       call cylinder(x,y,z,rx,ry,rz,dd,rd,h,r,g,b)
c
   20  continue
c
       return
       end
c-----------------------------------------------
       subroutine cylinderm2(p,n1,rd,r,g,b,h)
c-----------------------------------------------
       dimension p(1)
       rx=0.0
       ry=0.0
       rz=0.0
       dd=0.78539
c
       do 20 i=1,n1
c
       i1=3*(i-1)+1
       i2=i1+1
       i3=i1+2
       x=p(i1)
       y=p(i2)
       z=p(i3)
c
       x1=p(i1+3)
       y1=p(i2+3)
       z1=p(i3+3)
c
       x2=abs(x1-x)
       y2=abs(y1-y)
       z2=abs(z1-z)
c
c---<< x houkou >>
       if(x2.gt.y2.and.x2.gt.z2) then
          rz=-1
          if(z2.gt.y2) then
              rx=1
          end if
       end if
c---<< y houkou >>
c      if(y2.gt.z2.and.y2.gt.x2) then
c                rx=-1
c         else
c                rx=1
c         end if
c      end if
c---<< z houkou >>
       if(z2.gt.x2.and.z2.gt.y2) then
          if(x2.lt.y2) then
              rx=-1
          else
              rx=1
          end if  
       end if
c
       call cylinder(x,y,z,rx,ry,rz,dd,rd,h,r,g,b)
c      call def1(x,y,z,rx,ry,rz,dd,rd,h,r,g,b)
c
   20  continue
c
   99  continue
c
       return
       end
c-----------------------------------------------------------
       subroutine cylinder(cx,cy,cz,rx,ry,rz,dd,rd,t,r,g,b)
c-----------------------------------------------------------
       write(10,20)
   20  format("Transform {")
       write(10,21) rx,ry,rz,dd
   21  format("rotation ",f7.2,f7.2,f7.2,f7.4)
       write(10,22) cx,cy,cz
   22  format("translation ",f7.2,f7.2,f7.2)
       write(10,23)
   23  format("children Shape {")
       write(10,24)
   24  format("appearance Appearance {")
       write(10,25)
   25  format("material Material {")
       write(10,26)
   26  format("transparency 0.000000")
       write(10,27) r,g,b
   27  format("diffuseColor",f5.2,f5.2,f5.2)
       write(10,28)
   28  format("}",/"}")
       write(10,29)
   29  format("geometry Cylinder {")
       write(10,30) rd
   30  format("radius",f7.2)
       write(10,31) t
   31  format("height",f7.2)
       write(10,32)
   32  format("}",/"}",/"}")
c
       return
       end
c-------------------------------------------------------
       subroutine defl(cx,cy,cz,rx,ry,rz,dd,rd,t,r,g,b)
c-------------------------------------------------------
       write(10,20)
   20  format("Transform {")
       write(10,21) rx,ry,rz,dd
   21  format("rotation ",f9.2,f9.2,f9.2,f9.2)
       write(10,22) cx,cy,cz
   22  format("translation ",f9.2,f9.2,f9.2)
       write(10,23)
   23  format("children DEF def1 Shape {")
       write(10,24)
   24  format("appearance Appearance {")
       write(10,25)
   25  format("material Material {")
       write(10,26)
   26  format("transparency 0.000000")
       write(10,27) r,g,b
   27  format("diffuseColor",f5.2,f5.2,f5.2)
       write(10,28)
   28  format("}",/"}")
       write(10,29)
   29  format("geometry Cylinder {")
       write(10,30) rd
   30  format("radius",f8.2)
       write(10,31) t
   31  format("height",f8.2)
       write(10,32)
   32  format("}",/"}",/"}")
c
       return
       end
c------------------------------------------------
       subroutine use1(cx,cy,cz,rx,ry,rz,dd)
c------------------------------------------------
       write(10,20)
   20  format("Transform {")
       write(10,21) rx,ry,rz,dd
   21  format("rotation ",f9.2,f9.2,f9.2,f9.2)
       write(10,22) cx,cy,cz
   22  format("translation ",f9.2,f9.2,f9.2)
       write(10,23)
   23  format("children USE def1")
       write(10,32)
   32  format("}")
c
       return
       end
c      ---------------------
c      subroutine pixel2
c      ---------------------
       subroutine pixel2(nx,ny,xb,yb,xl,yl,ipx0,ico,icc,zcc,u)
c
c      nx: number in x-direction
c      nx: number in x-direction
c      x0: x-value at the left and bottom point
c      y0: y-value at the left and bottom point
c      xl: length in x-direction
c      yl: length in y-direction
c      ipx0: 3 for color, 4 for transparency
c      ico: color code
c      icc: selection of plane, i
c         : 1 for xy plane, s2 for xz plane, 3 for yz plane
c      zcc: position of selected plane
c      u(nx*ny): input data
c
       dimension u(1)
       nxy=nx*ny
c
       call pixebe(nx,ny,ipx0)
       if(ipx0.eq.3) call pixeco(nx,ny,ico,u)
       if(ipx0.eq.4) call pixect(nx,ny,ico,u)
       call pixeen(xb,yb,xl,yl,icc,zcc)
c
       return
       end
c      ---------------------
c      subroutine pixebe
c      ---------------------
       subroutine pixebe(nx,ny,ipx0)
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("appearance Appearance {")
       write(10,23)
   23  format("texture PixelTexture {")
       write(10,24) nx,ny,ipx0
   24  format("image ",2i5,i2)
       return
       end
c      ---------------------
c      subroutine pixeco
c      ---------------------
       subroutine pixeco(nx,ny,ico,u)
       dimension u(1),ipx(6)
       r=1.0
       g=1.0
       b=1.0
       do 100 j=1,ny
       do 100 i=1,nx
       i1=i+nx*(j-1)
       x1=u(i1)
c
       a0=1.20
       a0=1.00
       if(ico.eq.2) go to 32
       if(ico.eq.3) go to 33
       if(ico.eq.4) go to 34
       if(ico.eq.5) go to 35
       if(ico.eq.6) go to 36
       if(ico.eq.7) go to 37
       r=a0-(a0-0.5)*x1
       g=a0-(a0-0.5)*x1
       b=a0-(a0-0.5)*x1
       go to 40
   32  continue
       r=a0-a0*x1
       g=1.0
       b=a0-a0*x1
       go to 40
   33  continue
       r=1.0
       g=a0-a0*x1
       b=a0-a0*x1
       go to 40
   34  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       g=a0+3.0-y1
       b=a0+1.0-y1
       if(y1.ge.a0+2.0) b=y1-(a0+3.0)
       go to 40
   35  continue
       a0=2.0
       y1=x1*(a0+4.0)
       b=a0-y1
       if(y1.ge.a0) b=y1-(a0+1.0)
       r=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   36  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       b=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   37  continue
       x1=10.0*x1
       b=x1
       if(x1.ge.3.0) b=4.0-x1
       if(x1.ge.6.0) b=0.25*(x1-6.0)
       g=x1-2.0
       if(x1.ge.5.0) g=6.0-x1
       if(x1.ge.6.0) g=0.25*(x1-6.0)
       r=x1
       if(x1.ge.1.0) r=2.0-x1
       if(x1.ge.4.0) r=x1-4.0
       go to 40
c
   40  continue
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       ir=r*256.0
       ig=g*256.0
       ib=b*256.0
       ir=max(0,ir)
       ir=min(255,ir)
       ig=max(0,ig)
       ig=min(255,ig)
       ib=max(0,ib)
       ib=min(255,ib)
       ir1=ir/16
       ir2=ir-16*ir1
       ig1=ig/16
       ig2=ig-16*ig1
       ib1=ib/16
       ib2=ib-16*ib1
       ipx(1)=ir1
       ipx(2)=ir2
       ipx(3)=ig1
       ipx(4)=ig2
       ipx(5)=ib1
       ipx(6)=ib2
c
       write(10,20) ipx
c  20  format("0x",6z1,1x)
   20  format(" 0x",6z1)
c
  100  continue
       return
       end
c      ---------------------
c      subroutine pixect
c      ---------------------
       subroutine pixect(nx,ny,ico,u)
       dimension u(1),ipx(8)
       r=1.0
       g=1.0
       b=1.0
       do 100 j=1,ny
       do 100 i=1,nx
       i1=i+nx*(j-1)
       x1=u(i1)
c
       a0=1.20
       a0=1.00
       if(ico.eq.2) go to 32
       if(ico.eq.3) go to 33
       if(ico.eq.4) go to 34
       if(ico.eq.5) go to 35
       if(ico.eq.6) go to 36
       if(ico.eq.7) go to 37
       r=a0-(a0-0.5)*x1
       g=a0-(a0-0.5)*x1
       b=a0-(a0-0.5)*x1
       go to 40
   32  continue
       r=a0-a0*x1
       g=1.0
       b=a0-a0*x1
       go to 40
   33  continue
       r=1.0
       g=a0-a0*x1
       b=a0-a0*x1
       go to 40
   34  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       g=a0+3.0-y1
       b=a0+1.0-y1
       if(y1.ge.a0+2.0) b=y1-(a0+3.0)
       go to 40
   35  continue
       a0=2.0
       y1=x1*(a0+4.0)
       b=a0-y1
       if(y1.ge.a0) b=y1-(a0+1.0)
       r=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   36  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       b=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   37  continue
       x1=10.0*x1
       b=x1
       if(x1.ge.3.0) b=4.0-x1
       if(x1.ge.6.0) b=0.25*(x1-6.0)
       g=x1-2.0
       if(x1.ge.5.0) g=6.0-x1
       if(x1.ge.6.0) g=0.25*(x1-6.0)
       r=x1
       if(x1.ge.1.0) r=2.0-x1
       if(x1.ge.4.0) r=x1-4.0
       go to 40
c
   40  continue
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       ir=r*256.0
       ig=g*256.0
       ib=b*256.0
       ir=max(0,ir)
       ir=min(255,ir)
       ig=max(0,ig)
       ig=min(255,ig)
       ib=max(0,ib)
       ib=min(255,ib)
       ir1=ir/16
       ir2=ir-16*ir1
       ig1=ig/16
       ig2=ig-16*ig1
       ib1=ib/16
       ib2=ib-16*ib1
       ipx(1)=ir1
       ipx(2)=ir2
       ipx(3)=ig1
       ipx(4)=ig2
       ipx(5)=ib1
       ipx(6)=ib2
c
       it=u(i1)*256.0
       t=max(0,it)
       it=min(255,it)
       it1=it/16
       it2=it-16*it1
       ipx(7)=it1
       ipx(8)=it2
c
       write(10,20) ipx
c  20  format("0x",8z1,1x)
   20  format(" 0x",8z1)
c
  100  continue
       return
       end
c      ---------------------
c      subroutine pixeen
c      ---------------------
       subroutine pixeen(xb,yb,xl,yl,icc,zcc)
       x1=xb
       x2=xb+xl
       y1=yb
       y2=yb+yl
c
       write(10,21)
   21  format(" }}")
       write(10,22)
   22  format("geometry IndexedFaceSet {")
       write(10,23)
   23  format("coord Coordinate { point [")
c      write(10,24) nx,ny,ipx0
c  24  format("image ",2i5,i2)
c
       if (icc.eq.2) go to 52
       if (icc.eq.3) go to 53
c
       z=zcc
       write(10,31) x1,y1,z
       write(10,31) x2,y1,z
       write(10,31) x2,y2,z
       write(10,32) x1,y2,z
   31  format(3f6.2,",")
   32  format(3f6.2)
       go to 60
c
   52  continue
       x=zcc
       write(10,31) x,x1,y1
       write(10,31) x,x2,y1
       write(10,31) x,x2,y2
       write(10,32) x,x1,y2
       go to 60
c
   53  continue
       y=zcc
       write(10,31) x1,y,y1
       write(10,31) x2,y,y1
       write(10,31) x2,y,y2
       write(10,32) x1,y,y2
   60  continue
c
       write(10,41)
   41  format(" ]}")
       write(10,42)
   42  format("coordIndex [ 0 1 2 3 ]")
       write(10,43)
   43  format("solid FALSE"/"texCoord TextureCoordinate {")
       write(10,44)
   44  format("point ["/"0.0 0.0,"/"1.0 0.0,"/"1.0 1.0,")
       write(10,45)
   45  format("0.0 1.0"/"]}"/"texCoordIndex [ 0 1 2 3 ]")
       write(10,46)
   46  format("}}")
c
       return
       end