diff options
Diffstat (limited to 'Dragon/src/NXTPRA.f')
| -rw-r--r-- | Dragon/src/NXTPRA.f | 789 |
1 files changed, 789 insertions, 0 deletions
diff --git a/Dragon/src/NXTPRA.f b/Dragon/src/NXTPRA.f new file mode 100644 index 0000000..092b65a --- /dev/null +++ b/Dragon/src/NXTPRA.f @@ -0,0 +1,789 @@ +*DECK NXTPRA + FUNCTION NXTPRA(NFACES,POSCAR,POSANN,POSPIN,VOLINT) +* +*---------- +* +*Purpose: +* Compute the volume of intersection between +* a 2-D Cartesian region defined by N planes +* an 2-D annular region and +* an annular pin centered at the origin. +* +*Copyright: +* Copyright (C) 2005 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): G. Marleau. +* +*Parameters: input +* NFACES number of planes for Cartesian geometre (3 for triangles, +* 4 for rectangle s and 6 for hexagones). +* POSCAR Cartesian region corner definition: +* POSCAR(1,*) is X position; +* POSCAR(2,*) is Y position; +* POSCAR(*,IPLANE) is location of first corner of +* plane IPLANE; +* POSCAR(*,IPLANE+1) is location of second corner of +* plane IPLANE. +* For last plane, the position of the +* second corner is POSCAR(*,1) +* POSANN spatial description of the annular region with +* POSANN(0) the radius, POSANN(1) the $X$ position +* of center and POSANN(2) the $Y$ position +* of center. +* POSPIN spatial description of the annular pin region with +* POSPIN(0) the radius, POSPIN(1) the $X$ position +* of center and POSPIN(2) the $Y$ position +* of center. +* +*Parameters: output +* NXTPRA type of intersection between the three regions, where: +* = 0 means that the volume of intersection +* between the three regions vanishes; +* =-1 means that the volume of intersection +* between the three regions was computed. +* VOLINT 2-D volume of intersection (area) between the three regions. +* +*Reference: +* G. Marleau, +* New Geometries Processing in DRAGON: The NXT: Module, +* Report IGE-260, Polytechnique Montreal, +* Montreal, 2005. +* +*---------- +* + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + INTEGER NXTPRA + INTEGER NFACES + DOUBLE PRECISION POSCAR(2,NFACES),POSANN(0:2),POSPIN(0:2),VOLINT +*---- +* Functions +*---- + DOUBLE PRECISION XDRCST,PI,PIO2,TPIO2,FPIO2 +*---- +* Local parameters +*---- + INTEGER IOUT,IPRINT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,IPRINT=1,NAMSBR='NXTPRA') + DOUBLE PRECISION DZERO,DONE,DTWO,DHALF,CUTOFF + PARAMETER (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0,DHALF=0.5D0, + > CUTOFF=1.0D-10) +*---- +* Local variables +*---- + DOUBLE PRECISION PX,PY,RADAN2,RACEN2,RACEN,RADPI2, + > COSA,SINA,XIAPR,YIAPR,TRANX,TRANY, + > XYIAPD(2,2) + INTEGER IPOINT,IFACE,IPLIN,JPLIN,NPLIN,ICYL, + > INT,INTT + DOUBLE PRECISION XLOC,YLOC,XMIN,XMAX,COSR,SINR,ACARG + DOUBLE PRECISION POSCYL(2),RADCYL,R2CYL,YBOT,SOL, + > XINT(2),XLOCR,YLOCR + DOUBLE PRECISION XYBEG(2),RADBEG,THBEG,XYEND(2),RADEND,THEND, + > XYADD(2,2),HDT,FACT,DVOL + INTEGER NPOINT,NBPTS,NADSEG,TYADD(2,2),ISEG + DOUBLE PRECISION PNTINT(6),POINTS(2,24) + INTEGER TYPINT(6),TYPES(2,24) + DOUBLE PRECISION X1,Y1,X2,Y2,VLEN,DP1,DP2,DDP + INTEGER NPIN,IP1,IP2 + DOUBLE PRECISION ANGR + +*---- +* Initialize NXTPRA to no intersection +* Initialize PI and multiples +*---- + IF(IPRINT .GE. 200) THEN + WRITE(IOUT,6000) NAMSBR + write(IOUT,6100) + WRITE(IOUT,6101) (POSCAR(1,IP1),POSCAR(2,IP1),IP1=1,NFACES) + WRITE(IOUT,6102) POSCAR(1,1),POSCAR(2,1) + WRITE(IOUT,6103) POSANN(0), + > POSANN(1),POSANN(2) + WRITE(IOUT,6104) POSPIN(0), + > POSPIN(1),POSPIN(2) + ENDIF + COSR=DZERO + SINR=DZERO + XMAX=DZERO + XMIN=DZERO + YLOC=DZERO + NXTPRA=0 + PI=XDRCST('Pi',' ') + PIO2=PI/DTWO + TPIO2=3.0D0*PIO2 + FPIO2=5.0D0*PIO2 +*---- +* Locate annular region/annular pin intersection points +* and find transformation matrix to locate intersection points with +* respect to center of annular region/annular pin intersection. +* In this system of reference the intersection points are located at +* $(0,y_{i})$ and $(0,y_{i})$ or $\theta=\pm \pi/2$ respectively. +* This will become usefull when all the intersection points between +* the rectangular region, the annular region and the annular pins +* must be classified according to a counter-clockwise order. +* See validation in mathematica file: PRA.nb +*---- + RADAN2=POSANN(0)*POSANN(0) + PX=POSANN(1)-POSPIN(1) + PY=POSANN(2)-POSPIN(2) + RACEN2=PX*PX+PY*PY + RACEN=SQRT(RACEN2) + COSA=PX/RACEN + SINA=PY/RACEN + RADPI2=POSPIN(0)*POSPIN(0) + XIAPR=(RADPI2+RACEN2-RADAN2)/(DTWO*RACEN) + YIAPR=SQRT(RADPI2-XIAPR*XIAPR) + TRANX=-XIAPR + TRANY=DZERO +*---- +* Rotate points back to local frame of reference +* XYIAPD(*,1) is at $\theta = \pi/2$ +* XYIAPD(*,2) is at $\theta = -\pi/2$ +*---- + XYIAPD(1,1)=COSA*XIAPR-SINA*YIAPR+POSPIN(1) + XYIAPD(2,1)=SINA*XIAPR+COSA*YIAPR+POSPIN(2) + XYIAPD(1,2)=COSA*XIAPR+SINA*YIAPR+POSPIN(1) + XYIAPD(2,2)=SINA*XIAPR-COSA*YIAPR+POSPIN(2) +*---- +* For each of the faces of the Cartesian region, +* classify the intersection between the annular region, +* annular pins and corners in +* increasing order in a counter-clockwise fashion. +* A maximum of 6 intersection points can be found. +* Note that the first and last corners are identified by $\pm 3$ +* respectively, the first and last intersection point with the +* annular region are indicated by $\pm 1$ respectively +* and the first and last intersection point with the +* annular pin are indicated by $\pm 2$ respectively. +* Procedure to classify the intersection points: +* 1 - Rotate faces in such a way that they are parallel to the +* $X$ axis and below the Cartesian region. +* 2 - Fill in corner locations in increasing order +* 3 - Locate intersection with annular region (after rotation) and +* insert at adequate location in intersection point vector. +* 4 - Locate intersection with annular pin (after rotation) and +* insert at adequate location in intersection point vector. +*---- + IPOINT=0 + DO IFACE=1,NFACES + IP1=IFACE + IP2=MOD(IFACE,NFACES)+1 + ANGR=ATAN2(POSCAR(2,IP2)-POSCAR(2,IP1), + > POSCAR(1,IP2)-POSCAR(1,IP1)) + COSR=COS(-ANGR) + SINR=SIN(-ANGR) +*---- +* Left triangles +*---- +* write(6,'(A20,4F20.15)') +* > 'Angles de rotation ',ANGR,180.0*ANGR/PI,COSR,SINR + XMIN=COSR*POSCAR(1,IP1)-SINR*POSCAR(2,IP1) + XMAX=COSR*POSCAR(1,IP2)-SINR*POSCAR(2,IP2) + YLOC=SINR*POSCAR(1,IP1)+COSR*POSCAR(2,IP1) +* write(6,'(3F20.15)') XMIN,XMAX,YLOC +*---- +* Save corner location +*---- + IPLIN=1 + PNTINT(IPLIN)=XMIN + TYPINT(IPLIN)=4 + IPLIN=2 + PNTINT(IPLIN)=XMAX + TYPINT(IPLIN)=4 + NPLIN=IPLIN +*---- +* Loop over cylinder +* 1- annular region +* 2- annular pin +*---- + DO ICYL=1,2 +*---- +* Extract cylinder information +*---- + IF(ICYL .EQ. 1) THEN +*---- +* Cylinder is annular region +* Rotate as required. +*---- + POSCYL(1)=COSR*POSANN(1)-SINR*POSANN(2) + POSCYL(2)=SINR*POSANN(1)+COSR*POSANN(2) + RADCYL=POSANN(0) + ELSE +*---- +* Cylinder is annular pin +*---- + POSCYL(1)=COSR*POSPIN(1)-SINR*POSPIN(2) + POSCYL(2)=SINR*POSPIN(1)+COSR*POSPIN(2) + RADCYL=POSPIN(0) + ENDIF +*---- +* Find intersection points between Cartesian face and +* cylindrical region +*---- + R2CYL=RADCYL*RADCYL + YBOT=YLOC-POSCYL(2) + SOL=R2CYL-YBOT*YBOT + IF(SOL .GE. DZERO) THEN + XINT(1)=POSCYL(1)-SQRT(SOL) + XINT(2)=POSCYL(1)+SQRT(SOL) +*---- +* Classify intersection points per order of increasing +* x location for annular region and annular pin +*---- + DO INT=1,2 + DO IPLIN=1,NPLIN + IF(XINT(INT) .LE. PNTINT(IPLIN)) THEN + DO JPLIN=NPLIN,IPLIN,-1 + PNTINT(JPLIN+1)=PNTINT(JPLIN) + TYPINT(JPLIN+1)=TYPINT(JPLIN) + ENDDO + PNTINT(IPLIN)=XINT(INT) + TYPINT(IPLIN)=ICYL + GO TO 100 + ENDIF + ENDDO + IPLIN=NPLIN+1 + PNTINT(IPLIN)=XINT(INT) + TYPINT(IPLIN)=ICYL + 100 CONTINUE + NPLIN=NPLIN+1 + ENDDO + ENDIF + ENDDO +*---- +* All intersection points located and ordered for this face +* of the Cartesian region. +* Scan and locate those defining the intersection of three +* regions +* sum of TYPINT = 7 namely: +* +1 -> one annular region +* +2 -> one pin crossing +* +4 -> inside rectangle +* TYPES(1,*) is type of line segment before point +* TYPES(2,*) is type of line segment after point +* Here +* TYPES=1 means annular region, +* TYPES=2 means annular pin and +* TYPES=4 means rectangle side +*---- +* write(6,*) 'AVANT NPLIN',NPLIN +* write(6,'(F20.15,I10)') +* > (PNTINT(IPLIN),TYPINT(IPLIN),IPLIN=1,NPLIN) + INTT=0 + DO IPLIN=1,NPLIN +*---- +* Rotate back line to original location. +*---- + XLOC=PNTINT(IPLIN) + XLOCR=XLOC*COSR+YLOC*SINR + YLOCR=-XLOC*SINR+YLOC*COSR + IF(INTT .EQ. 7) THEN +*---- +* Already in 3 region intersection +* find the point at which one leaves this region +*---- + IPOINT=IPOINT+1 + POINTS(1,IPOINT)=XLOCR + POINTS(2,IPOINT)=YLOCR + TYPES(1,IPOINT)=4 + TYPES(2,IPOINT)=TYPINT(IPLIN) + ENDIF + INTT=INTT+TYPINT(IPLIN) + IF(INTT .EQ. 7) THEN + IPOINT=IPOINT+1 + POINTS(1,IPOINT)=XLOCR + POINTS(2,IPOINT)=YLOCR + TYPES(1,IPOINT)=TYPINT(IPLIN) + TYPES(2,IPOINT)=4 +*---- +* Test if new point is at the same location as previous point +* for rectangle corners and get rid of duplicates +*---- + IF(IPOINT .GE. 2) THEN +* write(6,'(A8,I10)') 'IPOINT ',IPOINT +* write(6,'(A8,2I10,2F20.15)') 'CURRENT ', +* > TYPES(1,IPOINT),TYPES(2,IPOINT), +* > POINTS(1,IPOINT),POINTS(2,IPOINT) +* write(6,'(A8,2I10,2F20.15)') 'PREVIOUS', +* > TYPES(1,IPOINT-1),TYPES(2,IPOINT-1), +* > POINTS(1,IPOINT-1),POINTS(2,IPOINT-1) + IF(TYPES(1,IPOINT) .EQ. 4) THEN + IF(TYPES(1,IPOINT-1) .EQ. 4 .AND. + > TYPES(2,IPOINT-1) .EQ. 4) THEN + DP1=POINTS(1,IPOINT-1)-POINTS(1,IPOINT) + DP2=POINTS(2,IPOINT-1)-POINTS(2,IPOINT) + DDP=SQRT(DP1*DP1+DP2*DP2) +* write(6,*) 'DP1,DP2,DDP',DP1,DP2,DDP + IF(DDP .LT. CUTOFF) THEN +* IF(POINTS(1,IPOINT-1) .EQ. POINTS(1,IPOINT) .AND. +* > POINTS(2,IPOINT-1) .EQ. POINTS(2,IPOINT) ) THEN + IPOINT=IPOINT-1 + ELSE + CALL XABORT(NAMSBR// + > ': Problem with corner position') + ENDIF + ELSE + CALL XABORT(NAMSBR// + > ': Problem with corner order') + ENDIF + ENDIF + ENDIF + ENDIF + ENDDO +* write(6,*) 'APRES NPLIN',NPLIN +* write(6,'(F20.15,I10)') +* > (PNTINT(IPLIN),TYPINT(IPLIN),IPLIN=1,NPLIN) + ENDDO + NPOINT=IPOINT +*---- +* Complete segment for geometry +*---- +* write(6,*) 'NPOINT',NPOINT + IF(NPOINT .LE. 1) THEN +*---- +* Path is empty or contain a single point for intersections with +* sides of the Cartesian region: +* A- Add annular/pin intersection points (XYIAPD) if both inside +* Cartesian region and create path with two arc segments +* 1- From pin to annulus with arc in annular region (1) +* 2- From annulus to pin with arc in pin region (2) +* 3- Closed loop (0) +* B- otherwise, there is no intersection +*---- + NPIN=0 + DO IPOINT=1,2 + DO IFACE=1,NFACES + IP1=IFACE + IP2=MOD(IP1,NFACES)+1 + X1=POSCAR(1,IP2)-POSCAR(1,IP1) + Y1=POSCAR(2,IP2)-POSCAR(2,IP1) + VLEN=SQRT(X1*X1+Y1*Y1) + X2=-Y1/VLEN + Y2=X1/VLEN + VLEN=(XYIAPD(1,IPOINT)-POSCAR(1,IP1))*X2+ + > (XYIAPD(2,IPOINT)-POSCAR(2,IP1))*Y2 + IF(VLEN. LT. DZERO) GO TO 101 + ENDDO + NPIN=NPIN+1 + 101 CONTINUE + ENDDO + IF(NPIN .EQ. 2) THEN +*---- +* TYPES(1,*) is type of line segment before point +* TYPES(2,*) is type of line segment after point +* where TYPES=1 means annular region and +* TYPES=2 means annular pin +*---- + NPOINT=NPOINT+3 + POINTS(1,1)=XYIAPD(1,1) + POINTS(2,1)=XYIAPD(2,1) + TYPES(1,1) =2 + TYPES(2,1) =1 + POINTS(1,2)=XYIAPD(1,2) + POINTS(2,2)=XYIAPD(2,2) + TYPES(1,2) =1 + TYPES(2,2) =2 + POINTS(1,3)=XYIAPD(1,1) + POINTS(2,3)=XYIAPD(2,1) + TYPES(1,2) =2 + TYPES(2,3) =1 + ELSE + NPOINT=0 + ENDIF + ELSE +*---- +* Test for cyclic track if first point is a corner +*---- + IPOINT=1 +* write(6,'(A8,I10)') 'IPOINT ',IPOINT +* write(6,'(A8,2I10,2F20.15)') 'CURRENT ', +* > TYPES(1,IPOINT),TYPES(2,IPOINT), +* > POINTS(1,IPOINT),POINTS(2,IPOINT) +* write(6,'(A8,2I10,2F20.15)') 'LAST ', +* > TYPES(1,NPOINT),TYPES(2,NPOINT), +* > POINTS(1,NPOINT),POINTS(2,NPOINT) + IF(TYPES(1,IPOINT) .EQ. 4 .AND. + > TYPES(2,IPOINT) .EQ. 4) THEN + IF(TYPES(1,NPOINT) .EQ. 4 .AND. + > TYPES(2,NPOINT) .EQ. 4) THEN + DP1=POINTS(1,NPOINT)-POINTS(1,IPOINT) + DP2=POINTS(2,NPOINT)-POINTS(2,IPOINT) + DDP=SQRT(DP1*DP1+DP2*DP2) +* write(6,*) 'DP1,DP2,DDP',DP1,DP2,DDP + IF(DDP .GE. CUTOFF) THEN +* IF(POINTS(1,NPOINT) .NE. POINTS(1,IPOINT) .OR. +* > POINTS(2,NPOINT) .NE. POINTS(2,IPOINT) ) THEN + CALL XABORT(NAMSBR// + > ': Problem with end corner position') + ENDIF + ELSE + CALL XABORT(NAMSBR// + > ': Problem with end corner order') + ENDIF + ELSE +*---- +* Duplicate first point for cyclic track +*---- + NPOINT=NPOINT+1 + POINTS(1,NPOINT)=POINTS(1,IPOINT) + POINTS(2,NPOINT)=POINTS(2,IPOINT) + TYPES(1,NPOINT)=TYPES(1,IPOINT) + TYPES(2,NPOINT)=TYPES(2,IPOINT) + ENDIF + IF(IPRINT .GE. 200) THEN + WRITE(IOUT,6015) + DO IPOINT=1,NPOINT + IF(IPOINT .EQ. NPOINT) THEN + WRITE(IOUT,6011) POINTS(1,IPOINT),POINTS(2,IPOINT), + > TYPES(1,IPOINT),TYPES(2,IPOINT) + ELSE + WRITE(IOUT,6012) POINTS(1,IPOINT),POINTS(2,IPOINT), + > TYPES(1,IPOINT),TYPES(2,IPOINT) + ENDIF + ENDDO + ENDIF +*---- +* Add missing arc segment if required +*---- + NBPTS=NPOINT + DO IPOINT=NPOINT,2,-1 + NADSEG=0 + IF(TYPES(1,IPOINT) .NE. 4) THEN +*---- +* This point finishes an arc segment +* previous point must begin an arc +*---- + IF(TYPES(2,IPOINT-1) .EQ. 4) CALL XABORT(NAMSBR// + > ': Starting point for arc not found') +*---- +* Find position of intersection points with respect to +* annular/pin center location and angular location +* Rotate to center annular region on $X_{+}$ axis (COSA,SINA) +* and translate by (-XIAPR,0) to center +* annular region/annular pin at $x=0$ +*---- + X1=POINTS(1,IPOINT-1)-POSPIN(1) + Y1=POINTS(2,IPOINT-1)-POSPIN(2) +* XYBEG(1)=COSA*POINTS(1,IPOINT-1)+SINA*POINTS(2,IPOINT-1) +* > -XIAPR +* XYBEG(2)=-SINA*POINTS(1,IPOINT-1)+COSA*POINTS(2,IPOINT-1) + XYBEG(1)=COSA*X1+SINA*Y1-XIAPR + XYBEG(2)=-SINA*X1+COSA*Y1 + RADBEG=SQRT(XYBEG(1)*XYBEG(1)+XYBEG(2)*XYBEG(2)) + X2=POINTS(1,IPOINT)-POSPIN(1) + Y2=POINTS(2,IPOINT)-POSPIN(2) +* XYEND(1)=COSA*POINTS(1,IPOINT)+SINA*POINTS(2,IPOINT) +* > -XIAPR +* XYEND(2)=-SINA*POINTS(1,IPOINT)+COSA*POINTS(2,IPOINT) + XYEND(1)=COSA*X2+SINA*Y2-XIAPR + XYEND(2)=-SINA*X2+COSA*Y2 + RADEND=SQRT(XYEND(1)*XYEND(1)+XYEND(2)*XYEND(2)) +*---- +* Find angular location of points +*---- + ACARG=XYBEG(1)/RADBEG + IF(ACARG .GE. 1.0D0) THEN + THBEG=ACOS(1.0D0) + ELSE IF(ACARG .LE. -1.0D0) THEN + THBEG=ACOS(-1.0D0) + ELSE + THBEG=ACOS(ACARG) + ENDIF + IF(XYBEG(2) .LT. DZERO) THBEG=-THBEG + ACARG=XYEND(1)/RADEND + IF(ACARG .GE. 1.0D0) THEN + THEND=ACOS(1.0D0) + ELSE IF(ACARG .LE. -1.0D0) THEN + THEND=ACOS(-1.0D0) + ELSE + THEND=ACOS(ACARG) + ENDIF + IF(XYEND(2) .LT. DZERO) THEND=-THEND + IF(THEND .LT. THBEG) THEND=DTWO*PI+THEND + IF(THBEG .LT. -PIO2) THEN +*---- +* For $\theta_{i}\le -\pi/2$ the segment must be of +* type 1 (annular region) +*---- + IF(TYPES(2,IPOINT-1) .NE. 1) CALL XABORT(NAMSBR// + >': Error -> Initial line segment must be an annular region') +*---- +* For $\theta_{f}\le -\pi/2$ the segment must be of +* type 1 (annular region) and there is no segment +* to add +* For $-\pi/2 < \theta_{f}\le \pi/2$ the segment must be of +* type 2 (annular region) and there is 1 segment +* to add +* For $\pi/2 < \theta_{f}$ the segment must be of +* type 1 (annular region) and there are 2 segments +* to add +*---- + IF(THEND .LT. -PIO2) THEN + IF(TYPES(1,IPOINT) .NE. 1) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be an annular region') + NADSEG=0 + ELSE IF(THEND .LT. PIO2) THEN + IF(TYPES(1,IPOINT) .NE. 2) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be a pin region') + NADSEG=1 + XYADD(1,1)=XYIAPD(1,2) + XYADD(2,1)=XYIAPD(2,2) + TYADD(1,1)=1 + TYADD(2,1)=2 + ELSE + IF(TYPES(1,IPOINT) .NE. 1) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be an annular region') + NADSEG=2 + XYADD(1,1)=XYIAPD(1,2) + XYADD(2,1)=XYIAPD(2,2) + TYADD(1,1)=1 + TYADD(2,1)=2 + XYADD(1,2)=XYIAPD(1,1) + XYADD(2,2)=XYIAPD(2,1) + TYADD(1,2)=2 + TYADD(2,2)=1 + ENDIF + ELSE IF(THBEG .LT. PIO2) THEN +*---- +* For $-\pi/2 < \theta_{i}\le \pi/2$ the segment must be of +* type 2 (pin region) +*---- + IF(TYPES(2,IPOINT-1) .NE. 2) CALL XABORT(NAMSBR// + >': Error -> Initial line segment must be a pin region') +*---- +* For $-\pi/2 < \theta_{f}\le \pi/2$ the segment must be of +* type 2 (pin region) and there is no segment +* to add +* For $\pi/2 < \theta_{f}\le 3\pi/2$ the segment must be of +* type 1 (annular region) and there is 1 segment +* to add +* For $3\pi/2< \theta_{f}$ the segment must be of +* type 2 (pin region) and there are 2 segments +* to add +*---- + IF(THEND .LT. PIO2) THEN + IF(TYPES(1,IPOINT) .NE. 2) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be a pin region') + NADSEG=0 + ELSE IF(THEND .LT. TPIO2) THEN + IF(TYPES(1,IPOINT) .NE. 1) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be an annular region') + NADSEG=1 + XYADD(1,1)=XYIAPD(1,1) + XYADD(2,1)=XYIAPD(2,1) + TYADD(1,1)=2 + TYADD(2,1)=1 + ELSE + IF(TYPES(1,IPOINT) .NE. 2) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be a pin region') + NADSEG=2 + XYADD(1,1)=XYIAPD(1,1) + XYADD(2,1)=XYIAPD(2,1) + TYADD(1,1)=2 + TYADD(2,1)=1 + XYADD(1,2)=XYIAPD(1,2) + XYADD(2,2)=XYIAPD(2,2) + TYADD(1,2)=1 + TYADD(2,2)=2 + ENDIF + ELSE +*---- +* For $\pi/2 < \theta_{i}$ the segment must be of +* type 1 (annular region) +*---- + IF(TYPES(2,IPOINT-1) .NE. 1) CALL XABORT(NAMSBR// + >': Error -> Initial line segment must be an annular region') +*---- +* For $\pi/2 < \theta_{f}\le 3\pi/2$ the segment must be of +* type 1 (annular region) and there is no segment +* to add +* For $3\pi/2< \theta_{f}\le 5*\pi/2$ the segment must be of +* type 2 (pin region) and there is 1 segment +* to add +* For $5*\pi/2 < \theta_{f}$ the segment must be of +* type 1 (annular region) and there are 2 segments +* to add +*---- + IF(THEND .LT. TPIO2) THEN + IF(TYPES(1,IPOINT) .NE. 1) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be an annular region') + NADSEG=0 + ELSE IF(THEND .LT. FPIO2) THEN + IF(TYPES(1,IPOINT) .NE. 2) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be a pin region') + NADSEG=1 + XYADD(1,1)=XYIAPD(1,2) + XYADD(2,1)=XYIAPD(2,2) + TYADD(1,1)=1 + TYADD(2,1)=2 + ELSE + IF(TYPES(1,IPOINT) .NE. 1) CALL XABORT(NAMSBR// + >': Error -> Final line segment must be an annular region') + NADSEG=2 + XYADD(1,1)=XYIAPD(1,2) + XYADD(2,1)=XYIAPD(2,2) + TYADD(1,1)=1 + TYADD(2,1)=2 + XYADD(1,2)=XYIAPD(1,1) + XYADD(2,2)=XYIAPD(2,1) + TYADD(1,2)=2 + TYADD(2,2)=1 + ENDIF + ENDIF + ENDIF +*---- +* Move end segments to create place for new segments +*---- + IF(NADSEG .GT. 0) THEN + DO ISEG=NBPTS,IPOINT,-1 + POINTS(1,ISEG+NADSEG)=POINTS(1,ISEG) + POINTS(2,ISEG+NADSEG)=POINTS(2,ISEG) + TYPES(1,ISEG+NADSEG)=TYPES(1,ISEG) + TYPES(2,ISEG+NADSEG)=TYPES(2,ISEG) + ENDDO +*---- +* Insert new segments +*---- + DO ISEG=NADSEG,1,-1 + POINTS(1,IPOINT+ISEG-1)=XYADD(1,ISEG) + POINTS(2,IPOINT+ISEG-1)=XYADD(2,ISEG) + TYPES(1,IPOINT+ISEG-1)=TYADD(1,ISEG) + TYPES(2,IPOINT+ISEG-1)=TYADD(2,ISEG) + ENDDO + NBPTS=NBPTS+NADSEG + ENDIF + ENDDO + NPOINT=NBPTS + ENDIF + IF(NPOINT .EQ. 0) THEN + NXTPRA=0 + VOLINT=DZERO + ELSE + VOLINT=DZERO + IF(IPRINT .GE. 200) THEN +*---- +* Print cell description if required +*---- + WRITE(IOUT,6010) + DO IPOINT=1,NPOINT + IF(IPOINT .EQ. NPOINT) THEN + WRITE(IOUT,6011) POINTS(1,IPOINT),POINTS(2,IPOINT), + > TYPES(1,IPOINT),TYPES(2,IPOINT) + ELSE + WRITE(IOUT,6012) POINTS(1,IPOINT),POINTS(2,IPOINT), + > TYPES(1,IPOINT),TYPES(2,IPOINT) + ENDIF + ENDDO + ENDIF + DO IPOINT=1,NPOINT-1 +*---- +* Add contribution under line segments +*---- + DVOL=(POINTS(1,IPOINT)-POINTS(1,IPOINT+1)) + > *(POINTS(2,IPOINT)+POINTS(2,IPOINT+1))/DTWO + VOLINT=VOLINT+DVOL + IF(TYPES(2,IPOINT) .EQ. 1) THEN +*---- +* Add annular region contribution (annular region is not centered) +* 1- Find angular width for two points +* 2- Compute volume above line joining the two points +*---- + XYBEG(1)=POINTS(1,IPOINT)-POSANN(1) + XYBEG(2)=POINTS(2,IPOINT)-POSANN(2) + XYEND(1)=POINTS(1,IPOINT+1)-POSANN(1) + XYEND(2)=POINTS(2,IPOINT+1)-POSANN(2) +*---- +* Find angular location of points +*---- + ACARG=XYBEG(1)/POSANN(0) + IF(ACARG .GE. 1.0D0) THEN + THBEG=ACOS(1.0D0) + ELSE IF(ACARG .LE. -1.0D0) THEN + THBEG=ACOS(-1.0D0) + ELSE + THBEG=ACOS(ACARG) + ENDIF + IF(XYBEG(2) .LT. DZERO) THBEG=-THBEG + ACARG=XYEND(1)/POSANN(0) + IF(ACARG .GE. 1.0D0) THEN + THEND=ACOS(1.0D0) + ELSE IF(ACARG .LE. -1.0D0) THEN + THEND=ACOS(-1.0D0) + ELSE + THEND=ACOS(ACARG) + ENDIF + IF(XYEND(2) .LT. DZERO) THEND=-THEND + IF(THEND .LT. THBEG) THEND=DTWO*PI+THEND + HDT=(THEND-THBEG)/DTWO + FACT=COS(HDT)*SIN(HDT) + DVOL=RADAN2*(HDT-FACT) + VOLINT=VOLINT+DVOL + ELSE IF (TYPES(2,IPOINT) .EQ. 2) THEN +*---- +* Add pin region contribution (pin is centered) +* 1- Find angular width for the two points +* 2- Compute volume above line joining the two points +*---- + XYBEG(1)=POINTS(1,IPOINT)-POSPIN(1) + XYBEG(2)=POINTS(2,IPOINT)-POSPIN(2) + XYEND(1)=POINTS(1,IPOINT+1)-POSPIN(1) + XYEND(2)=POINTS(2,IPOINT+1)-POSPIN(2) +*---- +* Find angular location of points +*---- + ACARG=XYBEG(1)/POSPIN(0) + IF(ACARG .GE. 1.0D0) THEN + THBEG=ACOS(1.0D0) + ELSE IF(ACARG .LE. -1.0D0) THEN + THBEG=ACOS(-1.0D0) + ELSE + THBEG=ACOS(ACARG) + ENDIF + IF(XYBEG(2) .LT. DZERO) THBEG=-THBEG + ACARG=XYEND(1)/POSPIN(0) + IF(ACARG .GE. 1.0D0) THEN + THEND=ACOS(1.0D0) + ELSE IF(ACARG .LE. -1.0D0) THEN + THEND=ACOS(-1.0D0) + ELSE + THEND=ACOS(ACARG) + ENDIF + IF(XYEND(2) .LT. DZERO) THEND=-THEND + IF(THEND .LT. THBEG) THEND=DTWO*PI+THEND + HDT=(THEND-THBEG)/DTWO + FACT=COS(HDT)*SIN(HDT) + DVOL=RADPI2*(HDT-FACT) + VOLINT=VOLINT+DVOL + ENDIF + ENDDO + ENDIF + IF(IPRINT .GE. 200) THEN + WRITE(IOUT,6020) VOLINT + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + 6010 FORMAT('FinalSegments={') + 6011 FORMAT('{',F20.10,',',F20.10,',',I10,',',I10,'}};') + 6012 FORMAT('{',F20.10,',',F20.10,',',I10,',',I10,'},') + 6015 FORMAT('OriginalSegments={') + 6020 FORMAT('Volint=',F20.10,';') + 6100 FORMAT('CartesianRegion={') + 6101 FORMAT(('{',F15.10,',',F15.10,'}',','/)) + 6102 FORMAT('{',F15.10,',',F15.10,'}','};') + 6103 FORMAT('RADAN = ',F15.10,';'/ + > 'POSANN={',F15.10,',',F15.10,'};') + 6104 FORMAT('RADIUS= ',F15.10,';'/ + > 'xypin={',F15.10,',',F15.10,'};') + END |