! !--------------------------------------------------------------------- ! !Purpose: ! Support module to compute a single track. ! !Copyright: ! Copyright (C) 2001 Ecole Polytechnique de Montreal. ! !Author(s): ! X. Warin ! !--------------------------------------------------------------------- ! MODULE SAL_TRAJECTORY_MOD USE PRECISION_AND_KINDS, ONLY : PDB,TWOPI,SMALL,PI USE SAL_TRACKING_TYPES, ONLY : NBER,LGOK,COSINE,AX,AY,EX,EY,ALPHA,F0,AT,BT, & LGTYPE,R,D0,TORIG USE SAL_NUMERIC_MOD, ONLY : SALACO CONTAINS SUBROUTINE SALTRA(DANGLT,NPERIM_MAC2,PERIM_MAC2,ISURF2_ELEM,IPAR,RPAR,PPERIM, & PERIM,IBC2_ELEM,IDATA_BC2,BCDATA,PPERIM_MAC2,DIST_AXIS) ! !--------------------------------------------------------------------- ! !Purpose: ! computes intersection of trajectory (T): R=A+D*E with a mesh composed ! of nodes and elements ! ! begin ! | ! first entrance (LGMORE=f), get entry point ! | ! get boundary condition at entry point and treat it ! | ! if there is a left domain,do left tracking,inverse the trajectory ! | ! if LGGEO1=t,test if there is a left re-entry point,keep it ! | ! do basic tracking ! | ! if LGGEO1=t,test if there is a right re-entry point,keep it ! | ! end ! !Parameters: input ! DANGLT angle cosines ! NPERIM_MAC2 number of elements composing perimeter of domain ! PERIM_MAC2 elements composing perimeter of domain ! ISURF2_ELEM relative 2D surf number per elem ! IPAR integer geometry descriptors ! RPAR REAL(PDB) geometry descriptors ! PPERIM array pointer to elements in the perimeter of nodes ! PERIM array of elements in perimeter of nodes ! !Parameters: input (optional data for cyclic tracking) ! IBC2_ELEM relative 2D boundary condition indices per element ! IDATA_BC2 position of data per 2D boundary condition ! BCDATA table of boundary condition descriptor ! PPERIM_MAC2 pointer to 'perim' and 'dist_axis': ! PPERIM_MAC2(1): beginning of elements on axis 1 ! PPERIM_MAC2(2): beginning of elements on axis 2 ! PPERIM_MAC2(3): beginning of elements not on axis ! DIST_AXIS distance of points on this axis to the center (0,0) ! !--------------------------------------------------------------------- ! USE SAL_GEOMETRY_TYPES, ONLY : ISPEC,TYPGEO USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2,ITRAC2,ANGTAB,ELMTAB,CNT,CNT0,NB_TOT,DNEW,DINIT, & NNEW,LNEW,IERR,LGMORE,DD0,NTRACK,EPS1,EX0,EY0,LGOK,IPART,DELX, & N_AXIS IMPLICIT NONE REAL(PDB), INTENT(IN), DIMENSION(:,:) :: DANGLT INTEGER, INTENT(IN) :: NPERIM_MAC2 INTEGER, INTENT(IN), DIMENSION(:) :: PERIM_MAC2,PPERIM,PERIM,ISURF2_ELEM INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR INTEGER, INTENT(IN), DIMENSION(:), OPTIONAL :: IBC2_ELEM,IDATA_BC2,PPERIM_MAC2 REAL(PDB), INTENT(IN), DIMENSION(:), OPTIONAL :: DIST_AXIS REAL(PDB), INTENT(IN), DIMENSION(:,:), OPTIONAL :: BCDATA !*** INTEGER :: IOUT,LEN,P1,P2,IPHI LOGICAL :: LGLEFT REAL(PDB) :: RADIA !*** ! initiate nber of sub-trajectories NB_TOT=0 IF(ISPEC == 0) THEN !* compute entry point for a basic trajectory CALL SAL240_3(PERIM_MAC2,NPERIM_MAC2,IPAR,RPAR) ELSE EX=EX0; EY=EY0 LGOK=.FALSE. ! initiate all elements status to untreated IPART(1,:)=-1 ! radia of the axis IF(PPERIM_MAC2(N_AXIS+1)==PPERIM_MAC2(N_AXIS)) THEN RADIA=0. ELSE RADIA=DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1) ENDIF IF(DELX=NMAX2) THEN CALL XABORT('SALTRA: Buffer overflow') ELSE ITRAC2(CNT+1:CNT+LEN)=ANGTAB(1:LEN) ITRAC2(CNT+NMAX2+1:CNT+NMAX2+LEN)=ELMTAB(1:LEN) CNT=CNT+LEN ENDIF ! END SUBROUTINE SALTRA SUBROUTINE SAL241_2(NPERIM,PERIM,DIST_AXIS,IPAR) ! !--------------------------------------------------------------------- ! !Purpose: ! compute an entry point on the axial element ! !Parameters: input ! NPERIM = number of elements on this axis ! PERIM = elements on this axis in perimeter ! DIST_AXIS = distance of points on this axis to the center (0,0) ! IPAR = integer descriptor of elements ! !--------------------------------------------------------------------- ! USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO USE SAL_TRACKING_TYPES, ONLY : IPART,N_AXIS,DNEW,DELX,NNEW,LNEW,COSINE,AX, & AY,HX,HY,BX,BY,EX,EY INTEGER, INTENT(IN) :: NPERIM INTEGER, INTENT(IN), DIMENSION(:) :: PERIM REAL(PDB), INTENT(IN), DIMENSION(:) :: DIST_AXIS INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR INTEGER :: I,J CHARACTER(LEN=131) :: HSMG !*** LNEW=0 !* compute crossed element DO I=1,NPERIM IF(DELX<=DIST_AXIS(I)) THEN LNEW=PERIM(I) EXIT ENDIF ENDDO IF(LNEW==0) THEN WRITE(HSMG,'(52HSAL241_2: Error of distances on the axis for typgeo=, & & i3,1h.)') TYPGEO CALL XABORT(HSMG) ENDIF !* get entered node NNEW=IPAR(2,LNEW) IF(NNEW<0) NNEW=IPAR(3,LNEW) IF(NNEW<0) CALL XABORT('SAL241_2: Error of element data') !* compute DNEW at entry point DNEW=DELX*(EX*HX(N_AXIS)+EY*HY(N_AXIS)) IF(N_AXIS>2) DNEW=DNEW+BX(N_AXIS)*EX+BY(N_AXIS)*EY !* compute COSINE COSINE=ABS(HX(N_AXIS)*EY-HY(N_AXIS)*EX) !* set all elements in this axis to be 'treated (0)' ! others are 'untreated' IPART(1,:)=-1 DO I=1,NPERIM J=PERIM(I) IPART(1,J)=0 ENDDO !* initial point AX=BX(N_AXIS)+DELX*HX(N_AXIS)-DNEW*EX AY=BY(N_AXIS)+DELX*HY(N_AXIS)-DNEW*EY ! END SUBROUTINE SAL241_2 ! SUBROUTINE SAL240_3(PERIM_MAC2,NPERIM_MAC2,IPAR,RPAR) ! !--------------------------------------------------------------------- ! !Purpose: ! compute entry point for a basic trajectory ! !Parameters: input ! PERIM_MAC2 elements composing perimeter of domain ! NPERIM_MAC2 number of elements composing perimeter of domain ! IPAR integer geometry descriptors ! RPAR floating point geometry descriptors ! !--------------------------------------------------------------------- ! USE SAL_TRACKING_TYPES, ONLY : IPART,EX0,EY0,EX,EY,DD0,LGOK,LGMORE,NBER, & DNEW,NNEW,LNEW,DINIT IMPLICIT NONE INTEGER, INTENT(IN) :: NPERIM_MAC2 INTEGER, INTENT(IN), DIMENSION(:) :: PERIM_MAC2 REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR !*** EX=EX0; EY=EY0 ! enter trajectory with vectors a and e and initial distance dd0 ! initiate all elements status to untreated IPART(1,:)=-1 CALL SAL241(PERIM_MAC2,NPERIM_MAC2,IPAR,RPAR,DD0,-100,DNEW,NNEW,LNEW) DINIT=DNEW ! if we have sevaral entry points LGMORE=NBER>3 !* if not succeed IF(.NOT.LGOK) CALL XABORT('SAL240_3: Could not enter domain') ! END SUBROUTINE SAL240_3 ! SUBROUTINE SAL245(ISURF2_ELEM,IPAR,RPAR,IOUT,LOLD,NOLD,DOLD,LGLEFT) ! !--------------------------------------------------------------------- ! !Purpose: ! boundary condition treatment at entry point ! !Parameters: input ! ISURF2_ELEM relative 2D surf nber per elem ! IPAR integer geometry descriptors ! RPAR floating point geometry descriptors ! !Parameters: input/output ! LOLD entry point ! NOLD node index ! DOLD distance ! LGLEFT flag set to TRUE if there is a left domain ! !Parameters: output ! IOUT (to load outside surface) = 5 (left trajectory) ! 6 (right trajectory) ! !--------------------------------------------------------------------- ! USE SAL_TRACKING_TYPES, ONLY : ITRAC2,RTRAC2,CNT0,COSINE,EX,EY USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE,ISPEC IMPLICIT NONE INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR INTEGER, INTENT(IN), DIMENSION(:) :: ISURF2_ELEM REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR INTEGER, INTENT(INOUT) :: LOLD,NOLD REAL(PDB), INTENT(INOUT) :: DOLD LOGICAL, INTENT(INOUT) :: LGLEFT INTEGER, INTENT(OUT) :: IOUT !*** INTEGER :: BCIN,SURF !*** ! initiate surface info ITRAC2(CNT0+5)=0 ITRAC2(CNT0+6)=0 !* get boundary condition at entry surface BCIN=IPAR(2,LOLD) IF(NOLD==BCIN)BCIN=IPAR(3,LOLD) ! add entry cosine anyhow (for characteristics) RTRAC2(CNT0+2)=COSINE LGLEFT=.FALSE. IF(ISPEC == 1) RETURN IF(BCIN>=G_BC_TYPE(0).OR.BCIN==G_BC_TYPE(-1).OR.BCIN==G_BC_TYPE(1))THEN !* trajectory enters through vacuum: ! - compute angle, trajectory-normal and store cos and sin ! (inverse vector E in order to pass outgoing direction) CALL SAL247_1(RPAR(:,LOLD),IPAR(:,LOLD),DOLD,RTRAC2(CNT0+3),RTRAC2(CNT0+5),-EX,-EY) !* get index of the entry surface SURF=ISURF2_ELEM(LOLD) IF(SURF/=0) ITRAC2(CNT0+5)=SURF ! store 2D surface index IOUT=6 ELSE CALL XABORT('SAL245: Reversed direction') ENDIF ! END SUBROUTINE SAL245 ! SUBROUTINE SAL240_4(PERIM,PPERIM,IPAR,RPAR,IPHI,ISURF2_ELEM,IOUT) ! !--------------------------------------------------------------------- ! !Purpose: ! track a trajectory until leaving the domain ! begin: having (nnew,dnew,lnew) ! | ! (1)write horizontal angle information ! | ! (2)compute successive crossed nodes ! | ! (3)boundary condition: ! re-entry: get re-entry point,go to (1) ! go out: compute surface number,end ! !Parameters: input ! PERIM array of elements in perimeter of nodes ! PPERIM array pointer to elements in the perimeter of nodes ! IPAR integer geometry descriptors ! RPAR floating point geometry descriptors ! IPHI angular index of the track ! ISURF2_ELEM relative 2D surf number per elem ! IOUT (to load outside surface) = 5 (left trajectory) ! !--------------------------------------------------------------------- ! USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2,ITRAC2,RTRAC2,ANGTAB,ELMTAB,PRTIND,CNT, & CNT0,EX,EY,LGOK,NB_TOT,NB_MAX,DNEW,NNEW,LNEW,IERR USE SAL_GEOMETRY_TYPES, ONLY : G_BC_TYPE IMPLICIT NONE ! IN VARIABLE INTEGER, INTENT(IN) :: IPHI INTEGER, INTENT(IN), DIMENSION(:) :: PPERIM,PERIM INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR INTEGER, INTENT(IN), DIMENSION(:) :: ISURF2_ELEM INTEGER, INTENT(IN) :: IOUT !*** INTEGER :: NOLD,LOLD,P1,P2,CNT1,SURF REAL(PDB) :: DOLD,LENGTH LOGICAL :: LGON INTEGER, PARAMETER :: FOUT =6 !*** LGON = .TRUE. ! initiate counter CNT1=CNT EXTERIOR : DO WHILE(LGON) NB_TOT=NB_TOT+1 IF(NB_TOT > NB_MAX) CALL XABORT('SAL240_4: NB_TOT overflow') ! keep horizontal phi nber in angtab IF(NB_TOT>1) CALL XABORT('SAL240_4: Angtab overflow') ANGTAB(2*NB_TOT)=IPHI ELMTAB(2*NB_TOT-1)=LNEW ; ELMTAB(2*NB_TOT)=0 ; ! !* track a sub-trajectory INTERIOR: DO WHILE(NNEW>0) ! UPDATE DATA TO COMPUTE NEXT NODE: DOLD=DNEW LOLD=LNEW NOLD=NNEW ! crossing NODE NOLD ! input: trajectory (T):R=A+D*E => A = (AX,AY), E = (EX,EY) ! DOLD = D at last intersection ! NOLD = NODE just entered P1=PPERIM(NOLD); P2=PPERIM(NOLD+1)-1 CALL SAL241(PERIM(P1:P2),P2-P1+1,IPAR,RPAR,DOLD,NOLD,DNEW,NNEW,LNEW) ! at return from SAL241: ! DNEW = D at point exiting node ! COSINE = cosine of trajectory with exiting normal ! NNEW = new node entered ! LNEW = element crossed when exiting node ! NBER = number of intersections with perimeter ! LGOK = .TRUE. if trajectory exits the node IF(.NOT.LGOK) THEN IERR=0 IF(PRTIND>0)WRITE(FOUT,'("SAL240_4 ==> couldnt exit node ",I5)') NOLD RETURN ENDIF ! store data LENGTH=DNEW-DOLD ! store new length in track arrays CNT=CNT+1 IF(CNT>=NMAX2) CALL XABORT('SAL240_4: NMAX2 overflow') RTRAC2(CNT)=LENGTH ITRAC2(CNT+NMAX2)=LNEW ITRAC2(CNT)=NOLD END DO INTERIOR ! !* exiting motif and analyzing bc condition LGON =(NNEW=G_BC_TYPE(5)) ! STORE NBER OF REGIONS TO ANGTAB ANGTAB(2*NB_TOT-1)=CNT-CNT1 CNT1=CNT IF(LGON) THEN CALL XABORT('SAL240_4: Lgon is true') ELSE ! compute leaving surface IF(NNEW<=0)THEN ! we got a surface: end of trajectory ! vacuum bd condition. put a marker (surface number) to allow ! psi and pss computation ! store exiting surface, cosphi and sinphi ! compute angle trajectory-normal CALL SAL247_1(RPAR(:,LNEW),IPAR(:,LNEW),DNEW, & RTRAC2(CNT0+IOUT-2),RTRAC2(CNT0+IOUT),EX,EY) ! store 2D surface-cone nber SURF=ISURF2_ELEM(LNEW) IF(SURF/=0) ITRAC2(CNT0+IOUT)=SURF ENDIF ENDIF ! ENDDO EXTERIOR ! set success flag IERR=1 ! END SUBROUTINE SAL240_4 ! SUBROUTINE SAL240_4_2(DANGLT,PPERIM,PERIM,IPAR,RPAR,IDATA_BC2,BCDATA, & PERIM_MAC2,PPERIM_MAC2,DIST_AXIS,IBC2_ELEM) ! !--------------------------------------------------------------------- ! !Purpose: ! track a cyclic trajectory until the track length equal to the predefined ! length of this trajectory ! begin: having (nnew,dnew,lnew) ! | ! (1)write horizontal angle information ! | ! (2)compute successive crossed nodes ! | ! (3)if abs(total_length*length_inv_cycl-1)>eps1 : ! yes:continue, get re-entry point,go to (1) ! no: end of tracking ! !Parameters: input ! DANGLT angle cosines ! PERIM array of elements in perimeter of nodes ! PPERIM array pointer to elements in the perimeter of nodes ! IPAR integer geometry descriptors ! RPAR real geometry descriptors ! IDATA_BC2 position of bc data per 2D boundary conditions ! PERIM_MAC2 elements composing perimeter of domain ! PPERIM_MAC2 pointer to 'PERIM' and 'DIST_AXIS': ! PPERIM_MAC2(1):beginning of elements on axis 1 ! PPERIM_MAC2(2):beginning of elements on axis 2 ! PPERIM_MAC2(3):beginning of elements not on axis ! DIST_AXIS distance of points on this axis to the center (0,0) ! BCDATA table of bc descriptor ! IBC2_ELEM relative 2D bc nber per elem ! !Parameters: output ! NB_TOT total number of sub-trajectories ! ANGTAB table of {N_K,ANGLE_K} (1 NB_MAX) CALL XABORT('SAL240_4_2: NB_TOT overflow') N_AXIS_KEEP(NB_TOT)=N_AXIS ! set horizontal angle index in angtab ANGTAB(2*NB_TOT)=ANGLE_TO_NUMBER(EX,EY,DANGLT) ELMTAB(2*NB_TOT-1)=LNEW ; ELMTAB(2*NB_TOT)=0 ; TORIG(1,NB_TOT)=AX+DNEW*EX ; TORIG(2,NB_TOT)=AY+DNEW*EY ; IF(IMPX > 4) WRITE(6,*) 'SAL240_4_2: beginning of track=',TORIG(:2,NB_TOT) ! !* track a sub-trajectory INTERIOR: DO WHILE(NNEW.GT.0) ! update data to compute next node: DOLD=DNEW LOLD=LNEW NOLD=NNEW ! crossing node NOLD ! input: trajectory (t):r=a+d*e => a = (ax,ay), e = (ex,ey) ! DOLD = d at last intersection ! NOLD = node just entered P1=PPERIM(NOLD); P2=PPERIM(NOLD+1)-1 CALL SAL241(PERIM(P1:P2),P2-P1+1,IPAR,RPAR,DOLD,NOLD,DNEW,NNEW,LNEW) ! at return from SAL241: ! DNEW = d at point exiting node ! COSINE = cosine of trajectory with exiting normal ! NNEW = new node entered ! LNEW = element crossed when exiting node ! NBER = nber of intersections with perimeter ! LGOK = .true. if trajectory exits the node IF(.NOT.LGOK) THEN IERR=0 WRITE(FOUT,'(" SAL240_4_2 ==> couldnt exit node ",I5)') NOLD RETURN ENDIF ! store data LENGTH=DNEW-DOLD ! add to total length LENGTH_TOT=LENGTH_TOT+LENGTH CNT=CNT+1 IF(CNT>=NMAX2) THEN ALLOCATE(ITRAC3(4*NMAX2),RTRAC3(2*NMAX2),STAT=OK) IF(OK/=0) CALL XABORT('SAL240_4_2: NMAX2 overflow.') RTRAC3(:NMAX2)=RTRAC2(:NMAX2) ITRAC3(:2*NMAX2)=ITRAC2(:2*NMAX2) DEALLOCATE(RTRAC2,ITRAC2) RTRAC2=>RTRAC3 ITRAC2=>ITRAC3 NMAX2=2*NMAX2 ENDIF RTRAC2(CNT)=LENGTH ITRAC2(CNT+NMAX2)=LNEW ITRAC2(CNT)=NOLD ENDDO INTERIOR ! !* exiting motif and analyzing bc condition LGON=NNEW<=G_BC_TYPE(1).AND.NNEW>=G_BC_TYPE(5).AND.(ABS(LENGTH_TOT*LENGTH_INV_CYCL-1.)>EPS1) ! store nber of regions to angtab ANGTAB(2*NB_TOT-1)=CNT-CNT1 CNT1=CNT IF(LGON)THEN ! treat boundary condition and get new entry point IF(PRESENT(IDATA_BC2).AND.PRESENT(DIST_AXIS)) THEN IDATA=IDATA_BC2(IBC2_ELEM(LNEW)) ! treat bondary condition CALL SAL247_3(BCDATA(:,IDATA)) ! at return from SAL247_3: ! AX AY = entry point at boundary ! EX EY = new direction ! compute entry point IF(NNEW==G_BC_TYPE(3).OR.NNEW==G_BC_TYPE(4).OR.NNEW==G_BC_TYPE(2)) THEN P1=PPERIM_MAC2(N_AXIS); P2=PPERIM_MAC2(N_AXIS+1)-1 ! re-entry point is on the axis CALL SAL241_2(P2-P1+1,PERIM_MAC2(P1:P2),DIST_AXIS(P1:P2),IPAR) ENDIF ELSE CALL XABORT('SAL240_4_2: missing IDATA_BC2 or DIST_AXIS argument') ENDIF ENDIF ! ENDDO EXTERIOR ! set success flag IERR=1 ! END SUBROUTINE SAL240_4_2 ! SUBROUTINE SAL247_3(BCDATA) ! !---------------------------------------------------------------------- ! !Purpose: ! Treatment of boundary conditions for the cyclic cases ! side 4 ! side 3 ------ ! ---------- /\ side 3 / \ side 5 ! | | / \ / \ ! side 2 | | side 4 side 2 / \ side 3 \ / ! | | / \ side 2 \ / side 6 ! ---------- -------- ------ ! side 1 side 1 side 1 ! !Parameters: input ! BCDATA boundary condition descriptor ! !---------------------------------------------------------------------- ! USE SAL_TRACKING_TYPES, ONLY : NNEW,DNEW,AX,AY,BX,BY,EX,EY,DELX,N_AXIS USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,EPS,LX=>LENGTHX IMPLICIT NONE REAL(PDB), INTENT(IN), DIMENSION(:) :: BCDATA !*** REAL(PDB) :: AUX,ACX,ACY,COSTHE,SINTHE,CX,CY,EX0 !*** ! boundary point: AX=AX+DNEW*EX AY=AY+DNEW*EY CX=BCDATA(1) ; CY=BCDATA(2) COSTHE=BCDATA(3) SINTHE=BCDATA(4) SELECT CASE(NNEW) CASE(-2,-3) ! translation/rotation for a segment (get displacement data) ! a = displacement vector ! new axis SELECT CASE(TYPGEO) CASE(5) IF (CY<0) THEN ! end axe 2 N_AXIS=1 DELX=AX ELSEIF (CY>0) THEN ! end axe 1 N_AXIS=3 DELX=AX ELSEIF (CX<0) THEN ! end axe 4 N_AXIS=2 DELX=AY ELSEIF (CX>0) THEN ! end axe 3 N_AXIS=4 DELX=AY ENDIF CASE(9) IF (ABS(BCDATA(1))0) THEN IF (CY>0) THEN N_AXIS=5 ELSE N_AXIS=6 ENDIF ELSEIF (CX<0) THEN IF (CY>0) THEN N_AXIS=3 ELSE N_AXIS=2 ENDIF ENDIF ACX=AX+CX ; ACY=AY+CY DELX=SQRT((ACX-BX(N_AXIS))**2+(ACY-BY(N_AXIS))**2) CASE(10) ACX=LX-AX EX0=EX IF(CX>0) THEN ! cx>0: axis 3 EX=0.5*EX+0.5*SQRT(3.)*EY EY=0.5*EY-0.5*SQRT(3.)*EX0 N_AXIS=2 DELX=SQRT((AX-CX)*(AX-CX)+AY*AY) ELSEIF(BCDATA(5)>0.) THEN ! cx=0, angle>0: axis 2 (axis of angle pi/3) EX=0.5*EX-0.5*SQRT(3.)*EY EY=0.5*EY+0.5*SQRT(3.)*EX0 N_AXIS=3 DELX=SQRT(AX*AX+AY*AY) ELSE ! cx=0, angle=0: axis 1 (axis X) EX=-EX EY=-EY N_AXIS=1 DELX=ACX ENDIF CASE(11) ACX=AX-CX ; ACY=AY-CY AUX=2._PDB*(ACX*COSTHE+ACY*SINTHE) AX=AUX*COSTHE-ACX+CX AY=AUX*SINTHE-ACY+CY DELX=LX-SQRT((AX-CX)*(AX-CX)+(AY-CY)*(AY-CY)) EX0=EX IF(ABS(BCDATA(5))0.) THEN N_AXIS=2 ! cy>0: exit element is on axis 3 (angle=0) ELSE N_AXIS=4 ! cy=0: exit element is on axis 1 (angle=0) ENDIF ELSE ! exit side is on axes 2 or 4 EX=-0.5*EX-0.5*SQRT(3.)*EY EY=-0.5*EY+0.5*SQRT(3.)*EX0 AX=AX-1.5*(LX-DELX) AY=AY-0.5*SQRT(3.)*(LX-DELX) IF(CX>0.) THEN N_AXIS=1 ! cx>0: exit element is on axis 4 (angle>0) ELSE N_AXIS=3 ! cx=0: exit element is on axis 2 (angle>0) ENDIF ENDIF CASE DEFAULT CALL XABORT('SAL247_3: option not available(1)') END SELECT CASE(-4) ! symmetry with respect an axis (get axis data) ! axis: r=c+t*f (f unit vector of angle theta) ACX=AX-CX ; ACY=AY-CY AUX=2._PDB*(ACX*COSTHE+ACY*SINTHE) AX=AUX*COSTHE-ACX+CX AY=AUX*SINTHE-ACY+CY AUX=2._PDB*(EX*COSTHE+EY*SINTHE) EX=AUX*COSTHE-EX EY=AUX*SINTHE-EY SELECT CASE(TYPGEO) CASE(6) IF(BCDATA(5)>0.) THEN ! vertical axes IF(CX>0.) THEN N_AXIS=4 ELSE N_AXIS=2 ENDIF DELX=AY ELSE ! horizontal axes IF(CY>0.) THEN N_AXIS=3 ELSE N_AXIS=1 ENDIF DELX=AX ENDIF CASE(7) IF(CX>0) THEN ! cx>0: axis 3 (vertical axis) N_AXIS=3 DELX=AY ELSEIF(BCDATA(5)>0.) THEN ! cx=0, angle>0: axis 2 (axis of angle pi/4) N_AXIS=2 DELX=SQRT(AX*AX+AY*AY) ELSE ! cx=0, angle=0: axis 1 (axis x) N_AXIS=1 DELX=SQRT(AX*AX+AY*AY) ENDIF CASE(8,12) IF(CX>0) THEN ! cx>0: axis 3 N_AXIS=3 DELX=SQRT((AX-CX)*(AX-CX)+AY*AY) ELSEIF(BCDATA(5)>0.) THEN ! cx=0, angle>0: axis 2 (axis of angle pi/3 or 2*pi/3) N_AXIS=2 DELX=SQRT(AX*AX+AY*AY) ELSE ! cx=0, angle=0: axis 1 (axis X) N_AXIS=1 DELX=AX ENDIF CASE DEFAULT CALL XABORT('SAL247_3: option not available(2)') END SELECT CASE DEFAULT CALL XABORT('SAL247_3: option not available(3)') END SELECT ! END SUBROUTINE SAL247_3 ! INTEGER FUNCTION ANGLE_TO_NUMBER(EX,EY,DANGLT) ! !--------------------------------------------------------------------- ! !Purpose: ! search order number of a horizontal angle in the angular quadrature ! formula set ! !Parameters: input ! EX angle cosine ! EY angle sine ! DANGLT angle cosines table ! !--------------------------------------------------------------------- ! IMPLICIT NONE REAL(PDB), INTENT(IN) :: EX,EY REAL(PDB), INTENT(IN), DIMENSION(:,:) :: DANGLT !*** INTEGER :: I,NPHI REAL(PDB) :: EXREF,EYREF CHARACTER(LEN=131) :: HSMG !*** NPHI=SIZE(DANGLT,2) ANGLE_TO_NUMBER=0 DO I=1,NPHI EXREF=DANGLT(1,I) ; EYREF=DANGLT(2,I) ; IF((ABS(EX-EXREF)<1.E-3).AND.(ABS(EY-EYREF)<1.E-3)) THEN ANGLE_TO_NUMBER=I GO TO 10 ELSE IF((ABS(EX-EXREF)<1.E-3).AND.(ABS(EY+EYREF)<1.E-3)) THEN ANGLE_TO_NUMBER=2*NPHI-I+1 GO TO 10 ENDIF ENDDO WRITE(6,'(/29H ANGLE_TO_NUMBER: QUADRATURE:)') DO I=1,2*NPHI WRITE(6,'(1X,I5,1P,2E12.4)') I,DANGLT(:2,I) ENDDO WRITE(HSMG,'(47HANGLE_TO_NUMBER: UNABLE TO FIND ANGULAR COSINES,1P,2E12.4, & & 26H INTO QUADRATURE SELECTED.)') EX,EY CALL XABORT(HSMG) 10 RETURN ! END FUNCTION ANGLE_TO_NUMBER ! SUBROUTINE SAL247_1(RPAR,IPAR,D,SINPHI,COSPHI,EX,EY) ! !--------------------------------------------------------------------- ! !Purpose: ! computes COSPHI and SINPHI for the intersection of the trajectory ! with element LNEW of descriptors RPAR and IPAR. SINPHI is ! computed with respect to directions outgoing with the trajectory ! !Parameters: input ! RPAR floating point geometry descriptors ! IPAR integer geometry descriptors ! D distance measured along the trajectory ! EX first exiting unit vector in along trajectory ! EY second exiting unit vector in along trajectory ! !Parameters: output ! SINPHI cosine at intersection ! COSPHI sine at intersection ! !--------------------------------------------------------------------- ! USE SAL_TRACKING_TYPES, ONLY : AX,AY IMPLICIT NONE INTEGER, INTENT(IN), DIMENSION(:) :: IPAR REAL(PDB), INTENT(IN), DIMENSION(:) :: RPAR REAL(PDB), INTENT(IN) :: D,EX,EY REAL(PDB), INTENT(OUT) :: SINPHI,COSPHI !> AX AY = components of origin of trajectory !*** INTEGER :: TYPE REAL(PDB) :: NX,NY !*** NX=0._PDB NY=0._PDB TYPE=IPAR(1) SELECT CASE(TYPE) CASE(1) ! TYPE=1=> segment (s): R=C+T*F with T in (0,1) ! RPAR(1),RPAR(2) = C = (CX,CY) ! RPAR(3),RPAR(4) = F = (FX,FY) ! components of normal NX=RPAR(4)/RPAR(5) NY=-RPAR(3)/RPAR(5) CASE(2,3) ! TYPE=2,3=> arc of circle NX=(AX+D*EX-RPAR(1))/RPAR(3) NY=(AY+D*EY-RPAR(2))/RPAR(3) END SELECT SINPHI=NX*EY-NY*EX IF((EX*NX+EY*NY)<0._PDB) SINPHI=-SINPHI COSPHI=SQRT(1._PDB-SINPHI*SINPHI) ! END SUBROUTINE SAL247_1 ! SUBROUTINE SAL241(PERIM,NPERIM,IPAR,RPAR,DOLD,NOLD,DNEW,NNEW,LNEW) ! !--------------------------------------------------------------------- ! !Purpose: ! computes intersection of trajectory (t): R=A+D*E with a perimeter ! composed of the elements given in array perim ! !Parameters: input ! PERIM elements in the perimeter of a node ! NPERIM number of elements in the perimeter ! RPAR floating point geometry descriptors ! IPAR integer geometry descriptors ! DOLD value of D at current position ! NOLD current entered node ! !Parameters: output ! DNEW value of D at the intersection ! NNEW node index that is enter after intersection ! LNEW element index that is intersected ! !--------------------------------------------------------------------- ! USE SAL_GEOMETRY_TYPES, ONLY : NRPAR,NIPAR USE SAL_TRACKING_TYPES, ONLY : NRPART,NIPART,RPART,IPART,EPS1 !*** IMPLICIT NONE INTEGER, INTENT(IN), DIMENSION(:) :: PERIM INTEGER, INTENT(IN), DIMENSION(:,:) :: IPAR REAL(PDB), INTENT(IN), DIMENSION(:,:) :: RPAR INTEGER, INTENT(IN) :: NPERIM,NOLD INTEGER, INTENT(OUT) :: NNEW,LNEW REAL(PDB), INTENT(IN) :: DOLD REAL(PDB), INTENT(OUT) :: DNEW !*** REAL(PDB) :: D INTEGER :: I,L,INDEX,N,TYPE,NBINTE ! INFTY is used to initialize search for minimum distance REAL(PDB), PARAMETER :: INFTY=1.E+10 INTEGER, PARAMETER :: FOUT =6 !**** ! initialize distance for intersection DNEW=INFTY NBER=0 ! intersection: DO I=1,NPERIM L=PERIM(I) ! get order nber of element in the perimeter NBINTE=IPART(1,L) IF(NBINTE<0)THEN ! compute and store intersections TYPE=IPAR(1,L) IF(TYPE==1)THEN ! segment CALL SAL242(RPAR(:,L),IPAR(:,L),RPART(:,L),IPART(2:,L),NBINTE) ELSEIF(TYPE<=3)THEN ! arc of circle or circle CALL SAL243(RPAR(:,L),IPAR(:,L),RPART(:,L),IPART(2:,L),NBINTE) ELSE CALL XABORT('SAL241: Not implemented') ENDIF IPART(1,L)=NBINTE ENDIF ! IF(NBINTE/=0)THEN DO INDEX=1,NBINTE D=RPART(INDEX,L) N=IPART(1+INDEX,L) ! analyzes feasability of intersection to eliminate conflicts due ! to concavities and crossing of mesh points IF(N/=NOLD.AND.D>(DOLD-EPS1))THEN NBER=NBER+1 IF(D<(DNEW-EPS1))THEN DNEW=D COSINE=RPART(INDEX+2,L) NNEW=N LNEW=L ELSEIF(D<(DNEW+EPS1))THEN ! case of two close intersections : liquidate smaller NBER=NBER-1 IF(D>DNEW)THEN DNEW=D COSINE=RPART(INDEX+2,L) NNEW=N LNEW=L ENDIF ENDIF ENDIF ENDDO ENDIF ENDDO ! LGOK=DNEW/=INFTY IF(LGOK)THEN ! eliminate intersection IPART(1,LNEW)=IPART(1,LNEW)-1 IF(IPART(1,LNEW)==1.AND.RPART(1,LNEW)==DNEW)THEN ! FOR AN ELEMENT WITH TWO INTERSECTIONS, MOVE 2ND ! INTERSECTION INTO FIRST IF FIRST HAS BEEN TAKEN RPART(1,LNEW)=RPART(2,LNEW) RPART(3,LNEW)=RPART(4,LNEW) IPART(2,LNEW)=IPART(3,LNEW) ENDIF ELSE ! print out problem IF(NOLD/=-100) THEN WRITE(FOUT,*)'Problem in SAL241: NOLD, DOLD, NBINTE ',NOLD,DOLD,NBINTE ENDIF ENDIF END SUBROUTINE SAL241 ! SUBROUTINE SAL242(RPAR,IPAR,D,NODE,NBINTE) ! !--------------------------------------------------------------------- ! !Purpose: ! analysis of the intersection of the trajectory (T):R=A+D*E and ! a segment ! !Parameters: input ! RPAR floating point geometry descriptors ! IPAR integer geometry descriptors ! !Parameters: output ! D value of D at intersection ! NODE order nber of node entered after intersection ! NBINTE number of intersections ! !--------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN), DIMENSION(:) :: IPAR INTEGER, INTENT(OUT), DIMENSION(:) :: NODE INTEGER, INTENT(OUT) :: NBINTE REAL(PDB), INTENT(IN), DIMENSION(:) :: RPAR REAL(PDB), INTENT(OUT), DIMENSION(:) :: D ! DIMENSION RPAR(*),IPAR(*),NODE(*),D(*) !*** REAL(PDB) :: CAX,CAY,FX,FY,A,DELTAM,DELTA REAL(PDB), PARAMETER :: EPS2=0. !*** ! TYPE=1=> segment (S): R=C+T*F with T in (0,1) ! RPAR(1),RPAR(2) = C = (CX,CY) ! RPAR(3),RPAR(4) = F = (FX,FY) ! components of vector F FX=RPAR(3) FY=RPAR(4) ! DELTA=F X E DELTA=FX*EY-FY*EX IF(DELTA/=0._PDB)THEN ! components of vector CA=C-A CAX=RPAR(1)-AX CAY=RPAR(2)-AY ! A=AC X E A=CAY*EX-CAX*EY DELTAM=DELTA-A IF(DELTAM>(-EPS2).AND.A>(-EPS2))THEN ! crossing into + halfspace NODE(1)=IPAR(3) ELSEIF(DELTAM<=EPS2.AND.A<=EPS2)THEN ! crossing into - halfspace NODE(1)=IPAR(2) ELSE ! out-of-range crossing NBINTE=0 RETURN ENDIF ! compute distance to intersection along trajectory ! D = (AC X F)/DELTA D(1)=(CAY*FX-CAX*FY)/DELTA D(3)=ABS(DELTA) NBINTE=1 ELSE ! A/=0 => out-of-range crossing (infinity) ! A==0 => trajectory coincides with segment ! in any case neglect intersection NBINTE=0 ENDIF END SUBROUTINE SAL242 ! SUBROUTINE SAL243(RPAR,IPAR,D,NODE,NBINTE) ! !--------------------------------------------------------------------- ! !Purpose: ! analysis of the intersection of the trajectory (T):R=A+D*E and ! a circle or arc of circle ! !Parameters: input ! RPAR floating point geometry descriptors ! IPAR integer geometry descriptors ! !Parameters: output ! D value of D at intersection ! NODE order nber of node entered after intersection ! NBINTE number of intersections ! !--------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN), DIMENSION(:) :: IPAR INTEGER, INTENT(OUT), DIMENSION(:) :: NODE INTEGER, INTENT(OUT) :: NBINTE REAL(PDB), INTENT(IN), DIMENSION(:) :: RPAR REAL(PDB), INTENT(OUT), DIMENSION(:) :: D ! DIMENSION RPAR(*),IPAR(*),NODE(*),D(*) !*** INTEGER :: TYPE,I REAL(PDB) :: CAX,CAY,RAD,RAD2,RHOMI2,DMIN,THETA,THETA1,THETA2, & COSTHE,DELTA REAL(PDB), PARAMETER :: EPS2=0._PDB !*** ! TYPE=2,3=> arc of circle (C): R=C+R*F(THETA), ! THETA in (THETA1,THETA2) ! RPAR(1),RPAR(2) = C = (CX,CY) ! RPAR(3) = R = RADIUS ! RPAR(4),RPAR(5) = (THETA1,THETA2) in (0,2PI) with THETA1 THETA1 > 0 ! = .FALSE. => THETA1 < 0 ! components of vector CA=C-A CAX=RPAR(1)-AX CAY=RPAR(2)-AY TYPE=IPAR(1) ! value of R2 RAD=RPAR(3) RAD2=RAD**2 ! RHOMI2=(CA X E)**2 RHOMI2=(CAX*EY-CAY*EX)**2 IF(RAD2>=RHOMI2)THEN DELTA=SQRT(RAD2-RHOMI2) ! tangent point = two very close points (to avoid infinite loop) if(delta==0._pdb) delta=small ! DMIN=CA*E DMIN=CAX*EX+CAY*EY ! D(1) D(2) = min and max distances for the two intersections D(1)=DMIN-DELTA D(2)=DMIN+DELTA D(3)=DELTA/RAD D(4)=D(3) IF(TYPE==2)THEN ! full circle. both intersections are possible NODE(1)=IPAR(2) NODE(2)=IPAR(3) NBINTE=2 ELSE ! analysis for arc of circle: THETA1=RPAR(4) THETA2=RPAR(5) NBINTE=0 ! compute angles for closest and farthest intersections ! and check feasability DO I=1,2 COSTHE=(D(I)*EX-CAX)/RAD THETA=SALACO(COSTHE,D(I)*EY-CAY) IF( (((THETA1-THETA)