c
C FILE NAME chamber.for ---------------------------------------------
C---------------------------------------------------------------------------
c
c 14 name.dat files are to be written, corresponding to the 14 objects
c of which the chamber is constructed.
c
parameter (pi=3.14159)
real bin(14),tin(14),angf(14),angl(14),zf(14),zl(14)
dimension icolor(14)
character*12 fname(14)
logical adjoutall
c Here the file names are defined.
data fname/'sol1t.dat', 'sol2t.dat', 'sol3t.dat', 'sol4t.dat',
1 's5cr1t.dat','s5cr2t.dat','s5cr3t.dat','s5cr4t.dat',
1 'gas1t.dat', 'gas2t.dat', 'gas3t.dat', 'gas4t.dat',
1 'sol6t.dat', 'sol7t.dat'/
c
c Here are defined bin and tin, the bottom and top inner radii of the
c cone for each of the 14 objects.
data bin /.2, .5, .8, .9, .91, .92, .9, .9,
1 .9, .9, .9, .9, .9, .89 /
data tin /.5, .8, .8, .91, .92, .9, .9, .9,
1 .9, .9, .9, .9, .89, .5 /
c
c Here are defined the first and last angles subtended by the objects at
c the chamber axis.
data angf /0., 0., 0., 0., .00, .50, 1.00, 1.50,
1 .40, .90, 1.40, 1.90, 0., 0. /
data angl /2., 2., 2., 2., .40, 0.90, 1.40, 1.90,
1 .50, 1.00, 1.50, 2.00, 2., 2. /
c
c Here are defined the first and last axial-coordinate values of the
c objects.
data zf /0., .09, .2, .25, .45, .45, .45, .45,
1 .45, .45, .45, .45, .55,.8 /
data zl /.09, .2, .25, .45, .55, .55, .55, .55,
1 .55, .55, .55, .55, .8, 1./
c Each object is given a different colour
data icolor/ 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
1 110, 120, 130, 140/
c
c The objects will be represented as "transparent"
isign=-1
c The number of facets per side is set to 32
nfside=32
nx=40
nxin=( .5/tan(angl(5)*pi))*nx
write(*,*) 'nxin=',nxin
angnew=0.5-asin(float(nxin)/float(nx/2)/bin(5))/pi
write(*,*) 'angnew=',angnew
do i=5,8
angl(i)=angnew+(i-5)*.5
enddo
do i=9,12
angf(i)=angl(i-4)
enddo
do i=1,14
adjoutall=i.ge.9.and.i.le.12
call make_square_cone(pi, bin(i), tin(i), nfside,
1 angf(i)*pi, angl(i)*pi, zf(i), zl(i),.true.,adjoutall)
call write_geometry
1 ('\phoenics\d_satell\d_vrgeom\myshape\'//fname(i),icolor(i),
1 isign)
enddo
end
C
c The following subroutine creates objects consisting of what is left
c when a symmetrically-placed conical cavity is carved out from a thick
c square slab, whereafter then all remaining material is removed which
c subtends an angle at the symmetry axis between:
c angfst & anglst.
c
c The slab has lower and upper surfaces at:
c z fst & zlst
c
subroutine make_square_cone(pi, binfr, tinfr, nfside,
1 angfst, anglst, zfst, zlst, adjoutxy,adjoutall)
logical adjoutxy,adjoutall,passcorn1,passcorn2,passcorn3,passcorn4
common /myshi1/np,nf,ipf(4,100000)
common /myshr1/xp(100000),yp(100000),zp(100000)
angfac=(anglst-angfst)/(2.*pi)
C
ang=angfac * 2.*pi/nfside
C
np=0
C Points at bottom inner rings
binf=binfr
do i=0,nfside
np=np+1
xp(np)=0.5+0.5*binf*cos(angfst + ang*i)
yp(np)=0.5+0.5*binf*sin(angfst + ang*i)
zp(np)=zfst
enddo
C Points at bottom outer edges
piby4=pi/4.
passcorn1=.false.
passcorn2=.false.
passcorn3=.false.
passcorn4=.false.
do i=0,nfside
np=np+1
angle=angfst + ang * i
if(angle.le.piby4.or.angle.ge.7.*piby4) then
xp(np)=1.
yp(np)=0.5+0.5*tan(angle)
if(angle.ge.7.*piby4.and..not.passcorn1) then
yp(np)=0.
passcorn1=.true.
endif
if(adjoutall.or.(adjoutxy.and.(i.eq.0.or.i.eq.nfside)))
1 yp(np)=yp(np-nfside-1)
elseif(angle.ge.3.*piby4.and.angle.le.5.*piby4) then
xp(np)=0.
yp(np)=0.5-0.5*tan(angle)
if(angle.ge.3.*piby4.and..not.passcorn2) then
yp(np)=1.
passcorn2=.true.
endif
if(adjoutall.or.(adjoutxy.and.(i.eq.0.or.i.eq.nfside)))
1 yp(np)=yp(np-nfside-1)
elseif(angle.ge.piby4.and.angle.le.3.*piby4) then
xp(np)=0.5-0.5*tan(angle-pi/2)
yp(np)=1.
if(angle.ge.piby4.and..not.passcorn3) then
xp(np)=1.
passcorn3=.true.
endif
if(adjoutall.or.(adjoutxy.and.(i.eq.0.or.i.eq.nfside)))
1 xp(np)=xp(np-nfside-1)
else
xp(np)=0.5+0.5*tan(angle-pi/2)
yp(np)=0.
if(angle.ge.5.*piby4.and..not.passcorn4) then
xp(np)=0.
passcorn4=.true.
endif
if(adjoutall.or.(adjoutxy.and.(i.eq.0.or.i.eq.nfside)))
1 xp(np)=xp(np-nfside-1)
endif
zp(np)=zfst
enddo
C Points at top inner rings
tinf=tinfr
do i=0,nfside
np=np+1
xp(np)=0.5+0.5*tinf*cos(angfst + ang*i)
yp(np)=0.5+0.5*tinf*sin(angfst + ang*i)
zp(np)=zlst
enddo
C Points at top outer edges
do i=0,nfside
np=np+1
xp(np)=xp(np-nfside*2-2)
yp(np)=yp(np-nfside*2-2)
zp(np)=zlst
enddo
C facets
nf=0
do i=1,nfside
c... bottom
nf=nf+1
ipf(1,nf)=i
ipf(2,nf)=i+1
ipf(3,nf)=i+nfside+2
ipf(4,nf)=i+nfside+1
c... outside
nf=nf+1
ipf(1,nf)=i+nfside+1
ipf(2,nf)=i+nfside+2
ipf(3,nf)=i+3*nfside+4
ipf(4,nf)=i+3*nfside+3
c... top
nf=nf+1
ipf(1,nf)=i+3*nfside+3
ipf(2,nf)=i+3*nfside+4
ipf(3,nf)=i+2*nfside+3
ipf(4,nf)=i+2*nfside+2
c... inside
nf=nf+1
ipf(1,nf)=i+2*nfside+2
ipf(2,nf)=i+2*nfside+3
ipf(3,nf)=i+1
ipf(4,nf)=i
enddo
if(angfac.lt..999) then
c... close the lower angle end
nf=nf+1
ipf(1,nf)=1
ipf(2,nf)=1+nfside+1
ipf(3,nf)=1+3*nfside+3
ipf(4,nf)=1+2*nfside+2
c... close the higher angle end
nf=nf+1
ipf(1,nf)=nfside+1
ipf(2,nf)=nfside+1+2*nfside+2
ipf(3,nf)=nfside+1+3*nfside+3
ipf(4,nf)=nfside+1+nfside+1
endif
end
C
subroutine write_geometry(filename,icol,isign)
common /myshi1/np,nf,ipf(4,100000)
common /myshr1/xp(100000),yp(100000),zp(100000)
character*(*) filename
open(53,file=filename,status='unknown')
write(53,'(i6)') np
do 30 i=1,np
write(53,'(1p3e12.5)') xp(i),yp(i),zp(i)
30 continue
write(53,'(i6)') nf
do 40 i=1,nf
write(53,'(5i6)') (ipf(j,i),j=1,4), icol * isign
40 continue
write(53,'(a)') ' 24 1 2 3 4 5 6 7 8 9 10 11 12'//
1 ' 13 14 15 16 17 18 19 20 21 22 23 24'
close(53)
end
c