path: root/Dragon/src/NXTPRA.f

*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