| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25
26 /// \file materials/src/G4LatticeLogical.cc
27 /// \brief Implementation of the G4LatticeLogi
28
29 #include "G4LatticeLogical.hh"
30
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33
34 #include <cmath>
35 #include <fstream>
36
37
38 G4LatticeLogical::G4LatticeLogical()
39 {
40 for (G4int i = 0; i < 3; i++) {
41 for (G4int j = 0; j < MAXRES; j++) {
42 for (G4int k = 0; k < MAXRES; k++) {
43 fMap[i][j][k] = 0.;
44 fN_map[i][j][k].set(0., 0., 0.);
45 }
46 }
47 }
48 }
49
50 ///////////////////////////////////////////
51 // Load map of group velocity scalars (m/s)
52 ////////////////////////////////////////////
53 G4bool G4LatticeLogical::LoadMap(G4int tRes, G
54 {
55 if (tRes > MAXRES || pRes > MAXRES) {
56 G4cerr << "G4LatticeLogical::LoadMap excee
57 << MAXRES << ". terminating." << G4
58 return false; // terminate if resolution
59 }
60
61 std::ifstream fMapFile(map.data());
62 if (! fMapFile.is_open()) {
63 return false;
64 }
65
66 G4double vgrp = 0.;
67 for (G4int theta = 0; theta < tRes; theta++)
68 for (G4int phi = 0; phi < pRes; phi++) {
69 fMapFile >> vgrp;
70 fMap[polarizationState][theta][phi] = vg
71 }
72 }
73
74 if (verboseLevel != 0) {
75 G4cout << "\nG4LatticeLogical::LoadMap(" <
76 << " (Vg scalars " << tRes << " x "
77 << ")." << G4endl;
78 }
79
80 fVresTheta = tRes; // store map dimensions
81 fVresPhi = pRes;
82 return true;
83 }
84
85
86 ////////////////////////////////////
87 // Load map of group velocity unit vectors
88 ///////////////////////////////////
89 G4bool G4LatticeLogical::Load_NMap(G4int tRes,
90 {
91 if (tRes > MAXRES || pRes > MAXRES) {
92 G4cerr << "G4LatticeLogical::LoadMap excee
93 << MAXRES << ". terminating." << G4
94 return false; // terminate if resolution
95 }
96
97 std::ifstream fMapFile(map.data());
98 if (! fMapFile.is_open()) {
99 return false;
100 }
101
102 G4double x, y, z; // Buffers to read coordi
103 G4ThreeVector dir;
104 for (G4int theta = 0; theta < tRes; theta++)
105 for (G4int phi = 0; phi < pRes; phi++) {
106 fMapFile >> x >> y >> z;
107 dir.set(x, y, z);
108 fN_map[polarizationState][theta][phi] =
109 }
110 }
111
112 if (verboseLevel != 0) {
113 G4cout << "\nG4LatticeLogical::Load_NMap("
114 << " (Vdir " << tRes << " x " << pR
115 << ")." << G4endl;
116 }
117
118 fDresTheta = tRes; // store map dimensions
119 fDresPhi = pRes;
120 return true;
121 }
122
123
124 // Given the phonon wave vector k, phonon phys
125 // and polarizationState(0=LON, 1=FT, 2=ST),
126 // returns phonon velocity in m/s
127
128 G4double G4LatticeLogical::MapKtoV(G4int polar
129 {
130 G4double theta, phi, tRes, pRes;
131
132 tRes = pi / fVresTheta;
133 pRes = twopi / fVresPhi;
134
135 theta = k.getTheta();
136 phi = k.getPhi();
137
138 if (phi < 0) {
139 phi = phi + twopi;
140 }
141 if (theta > pi) {
142 theta = theta - pi;
143 }
144
145 G4double Vg = fMap[polarizationState][int(th
146
147 if (Vg == 0) {
148 G4cout << "\nFound v=0 for polarization "
149 << phi << " translating to map coor
150 << "theta " << int(theta / tRes) <<
151 }
152
153 if (verboseLevel > 1) {
154 G4cout << "G4LatticeLogical::MapKtoV theta
155 << int(theta / tRes) << " " << int(
156 }
157
158 return Vg;
159 }
160
161
162 // Given the phonon wave vector k, phonon phys
163 // and polarizationState(0=LON, 1=FT, 2=ST),
164 // returns phonon propagation direction as dim
165
166 G4ThreeVector G4LatticeLogical::MapKtoVDir(G4i
167 {
168 G4double theta, phi, tRes, pRes;
169
170 tRes = pi / (fDresTheta - 1); // The summan
171 pRes = 2 * pi / (fDresPhi - 1);
172
173 theta = k.getTheta();
174 phi = k.getPhi();
175
176 if (theta > pi) {
177 theta = theta - pi;
178 }
179 // phi=[0 to 2 pi] in accordance with DMC //
180 if (phi < 0) {
181 phi = phi + 2 * pi;
182 }
183
184 auto iTheta = int(theta / tRes + 0.5);
185 auto iPhi = int(phi / pRes + 0.5);
186
187 if (verboseLevel > 1) {
188 G4cout << "G4LatticeLogical::MapKtoVDir th
189 << iTheta << " " << iPhi << " : dir
190 << G4endl;
191 }
192
193 return fN_map[polarizationState][iTheta][iPh
194 }
195
196
197 // Dump structure in format compatible with re
198
199 void G4LatticeLogical::Dump(std::ostream& os)
200 {
201 os << "dyn " << fBeta << " " << fGamma << "
202 << " decay " << fA << "\nLDOS " << fLDOS
203 << std::endl;
204
205 Dump_NMap(os, 0, "LVec.ssv");
206 Dump_NMap(os, 1, "FTVec.ssv");
207 Dump_NMap(os, 2, "STVec.ssv");
208
209 DumpMap(os, 0, "L.ssv");
210 DumpMap(os, 1, "FT.ssv");
211 DumpMap(os, 2, "ST.ssv");
212 }
213
214 void G4LatticeLogical::DumpMap(std::ostream& o
215 {
216 os << "VG " << name << " "
217 << (pol == 0 ? "L"
218 : pol == 1 ? "FT"
219 : pol == 2 ? "ST"
220 : "??")
221 << " " << fVresTheta << " " << fVresPhi <
222
223 for (G4int iTheta = 0; iTheta < fVresTheta;
224 for (G4int iPhi = 0; iPhi < fVresPhi; iPhi
225 os << fMap[pol][iTheta][iPhi] << std::en
226 }
227 }
228 }
229
230 void G4LatticeLogical::Dump_NMap(std::ostream&
231 {
232 os << "VDir " << name << " "
233 << (pol == 0 ? "L"
234 : pol == 1 ? "FT"
235 : pol == 2 ? "ST"
236 : "??")
237 << " " << fDresTheta << " " << fDresPhi <
238
239 for (G4int iTheta = 0; iTheta < fDresTheta;
240 for (G4int iPhi = 0; iPhi < fDresPhi; iPhi
241 os << fN_map[pol][iTheta][iPhi].x() << "
242 << fN_map[pol][iTheta][iPhi].z() << s
243 }
244 }
245 }
246