diff options
Diffstat (limited to 'Dragon/src/NXTTPS.f')
| -rw-r--r-- | Dragon/src/NXTTPS.f | 218 |
1 files changed, 218 insertions, 0 deletions
diff --git a/Dragon/src/NXTTPS.f b/Dragon/src/NXTTPS.f new file mode 100644 index 0000000..cc526b6 --- /dev/null +++ b/Dragon/src/NXTTPS.f @@ -0,0 +1,218 @@ +*DECK NXTTPS + SUBROUTINE NXTTPS(IPRINT,NPIN ,IDGPP ,ITSYM ,DRAPIN) +* +*---------- +* +*Purpose: +* To test if pins satisfy required symmetry. +* +*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 +* IPRINT print level. +* NPIN number of pins. +* IDGPP pin direction. +* ITSYM flag for symmetries to test. +* DRAPIN pin position/angle/height/radius. +* +*Reference: +* G. Marleau, +* New Geometries Processing in DRAGON: The NXT: Module, +* Report IGE-260, Polytechnique Montreal, +* Montreal, 2005. +* +*---------- +* + IMPLICIT NONE +*---- +* Subroutine arguments +*---- + INTEGER IPRINT,NPIN,IDGPP,ITSYM(4) + DOUBLE PRECISION DRAPIN(-1:4,NPIN) +*---- +* Functions +*---- + DOUBLE PRECISION XDRCST,PI +*---- +* Local parameters +*---- + INTEGER IOUT + CHARACTER NAMSBR*6 + PARAMETER (IOUT=6,NAMSBR='NXTTPS') + DOUBLE PRECISION DCUTOF,DZERO,DONE,DTWO + PARAMETER (DCUTOF=1.0D-6,DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0) +*---- +* Local variables +*---- + INTEGER IPLOC + INTEGER IS,IP,JP,NPIR + DOUBLE PRECISION DNAP,PIO2,TWOPI + INTEGER ITSYR(4) +*---- +* Processing starts: +* print routine openning output header if required +* and initialize various parameters. +*---- + IPLOC=IPRINT + IF(IPLOC .GE. 100) THEN + WRITE(IOUT,6000) NAMSBR + IF(IPLOC .GE. 1000) THEN + WRITE(IOUT,*) 'Symmetry =',ITSYM + WRITE(IOUT,*) 'Npin =',NPIN + DO IP=1,NPIN + WRITE(IOUT,'(6F20.12)') (DRAPIN(JP,IP),JP=-1,4) + ENDDO + ENDIF + ENDIF +*---- +* Rotate symmetry factors for cell directions +*---- + DO IS=1,4 + ITSYR(IS)=ITSYM(IS) + ENDDO + IF(IDGPP .EQ. 1) THEN + ITSYR(1)=ITSYM(2) + ITSYR(2)=ITSYM(3) + ITSYR(3)=ITSYM(1) + ITSYR(4)=ITSYM(4) + IF(ITSYR(4) .EQ. 1) CALL XABORT(NAMSBR// + > ': X=Y symmetry invalid for pin in direction X') + ELSE IF(IDGPP .EQ. 2) THEN + ITSYR(1)=ITSYM(3) + ITSYR(2)=ITSYM(1) + ITSYR(3)=ITSYM(2) + ITSYR(4)=ITSYM(4) + IF(ITSYR(4) .EQ. 1) CALL XABORT(NAMSBR// + > ': X=Y symmetry invalid for pin in direction Y') + ENDIF + PI=XDRCST('Pi',' ') + PIO2=PI/DTWO + TWOPI=DTWO*PI +*---- +* Scan over symmetrization options (in plane only) +*---- + DO 100 IS=1,3 + IF(ITSYR(IS) .EQ. 1) THEN +*---- +* Scan over symmetrized pin +*---- + IF(IPLOC .GE. 1000) THEN + WRITE(IOUT,'(A3,1X,I8,I8)') 'IS=',IS,ITSYR(IS) + ENDIF + DO 110 IP=1,NPIN +*---- +* Find location of pin after symmetrisation +*---- + IF(IS .EQ. 1) THEN +*---- +* X symmetry +* Symmetric pin should be at \pi-\varphi +*---- + DNAP=PI-DRAPIN(-1,IP) + ELSE IF(IS .EQ. 2) THEN +*---- +* Y symmetry +* Symmetric pin should be at -\varphi +*---- + DNAP=-DRAPIN(-1,IP) + ELSE IF(IS .EQ. 3) THEN +*---- +* Z symmetry +* Symmetric pin should be at \varphi +*---- + DNAP=DRAPIN(-1,IP) + ELSE IF(IS .EQ. 4) THEN +*---- +* X=Y symmetry +* Symmetric pin should be at \pi/2-\varphi +*---- + DNAP=PIO2-DRAPIN(-1,IP) + ENDIF +*---- +* Position angle in range 0 to 2*Pi +*---- + IF(ABS(DNAP) .LE. DCUTOF) THEN + DNAP=DZERO + ELSE IF(DNAP .GT. DCUTOF) THEN + NPIR=INT((DNAP+DCUTOF)/TWOPI) + DNAP=DNAP-DBLE(NPIR)*TWOPI + ELSE + NPIR=INT((DNAP-DCUTOF)/TWOPI) + DNAP=DNAP-DBLE(NPIR-1)*TWOPI + ENDIF + IF(IPLOC .GE. 1000) THEN + WRITE(IOUT,'(A3,1X,I8,2F20.12)') + > 'IP=',IP,DRAPIN(0,IP),DNAP + ENDIF + IF(DRAPIN(0,IP) .LT. DCUTOF) THEN +*---- +* For centered pin, test for radial position only +*---- + DO JP=1,NPIN +*---- +* Verify if pin coincide +* with symmetrized pin +*---- + IF(IPLOC .GE. 1000) THEN + WRITE(IOUT,'(A3,1X,I8,3F20.12)') + > 'JP=',JP,DRAPIN(0,JP),DRAPIN(-1,JP), + > ABS(DRAPIN(0,IP)-DRAPIN(0,JP)) + ENDIF + IF(ABS(DRAPIN(0,IP)-DRAPIN(0,JP)) .LT. DCUTOF) + > GO TO 115 + ENDDO +*---- +* no pin coincide, symmetry not satisfied, abort +*---- + CALL XABORT(NAMSBR//': Symmetric pin not found (C)') + ELSE +*---- +* For pin not centered, test for angular and radial position +*---- + DO JP=1,NPIN +*---- +* Verify if pin coincide +* with symmetrized pin +*---- + IF(IPLOC .GE. 1000) THEN + WRITE(IOUT,'(A3,1X,I8,4F20.8)') + > 'JP=',JP,DRAPIN(0,JP),DRAPIN(-1,JP), + > ABS(DRAPIN(0,IP)-DRAPIN(0,JP)), + > ABS(DNAP-DRAPIN(-1,JP)) + ENDIF + IF(ABS(DRAPIN(0,IP)-DRAPIN(0,JP)) .LT. DCUTOF) THEN + IF(ABS(DNAP-DRAPIN(-1,JP)) .LT. DCUTOF) GO TO 115 + ENDIF + ENDDO +*---- +* no pin coincide, symmetry not satisfied, abort +*---- + CALL XABORT(NAMSBR//': Symmetric pin not found (O-C)') + ENDIF + 115 CONTINUE + 110 CONTINUE + ENDIF + 100 CONTINUE +*---- +* Processing finished: +* print routine closing output header if required +* and return +*---- + IF(IPLOC .GE. 100) THEN + WRITE(IOUT,6001) NAMSBR + ENDIF + RETURN +*---- +* Output formats +*---- + 6000 FORMAT('(* Output from --',A6,'-- follows ') + 6001 FORMAT(' Output from --',A6,'-- completed *)') + END |