1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
|
*DECK TLMGEO
SUBROUTINE TLMGEO(IPTRK,IPMAT,IPRINT,ITGEO,MAXMDH,NDIM,NUCELL,
> DGMESH,XYZL)
*
*-----------------------------------------------------------------------
*
*Purpose:
* To generate the Matlab instruction for drawing the global geometry.
*
*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):
* C. Plamondon, G. Marleau
*
*Parameters: input
* IPTRK tracking data structure.
* IPMAT pointer to Matlab-m file.
* IPRINT print level.
* ITGEO type of geometry (0 for annular; 1 for Cartesian;
* 2 for hexagonal).
* MAXMDH maximum dimensions of DGMESH.
* NDIM number of dimensions for problem.
* NUCELL cell dimensions for each direction.
* DGMESH mesh of global beometry.
*
*Parameters: output
* XYZL mesh limits.
*
*----------
*
USE GANLIB
IMPLICIT NONE
*----
* Subroutine arguments
*----
TYPE(C_PTR) IPTRK
INTEGER IPMAT
INTEGER IPRINT,ITGEO,MAXMDH,NDIM
INTEGER NUCELL(NDIM)
DOUBLE PRECISION DGMESH(-1:MAXMDH,4),XYZL(2,3)
*----
* Local parameters
*----
INTEGER IOUT
CHARACTER NAMSBR*6
PARAMETER (IOUT=6,NAMSBR='TLMGEO')
DOUBLE PRECISION DZERO
PARAMETER (DZERO=0.0D0)
*----
* Other local variables
*----
INTEGER IDIR,NX,NY,NZ,IPT,IH
CHARACTER NAMREC*12
DOUBLE PRECISION POSHRX(7),POSHRY(7),POSHDX(7),POSHDY(7),
> ROTMAT(2,2)
*----
* Data
*----
CHARACTER CDIR(4)*1
SAVE CDIR
DATA CDIR /'X','Y','Z','R'/
*----
* Processing starts:
* print routine opening header if required
* and initialize various parameters.
*----
IF(IPRINT .GE. 1) THEN
WRITE(IOUT,6000) NAMSBR
ENDIF
DO IDIR=1,NDIM
NAMREC='G00000001SM'//CDIR(IDIR)
IF(ITGEO .EQ. 2) THEN
CALL LCMGET(IPTRK,NAMREC,DGMESH(0,IDIR))
ELSE
CALL LCMGET(IPTRK,NAMREC,DGMESH(0,IDIR))
XYZL(1,IDIR)=DGMESH(0,IDIR)
XYZL(2,IDIR)=DGMESH(NUCELL(IDIR),IDIR)
ENDIF
ENDDO
WRITE(IPMAT,7000)
NX=NUCELL(1)
NY=NUCELL(2)
* write(6,*) 'ITGEO,NDIM=',ITGEO,NDIM,NX,NY
* write(6,*) 'DGMESHx=',(DGMESH(IH,1),IH=-1,NX)
* write(6,*) 'DGMESHy=',(DGMESH(IH,2),IH=-1,NX)
*----
* For hexagon, find reference corner positions
*----
IF(ABS(ITGEO) .EQ. 2) THEN
* IF(ITGEO .EQ. 2) THEN
*----
* One side parallel to x-axis
*----
* POSHRX(1)=DGMESH(0,1)
* POSHRY(1)=DZERO
* ELSE
*----
* One side parallel to y-axis
*----
POSHRX(1)=DZERO
POSHRY(1)=DGMESH(0,2)
* ENDIF
ROTMAT(1,1)=0.5D0
ROTMAT(2,1)=SQRT(3.0D0)/2.0D0
ROTMAT(1,2)=-ROTMAT(2,1)
ROTMAT(2,2)=ROTMAT(1,1)
DO IPT=2,7
POSHRX(IPT)=ROTMAT(1,1)*POSHRX(IPT-1)
> +ROTMAT(1,2)*POSHRY(IPT-1)
POSHRY(IPT)=ROTMAT(2,1)*POSHRX(IPT-1)
> +ROTMAT(2,2)*POSHRY(IPT-1)
ENDDO
WRITE(IPMAT,7040) CDIR(1),NX
WRITE(IPMAT,7040) CDIR(2),NX
IF(NDIM .EQ. 3) THEN
WRITE(IPMAT,7040) CDIR(3),NZ
ENDIF
DO IH=1,NX
DO IPT=1,7
POSHDX(IPT)=POSHRX(IPT)+DGMESH(IH,1)
POSHDY(IPT)=POSHRY(IPT)+DGMESH(IH,2)
ENDDO
WRITE(IPMAT,7041) CDIR(1),IH,(POSHDX(IPT),IPT=1,7)
WRITE(IPMAT,7041) CDIR(2),IH,(POSHDY(IPT),IPT=1,7)
ENDDO
ENDIF
*----
* Print IPMAT header
*----
IF(NDIM .EQ. 2) THEN
IF(ITGEO .EQ. 2) THEN
DO IH=1,NX
WRITE(IPMAT,7042) IH,IH
ENDDO
ELSE
WRITE(IPMAT,7020)
> DGMESH(0,1),DGMESH(NX,1),
> DGMESH(0,2),DGMESH(NY,2)
WRITE(IPMAT,7021)
ENDIF
ELSE IF(NDIM .EQ. 3) THEN
NZ=NUCELL(3)
WRITE(IPMAT,7001)
IF(ITGEO .EQ. 2) THEN
WRITE(IPMAT,7043)
> DGMESH(1,3),DGMESH(NZ,3)
DO IH=1,NX
WRITE(IPMAT,7042) IH,IH
ENDDO
ELSE
WRITE(IPMAT,7030)
> DGMESH(1,1),DGMESH(NX,1),
> DGMESH(1,2),DGMESH(NY,2),
> DGMESH(1,3),DGMESH(NZ,3)
WRITE(IPMAT,7031)
ENDIF
ENDIF
*----
* Processing finished, return
*----
IF(IPRINT .GE. 1) THEN
WRITE(IOUT,6001) NAMSBR
ENDIF
RETURN
*----
* Output formats
*----
6000 FORMAT('(* Output from --',A6,'-- follows ')
6001 FORMAT(' Output from --',A6,'-- completed *)')
*----
* Matlab .m file format
*----
7000 FORMAT('figure;'/'hold on;'/'axis equal;'/
>12Hxlabel('x');/12Hylabel('y');)
7001 FORMAT(12Hzlabel('z');)
7020 FORMAT('xmin=',F18.10,';'/'xmax=',F18.10,';'/
> 'ymin=',F18.10,';'/'ymax=',F18.10,';'/)
7021 FORMAT('line([xmin,xmin],[ymin,ymax]);'/
> 'line([xmin,xmax],[ymin,ymin]);'/
> 'line([xmax,xmax],[ymin,ymax]);'/
> 'line([xmin,xmax],[ymax,ymax]);')
7030 FORMAT('xmin=',F18.10,';'/'xmax=',F18.10,';'/
> 'ymin=',F18.10,';'/'ymax=',F18.10,';'/
> 'zmin=',F18.10,';'/'zmax=',F18.10,';'/)
7031 FORMAT('line([xmin,xmin],[ymin,ymax],[zmin,zmin]);'/
> 'line([xmin,xmax],[ymin,ymin],[zmin,zmin]);'/
> 'line([xmax,xmax],[ymin,ymax],[zmin,zmin]);'/
> 'line([xmin,xmax],[ymax,ymax],[zmin,zmin]);'/
> 'line([xmin,xmin],[ymin,ymin],[zmin,zmax]);'/
> 'line([xmax,xmax],[ymin,ymin],[zmin,zmax]);'/
> 'line([xmax,xmax],[ymax,ymax],[zmin,zmax]);'/
> 'line([xmin,xmin],[ymax,ymax],[zmin,zmax]);'/
> 'line([xmin,xmin],[ymin,ymax],[zmax,zmax]);'/
> 'line([xmin,xmax],[ymin,ymin],[zmax,zmax]);'/
> 'line([xmax,xmax],[ymin,ymax],[zmax,zmax]);'/
> 'line([xmin,xmax],[ymax,ymax],[zmax,zmax]);')
7040 FORMAT('HexM',A1,'=zeros(',I4,',7);')
7041 FORMAT('HexM',A1,'(',I4,',:)=[',7F18.10,'];')
7042 FORMAT('line(HexMX(',I4,',:),HexMY(',I4,',:));')
7043 FORMAT('zmin=',F18.10,';'/'zmax=',F18.10,';'/)
END
|
