c
C FILE NAME frustrum.for ----------------------------- 180299
c This file creates shapes which form sectors of truncated
c thick-walled cones, the configurations of which are
c defined by the following parameters:
c boufr = bottom outside radius
c binfr = bottom inside radius
c toufr = top outside radius
c tinfr = top inside radius
c nfside = number of facets on the full conical surface
c angfac = fraction of 2 pi occupied by actual conical surface
c angf = angle at which the actual surface starts
c zpf = first z location, ie that of bottom
c zpl = last z location, ie that of top
c isign = 1 for opaque and -1 for transparent objects
C--------------------------------------------------------------
real xp(100000),yp(100000),zp(100000)
integer ipf(4,100000)
character*8 name
c
c...................default data input
c
boufr=1.0
binfr=0.7
toufr=0.49
tinfr=0.4
nfside=10
angfac=0.1
angf=0.1
zpf=0.2
zpl=0.4
name='frustrum'
isign=-1
c.... logical units for infile and outfile
luin=7
open(luin,file='frustrum')
c.... read data from shapdat file
read(luin,'(a)',err=10) name
write(*,*) name
read(luin,*,err=10) boufr
read(luin,*,err=10) binfr
read(luin,*,err=10) toufr
read(luin,*,err=10) tinfr
read(luin,*,err=10) nfside
read(luin,*,err=10) angfac
read(luin,*,err=10) angf
read(luin,*,err=10) zpf
read(luin,*,err=10) zpl
read(luin,*,err=10) isign
close(luin)
go to 11
10 write(6,*)'unable to read frustrum file. defaults taken'
c.... echo data to the screen
11 write(6,*)' name = ', name
write(6,*)' boufr = ', boufr
write(6,*)' binfr = ', binfr
write(6,*)' toufr = ', toufr
write(6,*)' tinfr = ', tinfr
write(6,*)' nfside= ', nfside
write(6,*)' angfac= ', angfac
write(6,*)' angf = ', angf
write(6,*)' zpf = ', zpf
write(6,*)' zpl = ', zpl
write(6,*)' isign = ', isign
C......................................... start calculation
C
pi=3.14159
c.................the angle subtended by a single side facet
ang=angfac * 2.*pi/nfside
C
C cartesian coordinates of points defining the bottom inner ring, of
c which, because of truncation, the radius is as on the next line.
binf=binfr + (tinfr - binfr) * zpf
bouf=boufr + (toufr - boufr) * zpf
tinf=binfr + (tinfr - binfr) * zpl
touf=boufr + (toufr - boufr) * zpl
c......................... do loop starts at 0 not 1
np=0
radius=binf
do i=0,nfside
np=np+1
angle= angf + ang*i
xp(np)=0.5 * ( 1.0 + radius*cos(angle) )
yp(np)=0.5 * ( 1.0 + radius*sin(angle) )
zp(np)=zpf
enddo
C similar coding for the bottom outer ring
radius=bouf
do i=0,nfside
np=np+1
angle= angf + ang*i
xp(np)=0.5 * ( 1.0 + radius*cos(angle) )
yp(np)=0.5 * ( 1.0 + radius*sin(angle) )
zp(np)=zpf
enddo
C Points at top inner rings
radius=tinf
do i=0,nfside
np=np+1
angle= angf + ang*i
xp(np)=0.5 * ( 1.0 + radius*cos(angle) )
yp(np)=0.5 * ( 1.0 + radius*sin(angle) )
zp(np)=zpl
enddo
C Points at top outer rings
radius=touf
do i=0,nfside
np=np+1
angle= angf + ang*i
xp(np)=0.5 * ( 1.0 + radius*cos(angle) )
yp(np)=0.5 * ( 1.0 + radius*sin(angle) )
zp(np)=zpl
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
C
c.... create and write the geometry file into the local directory
open(53,file=name//'.dat',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),isign*( 20+mod(i,200))
40 continue
c.... the following line indicates that the object can have any of
c the 24 possible orientations within the bounding box
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