summaryrefslogtreecommitdiff
path: root/Dragon/src/TLM.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/TLM.f')
-rw-r--r--Dragon/src/TLM.f425
1 files changed, 425 insertions, 0 deletions
diff --git a/Dragon/src/TLM.f b/Dragon/src/TLM.f
new file mode 100644
index 0000000..01f5863
--- /dev/null
+++ b/Dragon/src/TLM.f
@@ -0,0 +1,425 @@
+*DECK TLM
+ SUBROUTINE TLM(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Create Matlab procedure to trace the integration lines
+* generated with the NXT tracking module of DRAGON.
+*
+*Copyright:
+* Copyright (C) 2006 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
+* NENTRY number of data structures transfered to this module.
+* HENTRY name of the data structures.
+* IENTRY data structure type where:
+* =1 for LCM memory object;
+* =2 for XSM file;
+* =3 for sequential binary file;
+* =4 for sequential ASCII file.
+* JENTRY access permission for the data structure where:
+* =0 for a data structure in creation mode;
+* =1 for a data structure in modifications mode;
+* =2 for a data structure in read-only mode.
+* KENTRY data structure pointer.
+*
+*Comments:
+* Instructions for the use of the TLM: module:
+* M-file.m := TLM: [ M-file.m ] VOLTRK TRKFIL ::
+* [ EDIT [ iprint ] ]
+* [ NTPO nplots ]
+* (TLMget) ;
+* where
+* M-file.m : SEQ_ASCII file containing Matlab instructions.
+* VOLTRK : read-only tracking data structure
+* (signature L_TRACK).
+* TRKFIL : read-only sequential binary tracking file.
+* EDIT : keyword to specify print level.
+* iprint : print level. By default, iprint=1.
+* NTPO : keyword to specify number of plots generated
+* by this execution.
+* nplots : number of plots. By default, nplots=1.
+* (TLMget) : Processing options to select types of plots.
+* (read from input using the TLMGET routine).
+*
+*Reference:
+* G. Marleau,
+* New Geometries Processing in DRAGON: The NXT: and TLM: Modules,
+* Report IGE-260, Polytechnique Montreal,
+* Montreal, 2005.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+ IMPLICIT NONE
+*----
+* Subroutine arguments
+*----
+ INTEGER NENTRY,IENTRY(NENTRY),JENTRY(NENTRY)
+ TYPE(C_PTR) KENTRY(NENTRY)
+ CHARACTER HENTRY(NENTRY)*12
+*----
+* Local parameters
+*----
+ INTEGER IOUT
+ CHARACTER NAMSBR*6
+ PARAMETER (IOUT=6,NAMSBR='TLM ')
+ INTEGER ILCMUP,ILCMDN
+ PARAMETER (ILCMUP=1,ILCMDN=2)
+ INTEGER NSTATE,NIPLP
+ PARAMETER (NSTATE=40,NIPLP=6)
+ DOUBLE PRECISION DZERO,DONE
+ PARAMETER (DZERO=0.0D0,DONE=1.0D0)
+*----
+* Variables for input via REDGET
+*----
+ INTEGER ITYPLU,INTLIR
+ CHARACTER CARLIR*72,CARLST*72
+ REAL REALIR
+ DOUBLE PRECISION DBLLIR
+*----
+* Allocatable arrays
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: MATALB,NUMERG
+ INTEGER, ALLOCATABLE, DIMENSION(:,:) :: IPLP
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DPLP,DGMESH,DANGLT,
+ 1 VOMER1,VOMERG
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) ::DVNOR1,DVNOR
+*----
+* Local functions
+*----
+ INTEGER TLMVPL
+ INTEGER IVALID
+*----
+* Local variables
+*----
+ TYPE(C_PTR) IPTRK
+ INTEGER IMTRK,IFTRK,IMFTRK,IMMAT,IPMAT
+ INTEGER IEN
+ CHARACTER HSIGN*12,CMAT*4
+ INTEGER ISTATT(NSTATE),IEDIMG(NSTATE)
+ CHARACTER TITLE*72
+ INTEGER IPRINT,NPLOTS
+ CHARACTER CTYPE*4,COMNT*80
+ INTEGER NCOMNT,NBTR,ICOM
+ INTEGER NDIM,ISPEC,NREG,NSOUT,NALBG,NCOR,NANGL,NFSUR,
+ > NFREG,IFMT,MXSUB,MXSEG
+ INTEGER NBUCEL,NUCELL(3),MAXMSH,MAXMSP,MXGREG,MAXMDH
+ INTEGER II,KK,ITRKT,IRENOT,NBDR,NSUR,IPLOT
+ INTEGER NSKTRK
+ INTEGER ITGEO,ILONG,ITYLCM
+ DOUBLE PRECISION XYZL(2,3)
+ LOGICAL LMIX
+*----
+* Validate entry parameters
+*----
+ IF(NENTRY .NE. 3) CALL XABORT(NAMSBR//
+ > ': Three data structures required')
+ IPTRK=C_NULL_PTR
+ IMTRK=0
+ IFTRK=0
+ IMFTRK=0
+ IEN=1
+ IMMAT=JENTRY(IEN)
+ IPMAT=FILUNIT(KENTRY(IEN))
+ IF(IENTRY(IEN) .NE. 4 ) CALL XABORT(NAMSBR//
+ > ': Matlab .m file is not an ASCII file')
+ IF(IMMAT .NE. 0 .AND.
+ > IMMAT .NE. 1 ) CALL XABORT(NAMSBR//
+ > ': Matlab .m file not in create or update mode')
+*----
+* Scan data structure to determine type and mode
+*----
+ DO IEN=2,NENTRY
+ IF(IENTRY(IEN) .EQ. 1 .OR. IENTRY(IEN) .EQ. 2) THEN
+ IF(JENTRY(IEN) .EQ. 2) THEN
+ CALL LCMGTC(KENTRY(IEN),'SIGNATURE',12,HSIGN)
+ IF(HSIGN .EQ. 'L_TRACK') THEN
+ IPTRK=KENTRY(IEN)
+ IMTRK=IEN
+ CALL LCMGTC(KENTRY(IEN),'TRACK-TYPE',12,HSIGN)
+ IF((HSIGN .NE. 'EXCELL').AND.(HSIGN .NE. 'MCCG')) THEN
+ CALL XABORT(NAMSBR//
+ > ': Tracking data structure type is invalid')
+ ENDIF
+ ELSE
+ CALL XABORT(NAMSBR//
+ > ': Invalid signature for '//HENTRY(IEN))
+ ENDIF
+ ELSE
+ CALL XABORT(NAMSBR//
+ > ': Tracking data structure not in read-only mode')
+ ENDIF
+ ELSE IF(IENTRY(IEN) .EQ. 3) THEN
+ IF(JENTRY(IEN) .NE. 2) CALL XABORT(NAMSBR//
+ > ': Tracking file not in read-only mode')
+ IFTRK=FILUNIT(KENTRY(IEN))
+ IMFTRK=IEN
+ ELSE
+ CALL XABORT(NAMSBR//
+ > ': Invalid data structure format for '//HENTRY(IEN))
+ ENDIF
+ ENDDO
+ IF(IMFTRK .EQ. 0) CALL XABORT(NAMSBR//
+ >': No tracking file available')
+ IF(IMTRK .EQ. 0) CALL XABORT(NAMSBR//
+ >': No Tracking data structure available')
+*----
+* Recover EDIT level and number of track processing option
+* [ EDIT [ iprint ] ]
+* [ NTPO [ nplots ] ]
+* by default iprint=1 and nplots=1.
+*----
+ IPRINT=1
+ NPLOTS=1
+ LMIX=.FALSE.
+ 1010 CONTINUE
+ CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLLIR)
+ 1011 CONTINUE
+ IF(ITYPLU .EQ. 10) GO TO 1015
+ IF(ITYPLU .NE. 3) CALL XABORT(NAMSBR//
+ > ': Read error -- Character variable expected')
+ IF(CARLIR(1:4) .EQ. ';') THEN
+ GO TO 1015
+ ELSE IF(CARLIR(1:4) .EQ. 'EDIT') THEN
+ CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLLIR)
+ IF(ITYPLU .NE. 1) GO TO 1011
+ IPRINT=INTLIR
+ ELSE IF(CARLIR(1:7) .EQ. 'MIXTURE') THEN
+ LMIX=.TRUE.
+ ELSE IF(CARLIR(1:4) .EQ. 'NTPO') THEN
+ CALL REDGET(ITYPLU,INTLIR,REALIR,CARLIR,DBLLIR)
+ IF(ITYPLU .NE. 1) GO TO 1011
+ NPLOTS=INTLIR
+ ELSE
+ GO TO 1015
+ ENDIF
+ GO TO 1010
+ 1015 CONTINUE
+*----
+* Get Matlab plot options
+*----
+ CARLST=CARLIR
+ ALLOCATE(IPLP(NIPLP,NPLOTS),DPLP(4*NPLOTS))
+ CALL TLMGET(IPRINT,NPLOTS,NDIM,CARLST,IPLP,DPLP)
+*----
+* Read tracking file parameters
+*----
+ READ(IFTRK) CTYPE,NCOMNT,NBTR,IFMT
+ IF(CTYPE .NE. '$TRK') CALL XABORT(NAMSBR//
+ >': Binary file is not a valid NXT: tracking file')
+ IF(IFMT .NE. 1) CALL XABORT(NAMSBR//': IFMT.NE.1')
+ ITRKT=1
+ IRENOT=1
+ DO ICOM=1,NCOMNT
+ READ(IFTRK) COMNT
+ IF(COMNT(1:12) .EQ. 'TRKNOR ' ) THEN
+ IF(COMNT(15:26) .EQ. 'Directional ' ) THEN
+ IRENOT=-1
+ ELSE IF(COMNT(15:26) .EQ. 'Global ' ) THEN
+ IRENOT=0
+ ENDIF
+ ELSE IF(COMNT(1:12) .EQ. 'OPTION ' ) THEN
+ IF(COMNT(15:26) .EQ. 'Extended ' ) THEN
+ ITRKT=0
+ ENDIF
+ ENDIF
+ ENDDO
+ IF(ITRKT .NE .0) CALL XABORT(NAMSBR//
+ >': Insufficient information on tracking file'//
+ >' Use EDIT -1000 in NXT:')
+ READ(IFTRK) NDIM,ISPEC,NREG,NSOUT,NALBG,NCOR,NANGL,MXSUB,MXSEG
+ READ(IFTRK) KK
+ ALLOCATE(MATALB(-NSOUT:NREG))
+ READ(IFTRK) (MATALB(II),II=-NSOUT,NREG)
+ DO II=1,4
+ READ(IFTRK) KK
+ ENDDO
+ NSKTRK=NCOMNT+8
+ REWIND(IFTRK)
+*----
+* Initialize tracking parameters to 0
+*----
+ ISTATT(:NSTATE)=0
+ IEDIMG(:NSTATE)=0
+*----
+* Read state vectors
+*----
+ CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATT)
+ CALL LCMGTC(IPTRK,'TITLE',72,TITLE)
+ IF(ISTATT( 7) .NE. 4) CALL XABORT(NAMSBR//
+ >': Tracking data structure incompatible with current module')
+ NSUR=ISTATT(5)
+ IRENOT=ISTATT(8)
+ MXSEG=MAX(MXSEG,ISTATT(18))
+ CALL LCMSIX(IPTRK,'NXTRecords ',ILCMUP)
+ CALL LCMGET(IPTRK,'G00000001DIM',IEDIMG)
+ IF(ISTATT( 1) .NE. NREG .OR.
+ > ISTATT( 5) .NE. NSOUT .OR.
+ > ISTATT(21) .NE. NANGL .OR.
+ > IEDIMG( 1) .NE. NDIM ) THEN
+ WRITE(IOUT,9000) ISTATT( 1),NREG ,ISTATT( 5),NSOUT,
+ > ISTATT(21),NANGL,
+ > IEDIMG( 1),NDIM
+ CALL XABORT(NAMSBR//
+ >': Tracking data structure and file do mot match')
+ ENDIF
+ ITGEO=ABS(IEDIMG( 2))
+ IF(IPRINT .GE. 1) THEN
+ WRITE(IOUT,6000) NAMSBR,NREG,NSOUT,NANGL,NDIM,ITGEO
+ ENDIF
+ NBUCEL=IEDIMG( 5)
+ NUCELL(1)=IEDIMG(13)
+ NUCELL(2)=IEDIMG(14)
+ NUCELL(3)=IEDIMG(15)
+ MAXMSH =IEDIMG(16)
+ MAXMSP =IEDIMG(20)
+ NFSUR =IEDIMG(22)
+ NFREG =IEDIMG(23)
+ MXGREG =IEDIMG(25)
+ MAXMDH=MAX(MAXMSH,MAXMSP,MXGREG)
+ ALLOCATE(DGMESH((MAXMDH+2)*4))
+ CALL TLMGEO(IPTRK,IPMAT,IPRINT,ITGEO,MAXMDH,NDIM,NUCELL,
+ > DGMESH,XYZL)
+ ALLOCATE(DANGLT(NDIM*NANGL))
+ CALL LCMGET(IPTRK,'TrackingDirc',DANGLT)
+ NBDR=1
+ IF(IRENOT .EQ. -1) NBDR=NBDR+NANGL
+ ALLOCATE(DVNOR(NREG,NBDR),DVNOR1(NFREG,NBDR))
+ DVNOR1(:NFREG,:NBDR)=DONE
+ IF(IRENOT .EQ. -1) THEN
+ CALL LCMGET(IPTRK,'VTNormalizeD',DVNOR1(1,2))
+ ELSE IF(IRENOT .EQ. 0) THEN
+ CALL LCMGET(IPTRK,'VTNormalize ',DVNOR1(1,1))
+ ENDIF
+*----
+* Merge normalization factors
+*----
+ CALL LCMLEN(IPTRK,'NumMerge ',ILONG,ITYLCM)
+ IF(ILONG.EQ.NFREG) THEN
+ ALLOCATE(VOMER1(NFREG),NUMERG(NFREG),VOMERG(NREG))
+ CALL LCMGET(IPTRK,'VolMerge ',VOMER1)
+ CALL LCMGET(IPTRK,'NumMerge ',NUMERG)
+ VOMERG(:NREG)=DZERO
+ DVNOR(:NREG,:NBDR)=DZERO
+ DO II=1,NFREG
+ KK=NUMERG(II)
+ IF(KK.EQ.0) CYCLE
+ IF(KK.GT.NREG) CALL XABORT(NAMSBR//': DVNOR overflow')
+ VOMERG(KK)=VOMERG(KK)+VOMER1(II)
+ DVNOR(KK,:)=DVNOR(KK,:)+DVNOR1(II,:)*VOMER1(II)
+ ENDDO
+ DO II=1,NREG
+ DVNOR(II,:)=DVNOR(II,:)/VOMERG(II)
+ ENDDO
+ DEALLOCATE(VOMERG,NUMERG,VOMER1)
+ ELSE
+ IF(NREG.NE.NFREG) CALL XABORT(NAMSBR//': INVALID VALUE OF NREG')
+ DVNOR(:,:)=DVNOR1(:,:)
+ ENDIF
+ DEALLOCATE(DVNOR1)
+ CALL LCMSIX(IPTRK,'NXTRecords ',ILCMDN)
+*----
+* Read IPMAT to end-of-file and
+* insert pause if in update mode
+*----
+ IF(IMMAT .EQ. 1) THEN
+ 1000 CONTINUE
+ READ(IPMAT,'(A4)',END=1005) CMAT
+ GO TO 1000
+ 1005 CONTINUE
+ WRITE(IPMAT,7010)
+ ENDIF
+*----
+* Write execution comments on IPMAT
+*----
+ WRITE(IPMAT,7000) NAMSBR,HENTRY(IMTRK),HENTRY(IMFTRK),TITLE
+*----
+* Loop over PLOTS
+*----
+ DO IPLOT=1,NPLOTS
+ IF(ABS(IPLP(1,IPLOT)) .EQ. 1) THEN
+*----
+* POINTS
+*----
+ CALL TLMPNT(IPMAT ,IFTRK ,IPRINT,NSKTRK,NBTR ,NDIM ,
+ > NREG ,NSUR ,MXSUB ,MXSEG ,NANGL ,NBDR ,
+ > NPLOTS,IPLOT ,IPLP ,DANGLT,DVNOR )
+ ELSE IF(ABS(IPLP(1,IPLOT)) .EQ. 2) THEN
+*----
+* DIRECTIONS
+*----
+ CALL TLMDIR(IPMAT ,IFTRK ,IPRINT,ISPEC, NSKTRK,NBTR ,
+ > NDIM ,NSOUT, NREG ,MXSUB ,MXSEG ,NANGL ,
+ > NBDR ,NPLOTS,IPLOT ,IPLP ,DANGLT,DVNOR ,
+ > MATALB,LMIX )
+ ELSE IF(ABS(IPLP(1,IPLOT)) .EQ. 3) THEN
+*----
+* PLANA
+* Test if plane is valid
+*----
+ IVALID=TLMVPL(NDIM,NANGL,NPLOTS,IPLOT,IPLP,DPLP,DANGLT,XYZL)
+ IF(IVALID . GE. 0) THEN
+ CALL TLMPLA(IPMAT ,IFTRK ,IPRINT,NSKTRK,NBTR ,NDIM ,
+ > NREG ,MXSUB ,MXSEG ,NANGL ,NBDR ,NPLOTS,
+ > IPLOT ,IPLP ,DPLP ,DANGLT,DVNOR )
+ ENDIF
+ ELSE IF(ABS(IPLP(1,IPLOT)) .EQ. 4) THEN
+*----
+* PLANP
+* Test if plane is valid
+*----
+ IVALID=TLMVPL(NDIM,NANGL,NPLOTS,IPLOT,IPLP,
+ > DPLP,DANGLT,XYZL)
+ IF(IVALID . GE. 0) THEN
+ CALL TLMPLP(IPMAT ,IFTRK ,IPRINT,NSKTRK,NBTR ,NDIM ,
+ > NREG ,MXSUB ,MXSEG ,NANGL ,NBDR ,NPLOTS,
+ > IPLOT ,IPLP ,DPLP ,DANGLT,DVNOR )
+ ENDIF
+ ELSE IF(ABS(IPLP(1,IPLOT)) .EQ. 5) THEN
+*----
+* REGION
+*----
+ CALL TLMREG(IPMAT ,IFTRK ,IPRINT,NSKTRK,NBTR ,NDIM ,
+ > NSOUT ,NREG ,MXSUB ,MXSEG ,NANGL ,NBDR ,
+ > NPLOTS,IPLOT ,IPLP ,DANGLT,DVNOR ,MATALB,
+ > LMIX )
+ ENDIF
+ ENDDO
+*----
+* Release memory
+*----
+ DEALLOCATE(DPLP,IPLP,MATALB,DVNOR,DANGLT,DGMESH)
+*----
+* Processing finished, return
+*----
+ RETURN
+*----
+* Matlab .m file format
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows '/
+ > ' NREG =',I10/
+ > ' NSUR =',I10/
+ > ' NANGL=',I10/
+ > ' NDIM =',I10/
+ > ' ITGEO=',I10,10X,'*)')
+ 7000 FORMAT('%'/
+ > '% File generated using : ',6X,A6/
+ > '% Tracking structure name : ',A12/
+ > '% Tracking file name : ',A12/
+ > '% Title : ',A72/
+ > '%')
+ 7010 FORMAT('pause ;')
+ 9000 FORMAT('NREG =',2I10/
+ > 'NSUR =',2I10/
+ > 'NANGL=',2I10/
+ > 'NDIM =',2I10)
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:33:54 +0000