diff options
Diffstat (limited to 'Dragon/src/XCGROD.f')
| -rw-r--r-- | Dragon/src/XCGROD.f | 234 |
1 files changed, 234 insertions, 0 deletions
diff --git a/Dragon/src/XCGROD.f b/Dragon/src/XCGROD.f new file mode 100644 index 0000000..8ee931a --- /dev/null +++ b/Dragon/src/XCGROD.f @@ -0,0 +1,234 @@ +*DECK XCGROD + SUBROUTINE XCGROD(NRT,MSROD,NRODS,RODS,MATROD,RODR) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Check geometry and reorder rod clusters if necessary. +* +*Copyright: +* Copyright (C) 1994 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 +* NRT number of rod types. +* MSROD maximum number of subrods per rods. +* +*Parameters: input/output +* NRODS integer description of rod of a given type: +* NRODS(1,IRT) = number of rod; +* NRODS(2,IRT) = number of subrods in rod; +* NRODS(3,IRT) = first concentric region. +* RODS real description of rod of a given type: +* RODS(1,IRT) = rod center radius; +* RODS(2,IRT) = angular position of first rod. +* MATROD type of material for each subrod. +* RODR subrod radius. +* +*---------------------------------------------------------------------- +* + INTEGER IOUT + REAL PI + PARAMETER (IOUT=6,PI=3.1415926535898) + INTEGER NRT,NRODS(3,NRT),MATROD(MSROD,NRT) + REAL RODS(2,NRT),RODR(MSROD,NRT) + INTEGER, ALLOCATABLE, DIMENSION(:) :: IORD +*---- +* SCRATCH STORAGE ALLOCATION +* IORD : NEW ROD CLUSTER ORDER I(NRT) +*---- + ALLOCATE(IORD(NRT)) +*---- +* CLASSIFY ROD CLUSTER BY INCREASING DISTANCE OF CENTER AND ANGLE +*---- + DO 100 IRT=1,NRT + IORD(IRT)=IRT + 100 CONTINUE + DO 110 IRT=2,NRT + REFR=RODS(1,IRT) + REFA=RODS(2,IRT) + IPOS=IORD(IRT) + DO 111 JRT=IRT-1,1,-1 + KRT=JRT + IF(RODS(1,JRT).GT.REFR) THEN + RODS(1,JRT+1)=RODS(1,JRT) + RODS(2,JRT+1)=RODS(2,JRT) + IORD(JRT+1)=IORD(JRT) + ELSE IF(RODS(1,JRT).EQ.REFR) THEN + IPOS=-IPOS + GO TO 112 + ELSE + GO TO 112 + ENDIF + 111 CONTINUE + KRT=0 + 112 CONTINUE + RODS(1,KRT+1)=REFR + RODS(2,KRT+1)=REFA + IORD(KRT+1)=IPOS + IF(IPOS.LT.0) THEN + DO 113 JRT=KRT,1,-1 + LRT=JRT + IF((RODS(2,JRT).GT.REFA).AND. + > (RODS(1,JRT).EQ.REFR)) THEN + RODS(1,JRT+1)=RODS(1,JRT) + RODS(2,JRT+1)=RODS(2,JRT) + IORD(JRT+1)=IORD(JRT) + ELSE + GO TO 114 + ENDIF + 113 CONTINUE + LRT=0 + 114 CONTINUE + RODS(1,LRT+1)=REFR + RODS(2,LRT+1)=REFA + IORD(LRT+1)=-IPOS + ENDIF + 110 CONTINUE +*---- +* REORDER REMAINING VECTORS NRODS,MATROD,RODR +*---- + DO 140 IRT=1,NRT + JRT=IORD(IRT) + IF(JRT.NE.IRT) THEN + DO 141 IX=1,3 + NNR=NRODS(IX,IRT) + NRODS(IX,IRT)=NRODS(IX,JRT) + NRODS(IX,JRT)=NNR + 141 CONTINUE + DO 142 IS=1,MSROD + MATT=MATROD(IS,IRT) + MATROD(IS,IRT)=MATROD(IS,JRT) + MATROD(IS,JRT)=MATT + RROD=RODR(IS,IRT) + RODR(IS,IRT)=RODR(IS,JRT) + RODR(IS,JRT)=RROD + 142 CONTINUE + DO 143 KRT=IRT+1,NRT + IF(IORD(KRT).EQ.IRT) THEN + IORD(KRT)=JRT + IORD(IRT)=IRT + GO TO 144 + ENDIF + 143 CONTINUE + 144 CONTINUE + ENDIF + 140 CONTINUE +*---- +* FIND IF ROD OVERLAPP +*---- + DO 150 IRT=1,NRT + NRDB=NRODS(1,IRT) + NSBRB=NRODS(2,IRT) + RODRB=RODR(NSBRB,IRT) + RODRB2=RODRB*RODRB + RDPB=RODS(1,IRT) + XBOT=RDPB-RODRB + DANGB=2.*PI/FLOAT(NRDB) + ANGB=RODS(2,IRT) +*---- +* CHECK FOR ROD OVERLAPP INSIDE EACH CLUSTER +*---- + IF(NRDB.GT.1) THEN + IF(RODRB.GT.RDPB) THEN + WRITE(IOUT,'(1X,24HROD OVERLAP IN CLUSTER =,I10)') IRT + CALL XABORT('XCGROD: ROD OVERLAP IN A CLUSTER') + ELSE + ANGMIN=2.*ASIN(RODRB/RDPB) + IF(DANGB.LE.ANGMIN) THEN + WRITE(IOUT,'(1X,24HROD OVERLAP IN CLUSTER =,I10)') IRT + CALL XABORT('XCGROD: ROD OVERLAP IN A CLUSTER') + ENDIF + ENDIF + ENDIF +*---- +* CHECK FOR ROD OVERLAPP BETWEEN DIFFERENT CLUSTERS +*---- + DO 151 JRT=IRT-1,1,-1 + NRDT=NRODS(1,JRT) + NSBRT=NRODS(2,JRT) + RODRT=RODR(NSBRT,JRT) + RODRT2=RODRT*RODRT + RDPT=RODS(1,JRT) + XTOP=RDPT+RODRT + DANGT=2.*PI/FLOAT(NRDT) + ANGT=RODS(2,JRT) +*---- +* NO OVERLAPP +*---- + IF(XTOP.LT.XBOT) GO TO 152 +*---- +* SOME OVERLAPP POSSIBLE TEST FOR INTERSECTION +*---- + ANG1=ANGB + DO 160 IA1=1,NRDB +*---- +* FIND POSITION OF ROD (X0,Y0) +*---- + X01=RDPB*COS(ANG1) + Y01=RDPB*SIN(ANG1) + RRX=RODRB2-X01*X01 + RRY=RODRB2-Y01*Y01 + XY=X01*Y01 + RR1=(RRX-Y01*Y01) + ANG2=ANGT + DO 161 IA2=1,NRDT + X02=RDPT*COS(ANG2) + Y02=RDPT*SIN(ANG2) + RR2=(RODRT2-X02*X02-Y02*Y02) +*---- +* CHECK FOR ROD INSIDE ROD +*---- + DELX=X02-X01 + DELY=Y02-Y01 + DIST=SQRT(DELX**2+DELY**2) + IF(DIST.LT.RODRT+RODRB) THEN + WRITE(IOUT,'(1X,25HROD OVERLAP IN CLUSTERS =,2I10)') + > IRT,JRT + CALL XABORT('XCGROD: ROD OVERLAP IN 2 CLUSTERS') + ENDIF +*---- +* FIND IF CIRCLES +* (X-X01)**2+(Y-Y01)**2=RODRB*2 +* (X-X02)**2+(Y-Y02)**2=RODRT*2 +* INTERSECT +*---- + IF(X02.NE.X01) THEN + CCR=1./DELX + BBR=-DELY*CCR + AAR=0.5*CCR*(RR1-RR2) + ARGSQ=AAR*(2.*X01-2.*BBR*Y01-AAR) + > +BBR*(BBR*RRY+2.*XY)+RRX + ELSE + CCR=1./DELY + BBR=-DELX*CCR + AAR=0.5*CCR*(RR1-RR2) + ARGSQ=AAR*(2.*Y01-2.*BBR*X01-AAR) + > +BBR*(BBR*RRX+2.*XY)+RRY + ENDIF + IF(ARGSQ.GE.0.0) THEN + WRITE(IOUT,'(1X,25HROD OVERLAP IN CLUSTERS =,2I10)') + > IRT,JRT + CALL XABORT('XCGROD: ROD OVERLAP IN 2 CLUSTERS') + ENDIF + ANG2=ANG2+DANGT + 161 CONTINUE + ANG1=ANG1+DANGB + 160 CONTINUE + 151 CONTINUE + 152 CONTINUE + 150 CONTINUE +*---- +* SCRATCH STORAGE DEALLOCATION +*---- + DEALLOCATE(IORD) +*---- +* RETURN +*---- + RETURN + END |