c----------------------------------------------
c      ebmesh3d.f
c    
c      subrtn file : pk1subrtn
c       
c----------------------------------------------
c***********************************************************
c      main program  ebmesh3d.f                            *
c                   < z houkou 3D >                        *
c                                                          *
c***********************************************************
c      main program  ebmesh3d.f
c      triangm1.f    for triangle mesh program
       parameter (x0=0.0,y0=0.0,z0=0.0)
       parameter (xl=50.1,yl=20.1,zl=20.1)
       parameter (mx=500,mx2=mx+2,my=200,my2=my+2)
       parameter (mz=200,mz2=mz+2)
       parameter (mm1=mx2*my2,mm2=my2*my2)
c----( nx: x point, ny: y point )
c      parameter (nx=202,ny=202)
       parameter (nx=mx2,ny=my2,nz=10)
c      parameter (nx=my2,ny=my2)
c----( p: x,y,z no jun ni set   u: r,g,b no jun ni set )
       dimension pp1(mm1),pp2(mm2)
       dimension p(3*nx*ny)
       dimension u(3*nx*ny)
c
       open(21,file='pp1211380.data',
     1         access='sequential',form='unformatted')
       open(22,file='pp2211380.data',
     1         access='sequential',form='unformatted')
c
c
       read(21) pp1
       read(22) pp2
c
c
       dx=xl/float(mx2-1)
       dy=yl/float(my2-1)
c
       call initvrml
c
c################################( test data )
       pi=3.1415926
       umin=pp1(1)
       umax=pp1(1)
       do 10 i=1,mm1
       umin=amin1(umin,pp1(i))
       umax=amax1(umax,pp1(i))
   10  continue
c
       do 12 i=1,mm1
       pp1(i)=(pp1(i)-umin)/(umax-umin)
   12  continue
c
c
       nys=ny-1
       amp1=zl
c
       do 20  j=1,ny
       do 20  i=1,nx
       j1=i+mx2*(j-1)
       i1=1+3*(i-1+nx*(j-1))
       i2=i1+1
       i3=i1+2
       p(i1)=x0+dx*(i-1)
       p(i2)=y0+dy*(j-1)
       p(i3)=amp1*pp1(j1)
       cu1=4.0*pp1(j1)
c
       r=cu1
       g=1.0-cu1
       b=0.0
c
       cu2=4.0*cu1
       r=cu2-1.0
       g=cu2
       if(cu2.gt.1.0) g=2.0-cu2
       if(cu2.gt.2.0) g=cu2-2.0
       b=1.0-cu2
       if(cu2.gt.3.0) b=cu2-3.0
c
       r=amax1(r,0.0)
       r=amin1(r,1.0)
       g=amax1(g,0.0)
       g=amin1(g,1.0)
       b=amax1(b,0.0)
       b=amin1(b,1.0)
       u(i1)=r
       u(i2)=g
       u(i3)=b
   20  continue
c#############################################
c
c      do  30  k=1,nz
       call triangm1(nx,ny,u,p)
c  30  continue
c
       stop
       end
