| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25
26 /// \file materials/src/G4LatticeLogical.cc
27 /// \brief Implementation of the G4LatticeLogical class
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, G4int pRes, G4int polarizationState, G4String map)
54 {
55 if (tRes > MAXRES || pRes > MAXRES) {
56 G4cerr << "G4LatticeLogical::LoadMap exceeds maximum resolution of " << MAXRES << " by "
57 << MAXRES << ". terminating." << G4endl;
58 return false; // terminate if resolution out of bounds.
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] = vgrp * (m / s);
71 }
72 }
73
74 if (verboseLevel != 0) {
75 G4cout << "\nG4LatticeLogical::LoadMap(" << map << ") successful"
76 << " (Vg scalars " << tRes << " x " << pRes << " for polarization " << polarizationState
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, G4int pRes, G4int polarizationState, G4String map)
90 {
91 if (tRes > MAXRES || pRes > MAXRES) {
92 G4cerr << "G4LatticeLogical::LoadMap exceeds maximum resolution of " << MAXRES << " by "
93 << MAXRES << ". terminating." << G4endl;
94 return false; // terminate if resolution out of bounds.
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 coordinates from file
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] = dir.unit(); // Enforce unity
109 }
110 }
111
112 if (verboseLevel != 0) {
113 G4cout << "\nG4LatticeLogical::Load_NMap(" << map << ") successful"
114 << " (Vdir " << tRes << " x " << pRes << " for polarization " << polarizationState
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 physical volume Vol
125 // and polarizationState(0=LON, 1=FT, 2=ST),
126 // returns phonon velocity in m/s
127
128 G4double G4LatticeLogical::MapKtoV(G4int polarizationState, const G4ThreeVector& k) const
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(theta / tRes)][int(phi / pRes)];
146
147 if (Vg == 0) {
148 G4cout << "\nFound v=0 for polarization " << polarizationState << " theta " << theta << " phi "
149 << phi << " translating to map coords "
150 << "theta " << int(theta / tRes) << " phi " << int(phi / pRes) << G4endl;
151 }
152
153 if (verboseLevel > 1) {
154 G4cout << "G4LatticeLogical::MapKtoV theta,phi=" << theta << " " << phi << " : ith,iph "
155 << int(theta / tRes) << " " << int(phi / pRes) << " : V " << Vg << G4endl;
156 }
157
158 return Vg;
159 }
160
161
162 // Given the phonon wave vector k, phonon physical volume Vol
163 // and polarizationState(0=LON, 1=FT, 2=ST),
164 // returns phonon propagation direction as dimensionless unit vector
165
166 G4ThreeVector G4LatticeLogical::MapKtoVDir(G4int polarizationState, const G4ThreeVector& k) const
167 {
168 G4double theta, phi, tRes, pRes;
169
170 tRes = pi / (fDresTheta - 1); // The summant "-1" is required:index=[0:array length-1]
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 //if(phi>pi/2) phi=phi-pi/2;
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 theta,phi=" << theta << " " << phi << " : ith,iph "
189 << iTheta << " " << iPhi << " : dir " << fN_map[polarizationState][iTheta][iPhi]
190 << G4endl;
191 }
192
193 return fN_map[polarizationState][iTheta][iPhi];
194 }
195
196
197 // Dump structure in format compatible with reading back
198
199 void G4LatticeLogical::Dump(std::ostream& os) const
200 {
201 os << "dyn " << fBeta << " " << fGamma << " " << fLambda << " " << fMu << "\nscat " << fB
202 << " decay " << fA << "\nLDOS " << fLDOS << " STDOS " << fSTDOS << " FTDOS " << fFTDOS
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& os, G4int pol, const G4String& name) const
215 {
216 os << "VG " << name << " "
217 << (pol == 0 ? "L"
218 : pol == 1 ? "FT"
219 : pol == 2 ? "ST"
220 : "??")
221 << " " << fVresTheta << " " << fVresPhi << std::endl;
222
223 for (G4int iTheta = 0; iTheta < fVresTheta; iTheta++) {
224 for (G4int iPhi = 0; iPhi < fVresPhi; iPhi++) {
225 os << fMap[pol][iTheta][iPhi] << std::endl;
226 }
227 }
228 }
229
230 void G4LatticeLogical::Dump_NMap(std::ostream& os, G4int pol, const G4String& name) const
231 {
232 os << "VDir " << name << " "
233 << (pol == 0 ? "L"
234 : pol == 1 ? "FT"
235 : pol == 2 ? "ST"
236 : "??")
237 << " " << fDresTheta << " " << fDresPhi << std::endl;
238
239 for (G4int iTheta = 0; iTheta < fDresTheta; iTheta++) {
240 for (G4int iPhi = 0; iPhi < fDresPhi; iPhi++) {
241 os << fN_map[pol][iTheta][iPhi].x() << " " << fN_map[pol][iTheta][iPhi].y() << " "
242 << fN_map[pol][iTheta][iPhi].z() << std::endl;
243 }
244 }
245 }
246