diff options
Diffstat (limited to 'Dragon/src/XELTI3.f')
| -rw-r--r-- | Dragon/src/XELTI3.f | 433 |
1 files changed, 433 insertions, 0 deletions
diff --git a/Dragon/src/XELTI3.f b/Dragon/src/XELTI3.f new file mode 100644 index 0000000..9618bed --- /dev/null +++ b/Dragon/src/XELTI3.f @@ -0,0 +1,433 @@ +*DECK XELTI3 + SUBROUTINE XELTI3( IPRT,IFTEMP,NANGLE,DENUSR,ISYMM,ANGLES,DENSTY, + > NTOTCL,NEXTGE,MAXR,REMESH,LINMAX,RCUTOF, + > NSUR,NVOL,INDEL,MINDIM, + > MAXDIM,ICOORD,INCR,ICUR,TRKBEG,CONV,TRKDIR, + > LENGHT,NUMERO,DDENWT) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Construct the sequential tape that will contain tracks for +* isotropic BC in 3-D calculation. +* +*Copyright: +* Copyright (C) 1989 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): R. Roy +* +*Parameters: input +* IPRT intermediate printing level for output. +* IFTEMP tracking file number. +* NANGLE number of angles used in the tracking process. +* DENUSR density of tracks in the plane perpendicular +* to the tracking angles. +* ISYMM flag for symetry (1/0 for on/off): +* 2 reflection plane normal to X axis; +* 4 reflection plane normal to Y axis; +* 8 reflection plane normal to X and Y axis; +* 16 reflection plane normal to Z axis; +* 18 reflection plane normal to X and Z axis; +* 20 reflection plane normal to Y and Z axis; +* 24 reflection plane normal to X, Y and Z axis. +* ANGLES 3d angle values. +* DENSTY density of tracks angle by angle. +* NTOTCL number of cylindres of a type + 3. +* NEXTGE for tubez, nextge=1 +* MAXR max number of real mesh values in REMESH. +* REMESH real mesh values (rect/cyl). +* LINMAX max. number of track segments in a single track. +* RCUTOF cutof for corner tracking( 0.25 suggested ) +* NSUR number of surfaces. +* NVOL number of zones. +* INDEL numbering of surfaces & zones. +* MINDIM min index values for all axes (rect/cyl). +* MAXDIM max index values for all axes (rect/cyl). +* ICOORD principal axes direction (X/Y/Z) for meshes. +* ICUR current zonal location for a track segment. +* INCR increment direction for next track segment. +* TRKBEG position where a track begins. +* CONV segments of tracks. +* TRKDIR direction of a track in all axes. +* LENGHT relative lenght of each segment in a track. +* NUMERO material identification of each track segment. +* DDENWT density of tracks angle by angle. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +* +* DECLARE DUMMY ARGUMENTS + INTEGER IPRT,IFTEMP,NANGLE,NTOTCL,NEXTGE,MAXR,LINMAX, + > NSUR,NVOL + INTEGER MINDIM(NTOTCL), MAXDIM(NTOTCL), ICUR(NTOTCL), + > ICOORD(NTOTCL), INCR(NTOTCL), NUMERO(LINMAX), + > INDEL(4,*) + REAL DENUSR,REMESH(MAXR),TRKBEG(NTOTCL),TRKDIR(NTOTCL), + > CONV(NTOTCL),RCUTOF + DOUBLE PRECISION ANGLES(3,*),DENSTY(*),LENGHT(LINMAX), + > DDENWT(NANGLE) + INTEGER ISYMM + INTEGER IQUART(4) +* +* DECLARE LOCAL VARIABLES + REAL TRKEND(3), TRKORI(3), TRKOR2(3), PROJC2(3), + > ANGEQN(3,3), ANGLE2(3), ANGLE3(3), BARY(3), + > TRKCUT(3,2), TCUTOF(3,4), TORIC(3) + INTEGER NSBEG(4), NSEND(4), NSCUT(2) + LOGICAL LANGLE + EQUIVALENCE ( ANGEQN(1,2), ANGLE2 ), ( ANGEQN(1,3), ANGLE3 ) + CHARACTER TEDATA*13 + INTEGER NPAN, IOUT + PARAMETER ( NPAN=3, IOUT=6 ) +* + INTEGER NDIM, I, J, NPOINT, NPO2, NCUTOF, + > NOTRAK, NSOLAN, IANG, ISB, + > NANGLS, IPAN, NESTIM, JANG, IX, IY, IREF1, I2, I3, + > NSOLMX, NTTRK, IANG0, NDEBS, N, IZZ, IANGL, K, K3, + > NCROS, LINE, JC, ISE, IDIM + REAL ANORM2, A, B, DENLIN, RCIRC, R2CIRC, DP, + > SURTOT, VOLTOT, DEPART, X, Y, ANN, ODDNXT, + > TOTLEN, TOTXXX, DDENST + REAL DZ,ZMAX + DOUBLE PRECISION WEIGHT, WZ +* + ANORM2( A, B ) = A*A + B*B + NDIM= 3 +* +* ONE WEIGHT FOR ALL TRACKS +* DENLIN= # OF TRACKS / CM + DENLIN= SQRT(DENUSR) +* +* COMPUTE THE CIRCUMSCRIBED RADIUS AND +* THE COORDINATES FOR THE TRUE CENTER OF THE CELL + R2CIRC= 0.0 + DO 10 I = 1, 3 + BARY(I)= 0.5 * (REMESH(MAXDIM(I)) + REMESH(MINDIM(I))) + IF( NEXTGE.EQ.1 )THEN + CALL XABORT('XELTI3: TUBEZ NOT SUPPORTED') + ELSE + R2CIRC = R2CIRC + > + (REMESH(MAXDIM(I)) - REMESH(MINDIM(I)))**2 + ENDIF + 10 CONTINUE + ZMAX=MAX(ABS(REMESH(MAXDIM(3))),ABS(REMESH(MINDIM(3)))) + R2CIRC= 0.25 * R2CIRC + RCIRC = SQRT(R2CIRC) +* +* NPOINT= # OF TRACKS ALONG ONE PERPENDICULAR AXIS + NPOINT= INT( 2. * RCIRC * DENLIN ) +***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** +***** CHANGE THIS "NPOINT" PARAMETER HAS TREMENDOUS EFFECTS ON TRACKING +***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** BEWARE ***** +* +* OTHER POSSIBLE CHOICES (EXPLORED WITHOUT SUCCESS) ARE ==> +*1-) NPOINT= INT( 2. * RCIRC * DENLIN ) + 1 +*2-) NPOINT= NINT( 2. * RCIRC * DENLIN ) +*3-) NPOINT= NINT( 2. * RCIRC * DENLIN ) + 1 +* +* KEEP "NPOINT" ODD & CORRECT DENSITY + NPO2 = NPOINT / 2 + NPOINT= 2 * NPO2 + 1 + DP = 2. * RCIRC / NPOINT + DENLIN= 1. / DP + DENUSR= DENLIN**2 + ODDNXT= (2*NPO2+3) / (2.*RCIRC) + IF( RCUTOF.EQ.0.0 )THEN + NCUTOF= 1 + ELSE + NCUTOF= 4 + ENDIF + NOTRAK= 0 + SURTOT= 0.0 + VOLTOT= 0.0 + IQUART(1)=1 + IQUART(2)=1 + IQUART(3)=1 + IQUART(4)=1 + IF(NEXTGE .EQ. 1) THEN + IQUART(2)=0 + IQUART(3)=0 + IQUART(4)=0 + NSOLAN= (NANGLE * (NANGLE + 2)) / 2 + ELSE + IF( ISYMM .EQ. 8 .OR. ISYMM .EQ. 24 )THEN + NSOLAN= (NANGLE * (NANGLE + 2)) / 8 + IQUART(2)=0 + IQUART(3)=0 + IQUART(4)=0 + ELSE IF( ISYMM .EQ. 4 .OR. ISYMM .EQ. 20 )THEN + NSOLAN= (NANGLE * (NANGLE + 2)) / 4 + IQUART(2)=0 + IQUART(4)=0 + ELSE IF( ISYMM .EQ. 2 .OR. ISYMM .EQ. 18 )THEN + NSOLAN= (NANGLE * (NANGLE + 2)) / 4 + IQUART(3)=0 + IQUART(4)=0 + ELSE + NSOLAN= (NANGLE * (NANGLE + 2)) / 2 + ENDIF + ENDIF + DDENST= 1.0/(NPAN*DENUSR) + NANGLS= (NANGLE * (NANGLE + 2)) / 2 + CALL XELEQN( 3, 0, ANGEQN ) + IANG= 0 + DO 15 JANG= 1, NANGLS + DO 16 IPAN= 1, NPAN + CALL XELEQN( 3, NANGLE, ANGEQN ) + 16 CONTINUE + IF(IQUART(MOD(JANG-1,4)+1).NE.1 ) GO TO 15 + IANG= IANG+1 + DENSTY(IANG)= REAL(2*NSOLAN) + ANGLES(1,IANG)= ANGEQN(1,1) + ANGLES(2,IANG)= ANGEQN(2,1) + ANGLES(3,IANG)= ANGEQN(3,1) + 15 CONTINUE +* +* COPY ANGLES AND DENSITIES ON TEMPORARY TRACKING FILE + WRITE(IFTEMP) ((ANGLES(IDIM,IANG),IDIM=1,NDIM),IANG=1,NSOLAN) + WRITE(IFTEMP) (DENSTY(IANG) ,IANG=1,NSOLAN) +* +* TO REINITIATE THE EQN ANGLES + CALL XELEQN( 3, 0, ANGEQN ) + IF( NEXTGE.EQ.1 )THEN + DDENST= 12.0*DDENST + ENDIF + WEIGHT= 0.5*DDENST/DBLE(NSOLAN) + DO IANG= 1, NSOLAN + DDENWT(IANG)=WEIGHT + ENDDO + NSOLMX= 0 + NDEBS= 0 + IF( IPRT.GT.1 )THEN +* +* PREPARE & PRINT THE ESTIMATED NUMBER OF TRACKS + NESTIM= 0 + DEPART= - (NPO2+1) * DP + X = DEPART + DO 25 IX = 1, NPOINT + X = X + DP + Y = DEPART + DO 20 IY = 1, NPOINT + Y = Y + DP + IF( ANORM2( X, Y ) .LE. R2CIRC ) NESTIM= NESTIM + 1 + 20 CONTINUE + 25 CONTINUE + WRITE(IOUT,'(1H )') + WRITE(IOUT,'( 8H0ECHO = ,I8,20H TRACKS/AXIS/ANGLE )') + > NPOINT + WRITE(IOUT,'( 8H ECHO = ,I8,25H TRACKS/CIRCLE/ANGLE )') + > NESTIM + NESTIM= NESTIM * NPAN * NSOLAN + WRITE(IOUT,'(30H ECHO = NEXT ODD DENSITY > ,F15.7,4H/CM2)') + > ODDNXT**2 + WRITE(IOUT,'( 8H ECHO = ,28H ESTIMATED NUMBER OF TRACKS= ,I8 )') + > NESTIM +* +* PREPARE PRINTING WITH VARIABLE FORMAT + WRITE(IOUT,'(1H )') + WRITE(IOUT,'( 8H0ECHO = ,I3,27H SOLID ANGLES TO BE TRACKED )') + > NANGLS + NSOLMX= MIN(9, NANGLS/10) + IREF1 = 0 + WRITE(IOUT,'( 1H0,10(I1,9X))') (IREF1, IZZ=0,NSOLMX) + WRITE(IOUT,'( 1H ,10(I1,9X))') (MOD(IZZ,10), IZZ=0,NSOLMX) + WRITE(IOUT,'( 2H 0)') + TEDATA= '(1H+,TXXX,I1)' + ENDIF + IANG = 0 + IANG0 = 0 + NTTRK = 0 + DO 290 IANGL= 1, NANGLS + IANG=IANG+1 + IF(IQUART(MOD(IANG-1,4)+1).NE.1)THEN +*---- +* Do not track this angle because of the problem symmetry +*---- + LANGLE= .FALSE. + ELSE +*---- +* Track this angle +*---- + IANG0= IANG0+1 + LANGLE=.TRUE. + ENDIF + IF( IPRT.GT.1) THEN + IF( MOD(IANGL,100) .EQ. 0 )THEN + IREF1=IREF1+1 + NDEBS= NSOLMX+1 + NSOLMX=MIN(NDEBS+9, NANGLS/10) + WRITE(IOUT,'( 1H0,10(I1,9X))')(IREF1,IZZ=NDEBS,NSOLMX) + WRITE(IOUT,'( 1H ,10(I1,9X))') + > (MOD(IZZ,10),IZZ=NDEBS,NSOLMX) + WRITE(IOUT,'( 2H 0)') + ELSE + IF( IPRT.GT.10000.AND.MOD(IANGL,100).NE.0 )THEN + WRITE(IOUT,'( 1H ,10(I1,9X))')(IREF1,IZZ=NDEBS,NSOLMX) + WRITE(IOUT,'( 1H ,10(I1,9X))') + > (MOD(IZZ,10),IZZ=NDEBS,NSOLMX) + WRITE(IOUT,'( 2H )') + ENDIF + WRITE(TEDATA(7:9),'(I3.3)') MOD(IANGL,100) + 2 + WRITE(IOUT,TEDATA) MOD(IANGL,10) + ENDIF + ENDIF +* +* NPAN AXES DESCRIPTION (X=0.0, Y=0.0 & Z=0.0) + DO 250 IPAN= 1, NPAN +*---- +* Start tesp print +* WRITE(IOUT,7001) IANGL,IPAN +* 7001 FORMAT(' ANGLE = ',I8,5X,'PLAN =',I3) +* Finish test print +*---- + CALL XELEQN( 3, NANGLE, ANGEQN ) + IF(.NOT.LANGLE) GO TO 250 + IF( NEXTGE.EQ.1 )THEN + IF( IPAN.NE.2 ) GO TO 250 + ENDIF + DO 30 I = 1, 3 + N = ICOORD(I) + TRKDIR(N)= ANGEQN(N,1) + INCR(I)= +1 + IF( TRKDIR(N) .LT. 0.0 ) INCR(I)= -1 +* +* MODIFY ANGLES TO TAKE INTO ACCOUNT DP + ANGLE2(I)= DP * ANGLE2(I) + ANGLE3(I)= DP * ANGLE3(I) + IF( NCUTOF.NE.1 )THEN + TCUTOF(I,1)= RCUTOF*( ANGLE2(I)+ ANGLE3(I) ) + TCUTOF(I,2)= RCUTOF*( ANGLE2(I)- ANGLE3(I) ) + TCUTOF(I,3)= -TCUTOF(I,2) + TCUTOF(I,4)= -TCUTOF(I,1) + ENDIF +* +* DETERMINE THE ORIGIN OF ALL TRACKS + TRKOR2(I)= BARY(I) - (NPO2+1)*(ANGLE2(I)+ANGLE3(I)) + 30 CONTINUE + DO 45 I = 1, 3 + PROJC2(I)= 0.0 + DO 40 J = 1, 3 + IF( I.EQ.J ) GO TO 40 + PROJC2(I)= PROJC2(I) + TRKDIR(J) * TRKDIR(J) + 40 CONTINUE + 45 CONTINUE +* +* SCAN ALL TRACKS IN THE PERPENDICULAR PLANE + DO 180 I2 = 1, NPOINT + DO 50 J = 1, 3 + TRKOR2(J)= TRKOR2(J) + ANGLE2(J) + TRKORI(J)= TRKOR2(J) + 50 CONTINUE + DO 170 I3 = 1, NPOINT + ANN = 0.0 + DO 60 J = 1, 3 + TRKORI(J)= TRKORI(J) + ANGLE3(J) + ANN= ANN + (TRKORI(J)-BARY(J))**2 + 60 CONTINUE +*---- +* Start tesp print +* WRITE(IOUT,7002) I2,I3,(TRKORI(K)-BARY(K),K=1,NDIM) +*7002 FORMAT(' ORIGINE MESH:',I10,5X,I10,5X,3(F11.5)) +* Finish test print +*---- + WZ=1.0D0 + DZ=(TRKORI(3)-BARY(3))/ZMAX +*---- +* Start Z reflection symmetry +* +* IF(ISYMM .GE. 16) THEN +* IF (ABS(DZ) .LT. 1.0E-6) THEN +* WRITE(IOUT,'(A10)') 'ZERO Z ' +* ELSE IF(DZ .LT. 0.0) THEN +* WRITE(IOUT,'(A10)') 'NEGATIVE Z' +* GO TO 170 +* ELSE +* WRITE(IOUT,'(A10)') 'POSITIVE Z' +* WZ=2.0 +* ENDIF +* ENDIF +* Finish Z reflection symmetry +*---- +* +* ELIMINATE TRACKS OUTSIDE CIRCUMSCRIBED CIRCLE + IF( ANN.GT.R2CIRC ) GO TO 170 +* +* WHICH EXTERNAL SURFACES DO THIS TRACK CROSS ? + NTTRK=NTTRK+1 + CALL XELLSR( NDIM, NTOTCL, NSUR, MAXR, REMESH, + > INDEL, MINDIM, MAXDIM, ICOORD, ICUR, INCR, + > TRKORI, TRKDIR, TRKCUT, NSCUT, NCROS, + > TOTLEN) +* +* WHEN NOT SURFACES ARE CROSSED, ELIMINATE THE TRACK + IF(NCROS.LT.2) GO TO 170 + DO 70 K= 1, NDIM + TRKBEG(K)= TRKCUT(K,1) + TRKEND(K)= TRKCUT(K,2) + 70 CONTINUE + DO 75 K= 1, 4 + NSBEG(K)= NSCUT(1) + NSEND(K)= NSCUT(2) + 75 CONTINUE + IF( NCUTOF.NE.1 )THEN + DO 77 K= 1, 4 + DO 76 K3= 1, 3 + TORIC(K3)= TRKORI(K3)+TCUTOF(K3,K) + 76 CONTINUE + CALL XELLSR( NDIM, NTOTCL, NSUR, MAXR, REMESH, + > INDEL, MINDIM, MAXDIM, ICOORD, ICUR, INCR, + > TORIC, TRKDIR, TRKCUT, NSCUT, NCROS, + > TOTXXX) + IF(NSCUT(1).NE.0) NSBEG(K)= NSCUT(1) + IF(NSCUT(2).NE.0) NSEND(K)= NSCUT(2) + 77 CONTINUE + ENDIF + CALL XELLIN( NDIM, NTOTCL, MAXR, REMESH, + > NSUR, NVOL, INDEL, MINDIM, MAXDIM, + > ICOORD, ICUR, INCR, TRKBEG, TRKEND, TRKDIR, + > PROJC2, TOTLEN, + > CONV, LINMAX, LENGHT, NUMERO, LINE) + NOTRAK= NOTRAK+1 +* + WRITE(IFTEMP) 1,LINE+2*NCUTOF,WEIGHT*WZ,IANG0, + > (NSBEG(ISB),ISB=1,NCUTOF), + > (NUMERO(I),I=1,LINE), + > (NSEND(ISE),ISE=1,NCUTOF), + > ( DBLE(1.0/NCUTOF),ISB=1,NCUTOF), + > (LENGHT(I),I=1,LINE), + > ( DBLE(1.0/NCUTOF),ISE=1,NCUTOF) + IF( IPRT.GT.10000)THEN + WRITE(IOUT,6001) NOTRAK, + > NCUTOF,(TRKBEG(JC),JC=1,3), + > NCUTOF,(TRKEND(JC),JC=1,3), + > (TRKDIR(JC),JC=1,3) + WRITE(IOUT,6002) (1.0/NCUTOF,NSBEG(ISB),ISB=1,NCUTOF), + > (LENGHT(I), NUMERO(I),I=1,LINE), + > (1.0/NCUTOF,NSEND(ISE),ISE=1,NCUTOF) + ENDIF + 170 CONTINUE + 180 CONTINUE + 250 CONTINUE + 290 CONTINUE + IF( IPRT.GT.1 )THEN + WRITE(IOUT,'(1H )') + WRITE(IOUT,'(27H0ECHO = TRACKING PROPERTIES )') + WRITE(IOUT,'( 8H0ECHO = ,I3,20H ANGLES AND DENSITY:, + > F9.3,4H/CM2)') + > NANGLE, DENUSR + WRITE(IOUT,'( 8H0ECHO = ,I10,3H / ,I10, + > 23H TRACKS STORED ON TAPE ,I2/)') + > NOTRAK,NTTRK,IFTEMP + ENDIF +* + RETURN + 6001 FORMAT(' #',I8,1P,' B',I1,'(',2(E10.2,','),E10.2,')', + > ' E',I1,'(',2(E10.2,','),E10.2,')', + > ' D(',2(E10.2,','),E10.2,')' ) + 6002 FORMAT(1P,5(1X,E15.7,1X,I6)) + END |