program triplot
type surface
integer nsurface
real*4,dimension(:),pointer::r,s,x,y,z,xn,yn,zn
end type surface
integer window(4),red(230),green(230),blue(230)
integer xchoice
external xchoice
real*4 xr(8),yr(8),zr(8),tr(8),psi(8)
real*4 ur(8),vr(8),wr(8)
real*4 xc(7),yc(7),zc(7),tc(7)
real*4 uc(7),vc(7),wc(7)
real*4,dimension(:,:),allocatable::lsfc,lsfr
character fnme*256
real*4,dimension(:),allocatable::x,y,z,u,v,w,t,p,s
real*4,dimension(:,:),allocatable::lsf
integer,dimension(:,:),allocatable::icon,iface
integer,dimension(:),allocatable::octree,leaf,node,ftr,rtf,face
integer,dimension(:),allocatable::kfixx,kfixy,kfixz,kfixt
logical,dimension(:),allocatable::influid
type (surface),dimension(:),allocatable::surf
data xr/0.,1.,0.,1.,0.,1.,0.,1./
data yr/0.,0.,1.,1.,0.,0.,1.,1./
data zr/0.,0.,0.,0.,1.,1.,1.,1./
window(1)=100
window(2)=100
window(3)=700
window(4)=700
call xopen (window,char(0))
red=0
green=0
blue=0
red(1)=180
green(1)=120
green(3)=255
ncol=20
do i=1,ncol
red(100+i)=int(255.*float(i)/ncol)
blue(100+i)=255-int(255.*float(i)/ncol)
enddo
call xcmap (red,green,blue)
call xscale (0,window(3),0,window(4),-0.05,1.05,1.05,-0.05,1)
! call xscale (0,window(3),0,window(4),-0.01,0.01,0.01,-0.01,1)
! call xscale (0,window(3),0,window(4),-0.05,.25,.85,.55,1)
! call xscale (0,window(3),0,window(4),.75,1.05,.85,.55,1)
call xpen (1)
open (7,file='triplot.txt',status='old')
read (7,*) xn,yn,zn,d
read (7,'(a)') fnme
nfnme=len_trim(fnme)
close (7)
xyzn=sqrt(xn**2+yn**2+zn**2)
xn=xn/xyzn
yn=yn/xyzn
zn=zn/xyzn
d=d/xyzn
open (7,file=fnme(1:nfnme),status='old')
read (7,*) noctree,nnode,nleaves,nface,nlsf
allocate (x(nnode),y(nnode),z(nnode),u(nnode),v(nnode),w(nnode))
allocate (t(nnode),s(nnode),p(nnode),lsf(nnode,nlsf),lsfr(8,nlsf),lsfc(7,nlsf))
allocate (kfixt(nnode),kfixx(nnode),kfixy(nnode),kfixz(nnode))
allocate (surf(nlsf))
do k=1,nnode
read (7,*) x(k),y(k),z(k),u(k),v(k),w(k), &
(lsf(k,j),j=1,nlsf),t(k),s(k),p(k),kfixx(k), &
kfixy(k),kfixz(k),kfixt(k)
enddo
allocate (icon(8,nleaves))
read (7,*) ((icon(k,i),k=1,8),i=1,nleaves)
allocate (octree(noctree))
read (7,*) (octree(k),k=1,noctree)
allocate (iface(9,nface))
do i=1,nface
read (7,*) (iface(k,i),k=1,9)
enddo
allocate (node(nnode),leaf(nnode),ftr(nnode))
allocate (rtf(nnode),influid(nnode))
do k=1,nnode
read (7,*) node(k),leaf(k),ftr(k),rtf(k),influid(k)
enddo
allocate (face(nface))
do k=1,nface
read (7,*) face(k)
enddo
do k=1,nlsf
read (7,*) nsurface
surf(k)%nsurface=nsurface
allocate (surf(k)%x(nsurface),surf(k)%y(nsurface),surf(k)%z(nsurface), &
surf(k)%xn(nsurface),surf(k)%yn(nsurface),surf(k)%zn(nsurface), &
surf(k)%r(nsurface),surf(k)%s(nsurface))
do i=1,nsurface
read (7,*) rs,ss,xs,ys,zs,xns,yns,zns
surf(k)%x(i)=xs
surf(k)%y(i)=ys
surf(k)%z(i)=zs
surf(k)%xn(i)=xns
surf(k)%yn(i)=yns
surf(k)%zn(i)=zns
surf(k)%r(i)=rs
surf(k)%s(i)=ss
enddo
enddo
close (7)
do k=1,nlsf
print*,'Surface ',k,surf(k)%nsurface
print*,minval(surf(k)%x(:)),maxval(surf(k)%x(:))
print*,minval(surf(k)%y(:)),maxval(surf(k)%y(:))
print*,minval(surf(k)%z(:)),maxval(surf(k)%z(:))
print*,minval(surf(k)%xn(:)),maxval(surf(k)%xn(:))
print*,minval(surf(k)%yn(:)),maxval(surf(k)%yn(:))
print*,minval(surf(k)%zn(:)),maxval(surf(k)%zn(:))
print*,minval(surf(k)%r(:)),maxval(surf(k)%r(:))
print*,minval(surf(k)%s(:)),maxval(surf(k)%s(:))
enddo
fac=maxval(sqrt(u**2+v**2+w**2))*20.
do k=1,nnode
! if (abs(z(k)-.8125).lt.1.e-6) print*,x(k),y(k),z(k),lsf(k,1),influid(k),t(k)
enddo
do ileave=1,nleaves
do k=1,8
xr(k)=x(icon(k,ileave))
yr(k)=y(icon(k,ileave))
zr(k)=z(icon(k,ileave))
tr(k)=t(icon(k,ileave))
ur(k)=u(icon(k,ileave))
vr(k)=v(icon(k,ileave))
wr(k)=w(icon(k,ileave))
do j=1,nlsf
lsfr(k,j)=lsf(icon(k,ileave),j)
enddo
psi(k)=xn*xr(k)+yn*yr(k)+zn*zr(k)-d
enddo
ncross=0
do k=1,2
do j=1,2
do i=1,2
ijk=i+(j-1)*2+(k-1)*4
if (i.ne.2) then
ijkp=(i+1)+(j-1)*2+(k-1)*4
if (psi(ijk)*psi(ijkp).le.0) then
ncross=ncross+1
rat=abs(psi(ijk))/(abs(psi(ijk))+abs(psi(ijkp)))
xc(ncross)=xr(ijk)+(xr(ijkp)-xr(ijk))*rat
yc(ncross)=yr(ijk)+(yr(ijkp)-yr(ijk))*rat
zc(ncross)=zr(ijk)+(zr(ijkp)-zr(ijk))*rat
tc(ncross)=tr(ijk)+(tr(ijkp)-tr(ijk))*rat
uc(ncross)=ur(ijk)+(ur(ijkp)-ur(ijk))*rat
vc(ncross)=vr(ijk)+(vr(ijkp)-vr(ijk))*rat
wc(ncross)=wr(ijk)+(wr(ijkp)-wr(ijk))*rat
do l=1,nlsf
lsfc(ncross,l)=lsfr(ijk,l)+(lsfr(ijkp,l)-lsfr(ijk,l))*rat
enddo
endif
endif
if (j.ne.2) then
ijkp=i+j*2+(k-1)*4
if (psi(ijk)*psi(ijkp).le.0) then
ncross=ncross+1
xc(ncross)=xr(ijk)+(xr(ijkp)-xr(ijk))*rat
yc(ncross)=yr(ijk)+(yr(ijkp)-yr(ijk))*rat
zc(ncross)=zr(ijk)+(zr(ijkp)-zr(ijk))*rat
tc(ncross)=tr(ijk)+(tr(ijkp)-tr(ijk))*rat
uc(ncross)=ur(ijk)+(ur(ijkp)-ur(ijk))*rat
vc(ncross)=vr(ijk)+(vr(ijkp)-vr(ijk))*rat
wc(ncross)=wr(ijk)+(wr(ijkp)-wr(ijk))*rat
do l=1,nlsf
lsfc(ncross,l)=lsfr(ijk,l)+(lsfr(ijkp,l)-lsfr(ijk,l))*rat
enddo
endif
endif
if (k.ne.2) then
ijkp=i+(j-1)*2+k*4
if (psi(ijk)*psi(ijkp).le.0) then
ncross=ncross+1
xc(ncross)=xr(ijk)+(xr(ijkp)-xr(ijk))*rat
yc(ncross)=yr(ijk)+(yr(ijkp)-yr(ijk))*rat
zc(ncross)=zr(ijk)+(zr(ijkp)-zr(ijk))*rat
tc(ncross)=tr(ijk)+(tr(ijkp)-tr(ijk))*rat
uc(ncross)=ur(ijk)+(ur(ijkp)-ur(ijk))*rat
vc(ncross)=vr(ijk)+(vr(ijkp)-vr(ijk))*rat
wc(ncross)=wr(ijk)+(wr(ijkp)-wr(ijk))*rat
do l=1,nlsf
lsfc(ncross,l)=lsfr(ijk,l)+(lsfr(ijkp,l)-lsfr(ijk,l))*rat
enddo
endif
endif
enddo
enddo
enddo
if (ncross.gt.0) then
call order (xc,yc,zc,tc,uc,vc,wc,lsfc,nlsf,ncross,xn,yn,zn)
xc(ncross+1)=xc(1)
yc(ncross+1)=yc(1)
zc(ncross+1)=zc(1)
tc(ncross+1)=tc(1)
uc(ncross+1)=uc(1)
vc(ncross+1)=vc(1)
wc(ncross+1)=wc(1)
do k=1,nlsf
lsfc(ncross+1,k)=lsfc(1,k)
enddo
do k=1,ncol
val1=float(k-1)/float(ncol)
if (k.eq.1) val1=val1+1.e-5
val2=float(k)/float(ncol)
if (k.eq.ncol) val2=val2+1.e-5
icol=100+k
! call cont (xc,yc,zc,tc,ncross,val1,val2,icol,xn,yn,zn)
! call contl (xc,yc,zc,tc,ncross,val1,1,xn,yn,zn)
enddo
call xpen (3)
call project (xc(ncross),yc(ncross),zc(ncross),xn,yn,zn,xcp,ycp)
call xplot (xcp,ycp,3)
do k=1,ncross
call project (xc(k),yc(k),zc(k),xn,yn,zn,xcp,ycp)
call xplot (xcp,ycp,2)
enddo
call xpen (4)
rad=1./4000.
rad=1./400.
do k=1,8
if (kfixt(icon(k,ileave)).eq.1) then
call project (xr(k),yr(k),zr(k),xn,yn,zn,xrp,yrp)
! call xcircle (xrp,yrp,rad)
endif
enddo
call xpen (5)
xm=sum(xc(1:ncross))/ncross
ym=sum(yc(1:ncross))/ncross
zm=sum(zc(1:ncross))/ncross
um=sum(uc(1:ncross))/ncross
vm=sum(vc(1:ncross))/ncross
wm=sum(wc(1:ncross))/ncross
call project (xm,ym,zm,xn,yn,zn,xmp,ymp)
call project (xm+um/fac,ym+vm/fac,zm+wm/fac,xn,yn,zn,xmpp,ympp)
call arrow (xmp,ymp,xmpp,ympp)
do k=1,nlsf
call contl (xc,yc,zc,lsfc(1,k),ncross,0.0001,2,xn,yn,zn)
enddo
endif
enddo
call xpen (6)
do k=1,nlsf
do i=1,surf(k)%nsurface
call project (surf(k)%x(i),surf(k)%y(i),surf(k)%z(i),xn,yn,zn,xmp,ymp)
! if (k.eq.2) call xcirclef (xmp,ymp,rad)
! if (k.eq.2) call xcirclef (surf(k)%r(i),surf(k)%s(i),rad)
enddo
enddo
call xbox (xn,yn,zn)
call xsave ()
ichoice=xchoice ('Continue ',1)
call xclose ()
end
! --------------------------
subroutine order (x,y,z,t,u,v,w,lsf,nlsf,n,xn,yn,zn)
real*4,dimension(:),allocatable::xp,yp,zp,tp,up,vp,wp,an
real*4,dimension(:,:),allocatable::lsfp
integer,dimension(:),allocatable::indx
real*4 x(7),y(7),z(7),t(7),u(7),v(7),w(7),lsf(7,nlsf)
allocate (xp(n),yp(n),zp(n),tp(n),up(n),vp(n),wp(n), &
lsfp(n,nlsf),an(n),indx(n))
xp=x(1:n)
yp=y(1:n)
zp=z(1:n)
tp=t(1:n)
up=u(1:n)
vp=v(1:n)
wp=w(1:n)
lsfp=lsf(1:n,:)
xc=sum(x(1:n))/n
yc=sum(y(1:n))/n
zc=sum(z(1:n))/n
do i=1,n
an(i)=atan2(((y(i)-yc)*(z(1)-zc)-(z(i)-zc)*(y(1)-yc))*xn+ &
((z(i)-zc)*(x(1)-xc)-(x(i)-xc)*(z(1)-zc))*yn+ &
((x(i)-xc)*(y(1)-yc)-(y(i)-yc)*(x(1)-xc))*zn, &
(x(i)-xc)*(x(1)-xc)+ &
(y(i)-yc)*(y(1)-yc)+ &
(z(i)-zc)*(z(1)-zc))
enddo
call indexx (n,an,indx)
do i=1,n
x(i)=xp(indx(i))
y(i)=yp(indx(i))
z(i)=zp(indx(i))
t(i)=tp(indx(i))
u(i)=up(indx(i))
v(i)=vp(indx(i))
w(i)=wp(indx(i))
do k=1,nlsf
lsf(i,k)=lsfp(indx(i),k)
enddo
enddo
deallocate (xp,yp,zp,tp,up,vp,wp,lsfp,an,indx)
return
end
! -------------------------
SUBROUTINE INDEXX(N,ARRIN,INDX)
DIMENSION ARRIN(N),INDX(N)
DO 11 J=1,N
INDX(J)=J
11 CONTINUE
L=N/2+1
IR=N
10 CONTINUE
IF(L.GT.1)THEN
L=L-1
INDXT=INDX(L)
Q=ARRIN(INDXT)
ELSE
INDXT=INDX(IR)
Q=ARRIN(INDXT)
INDX(IR)=INDX(1)
IR=IR-1
IF(IR.EQ.1)THEN
INDX(1)=INDXT
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF(J.LE.IR)THEN
IF(J.LT.IR)THEN
IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
ENDIF
IF(Q.LT.ARRIN(INDX(J)))THEN
INDX(I)=INDX(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
INDX(I)=INDXT
GO TO 10
END
! ---------------------------------
subroutine cont (x,y,z,t,n,v1,v2,ic,xn,yn,zn)
real x(n+1),y(n+1),z(n+1),t(n+1)
real xpp(10),ypp(10),zpp(10),xppp(10),yppp(10)
nplot=0
do k=1,n
if ((t(k)-v1)*(t(k)-v2).le.0.) then
nplot=nplot+1
xpp(nplot)=x(k)
ypp(nplot)=y(k)
zpp(nplot)=z(k)
endif
ratio1=(v1-t(k))/(t(k+1)-t(k))
ratio2=(v2-t(k))/(t(k+1)-t(k))
ratio=min(ratio1,ratio2)
if (ratio*(ratio-1.).le.0.) then
nplot=nplot+1
xpp(nplot)=x(k)+ratio*(x(k+1)-x(k))
ypp(nplot)=y(k)+ratio*(y(k+1)-y(k))
zpp(nplot)=z(k)+ratio*(z(k+1)-z(k))
endif
ratio=max(ratio1,ratio2)
if (ratio*(ratio-1.).le.0.) then
nplot=nplot+1
xpp(nplot)=x(k)+ratio*(x(k+1)-x(k))
ypp(nplot)=y(k)+ratio*(y(k+1)-y(k))
zpp(nplot)=z(k)+ratio*(z(k+1)-z(k))
endif
enddo
if (nplot.ne.0) then
call xpen(ic)
do k=1,nplot
call project (xpp(k),ypp(k),zpp(k),xn,yn,zn,xppp(k),yppp(k))
enddo
call xfill (xppp,yppp,nplot)
endif
return
end
! -----------------------------
subroutine contl (x,y,z,t,n,v,ic,xn,yn,zn)
real x(n+1),y(n+1),z(n+1),t(n+1)
ipen=1
call xpen (ic)
do k=1,n
if ((t(k)-v)*(t(k+1)-v).le.0.) then
ipen=-ipen
rat=.5
if (t(k+1).ne.t(k)) rat=(v-t(k))/(t(k+1)-t(k))
xp=x(k)+rat*(x(k+1)-x(k))
yp=y(k)+rat*(y(k+1)-y(k))
zp=z(k)+rat*(z(k+1)-z(k))
call project (xp,yp,zp,xn,yn,zn,xpp,ypp)
if (ipen.gt.0) call xplot (xpp,ypp,2)
if (ipen.lt.0) call xplot (xpp,ypp,3)
endif
enddo
return
end
! -------------------------------
subroutine xbox (xn,yn,zn)
call xpen (1)
call project (0.,0.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,3)
call project (1.,0.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (1.,1.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (0.,1.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (0.,0.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (0.,0.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (1.,0.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (1.,1.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (0.,1.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (0.,0.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (1.,0.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,3)
call project (1.,0.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (1.,1.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,3)
call project (1.,1.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
call project (0.,1.,0.,xn,yn,zn,xp,yp)
call xplot (xp,yp,3)
call project (0.,1.,1.,xn,yn,zn,xp,yp)
call xplot (xp,yp,2)
return
end
! -----------------------------
subroutine project (x,y,z,xn,yn,zn,xp,yp)
if (abs(xn).ge.abs(yn) .and. abs(xn).ge.abs(zn)) then
xp=-yn*x+xn*y
yp=-zn*x+xn*z
return
elseif (abs(yn).ge.abs(zn)) then
xp=yn*x-xn*y
yp=-zn*y+yn*z
return
else
xp=zn*x-xn*z
yp=zn*y-yn*z
return
endif
end
!-----
subroutine arrow (x1,y1,x2,y2)
call xplot (x1,y1,3)
call xplot (x2,y2,2)
if (y2-y1.eq.0. .and. x2-x1.eq.0.) then
ang=0.
else
ang=atan2(y2-y1,x2-x1)
endif
ang1=ang+3.141592654+.5
ang2=ang+3.141592654-.5
dl=sqrt((x2-x1)**2+(y2-y1)**2)/3.
call xplot (x2+dl*cos(ang1),y2+dl*sin(ang1),3)
call xplot (x2,y2,2)
call xplot (x2+dl*cos(ang2),y2+dl*sin(ang2),2)
return
end