diff options
Diffstat (limited to 'Dragon/src/PSPRAI.f')
| -rw-r--r-- | Dragon/src/PSPRAI.f | 185 |
1 files changed, 185 insertions, 0 deletions
diff --git a/Dragon/src/PSPRAI.f b/Dragon/src/PSPRAI.f new file mode 100644 index 0000000..6dd379e --- /dev/null +++ b/Dragon/src/PSPRAI.f @@ -0,0 +1,185 @@ +*DECK PSPRAI + SUBROUTINE PSPRAI(MXSEG ,NPTS ,XYPOS ,CENTER,RCIRC , + > NSEG ,IORDER,RADANG) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Find general closed Cartesian region intersections +* with annular region and order points for plotting. +* +*Copyright: +* Copyright (C) 2014 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 +* MXSEG maximum number of segments. +* NPTS number of corners. +* XYPOS X and Y position of corners. +* CENTER X and Y position of annulus center. +* RCIRC annulus radius. +* +*Parameters: output +* NSEG number of region intersection. +* number of segments is NSEG-1 +* IORDER type of region: +* = -2 arc segment begins; +* = -1 arc segment ends; +* = 0 close path; +* > 0 corner. +* RADANG segments intersection points +* with respect to annular region center: +* RADANG(1) = radial position; +* RADANG(2) = angular position. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE + INTEGER IOUT + CHARACTER NAMSBR*6 + REAL PI + PARAMETER (IOUT=6,PI=3.1415926535897932,NAMSBR='PSPRAI') +*---- +* ROUTINE PARAMETERS +*---- + INTEGER MXSEG,NPTS,NSEG + INTEGER IORDER(MXSEG) + REAL XYPOS(2,NPTS),CENTER(2),RCIRC, + > RADANG(2,MXSEG) +*---- +* LOCAL PARAMETERS +*---- + INTEGER IPT,ICUR,INXT + REAL XYCUR(2),XYNXT(2),RADCUR,RADNXT, + > XANN,XANNI,XANNF,XANNT,YANNI,YANNF,YANNT, + > DELX,DELY,DELL,RCIRCM +*---- +* SCAN OVER CORNERS AND SIDES +*---- +* write(6,*) ' Circle ',MXSEG ,NPTS ,CENTER(1),CENTER(2),RCIRC + NSEG=1 + IORDER(NSEG)=0 + RCIRCM=0.0 + DO IPT=1,NPTS + ICUR=IPT + INXT=MOD(IPT,NPTS)+1 + XYCUR(1)=XYPOS(1,ICUR)-CENTER(1) + XYCUR(2)=XYPOS(2,ICUR)-CENTER(2) + XYNXT(1)=XYPOS(1,INXT)-CENTER(1) + XYNXT(2)=XYPOS(2,INXT)-CENTER(2) + RADCUR=SQRT(XYCUR(1)*XYCUR(1)+XYCUR(2)*XYCUR(2)) + RADNXT=SQRT(XYNXT(1)*XYNXT(1)+XYNXT(2)*XYNXT(2)) +* write(6,*) ' Point ',IPT ,ICUR ,XYCUR(1),XYCUR(2), +* > XYNXT(1),XYNXT(2),RADCUR,RADNXT,RADCUR-RCIRC + RCIRCM=RCIRC + IF(RADCUR .EQ. RCIRC) RCIRCM=RCIRCM+0.00001 +*---- +* WRITE(6,6002) IPT,XYCUR,XYNXT +* 6002 FORMAT('Line ',I5,5X,'Starts =',2F20.10,5X,'Ends =',2F20.10) +*---- +* CHECK IF CURRENT CORNER IS LOCATED INSIDE +* ANNULAR REGIONS +*---- + IF(RADCUR .LE. RCIRCM) THEN + IF(IORDER(NSEG) .NE. ICUR) THEN +*---- +* IT IS LOCATED INSIDE +* SET IORDER TO IPT TO SPECIFY THIS POINT TO CORRESPOND TO +* CORNER IPT +*---- + NSEG=NSEG+1 + IORDER(NSEG)=ICUR + RADANG(1,NSEG)=RADCUR + IF(RADCUR .EQ. 0.0) THEN + RADANG(2,NSEG)=0.0 + ELSE + RADANG(2,NSEG)=ATAN2(XYCUR(2),XYCUR(1)) + ENDIF + ENDIF + ENDIF +*---- +* Find line direction +*---- + DELY=XYNXT(2)-XYCUR(2) + DELX=XYNXT(1)-XYCUR(1) + DELL=SQRT(DELY*DELY+DELX*DELX) + DELY=DELY/DELL + DELX=DELX/DELL + XANNI=XYCUR(1)*DELX+XYCUR(2)*DELY + YANNI=-XYCUR(1)*DELY+XYCUR(2)*DELX + XANNF=XYNXT(1)*DELX+XYNXT(2)*DELY + YANNF=-XYNXT(1)*DELY+XYNXT(2)*DELX +*---- +* WRITE(6,6003) DELX,DELY,XANNI,YANNI,XANNF,YANNF +* 6003 FORMAT('Rotation ',2F20.10,5X, +* > 'Starts =',2F20.10,5X,'Ends =',2F20.10) +*---- + IF(YANNI .GE. -RCIRCM .AND. YANNI .LE. RCIRCM) THEN + XANN=-SQRT(RCIRCM*RCIRCM-YANNI*YANNI) + IF(XANN .GE. XANNI .AND. + > XANN .LE. XANNF) THEN + NSEG=NSEG+1 + IORDER(NSEG)=-1 + RADANG(1,NSEG)=RCIRCM + XANNT=XANN*DELX-YANNI*DELY + YANNT=XANN*DELY+YANNI*DELX + RADANG(2,NSEG)=ATAN2(YANNT,XANNT) + ENDIF + XANN=-XANN + IF(XANN .GE. XANNI .AND. + > XANN .LE. XANNF) THEN + NSEG=NSEG+1 + IORDER(NSEG)=-2 + RADANG(1,NSEG)=RCIRCM + XANNT=XANN*DELX-YANNI*DELY + YANNT=XANN*DELY+YANNI*DELX + RADANG(2,NSEG)=ATAN2(YANNT,XANNT) + ENDIF + ENDIF +*---- +* CHECK IF NEXT CORNER OF THE RECTANGLE IS LOCATED INSIDE +* ANNULAR REGIONS +*---- + IF(RADNXT .LE. RCIRCM) THEN +*---- +* IT IS LOCATED INSIDE +* SET IORDER TO IPT TO SPECIFY THIS POINT TO CORRESPOND TO +* CORNER IPT +*---- + NSEG=NSEG+1 + IORDER(NSEG) =INXT + RADANG(1,NSEG)=RADNXT + IF(RADNXT .EQ. 0.0) THEN + RADANG(2,NSEG)=0.0 + ELSE + RADANG(2,NSEG)=ATAN2(XYNXT(2),XYNXT(1)) + ENDIF + ENDIF + ENDDO +*---- +* STORE LAST SEGMENT ALSO AT FIRST POSITION +* FOR CYCLIC TRACKING +*---- + IF(NSEG .EQ. 1) THEN + NSEG=2 + IORDER(1)=-2 + RADANG(1,1)=RCIRCM + RADANG(2,1)=0.0 + IORDER(2)=-1 + RADANG(1,2)=RCIRCM + RADANG(2,2)=2.0*PI + ELSE + IORDER(1) =IORDER(NSEG) + RADANG(1,1)=RADANG(1,NSEG) + RADANG(2,1)=RADANG(2,NSEG) + ENDIF +* write(6,'(I5,2F20.10)') (IORDER(IPT),RADANG(1,IPT),RADANG(2,IPT), +* > IPT=1,NSEG) + RETURN + END |