| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 ////////////////////////////////////////////// 26 ////////////////////////////////////////////////////////////////////////
27 // Optical Photon Boundary Process Class Imple 27 // Optical Photon Boundary Process Class Implementation
28 ////////////////////////////////////////////// 28 ////////////////////////////////////////////////////////////////////////
29 // 29 //
30 // File: G4OpBoundaryProcess.cc 30 // File: G4OpBoundaryProcess.cc
31 // Description: Discrete Process -- reflection 31 // Description: Discrete Process -- reflection/refraction at
32 // optical in 32 // optical interfaces
33 // Version: 1.1 33 // Version: 1.1
34 // Created: 1997-06-18 34 // Created: 1997-06-18
35 // Modified: 1998-05-25 - Correct parallel 35 // Modified: 1998-05-25 - Correct parallel component of polarization
36 // (thanks to: Stefa 36 // (thanks to: Stefano Magni + Giovanni Pieri)
37 // 1998-05-28 - NULL Rindex point 37 // 1998-05-28 - NULL Rindex pointer before reuse
38 // (thanks to: Stefa 38 // (thanks to: Stefano Magni)
39 // 1998-06-11 - delete *sint1 in 39 // 1998-06-11 - delete *sint1 in oblique reflection
40 // (thanks to: Giova 40 // (thanks to: Giovanni Pieri)
41 // 1998-06-19 - move from GetLoca << 41 // 1998-06-19 - move from GetLocalExitNormal() to the new
42 // method: GetLocalE 42 // method: GetLocalExitNormal(&valid) to get
43 // the surface norma 43 // the surface normal in all cases
44 // 1998-11-07 - NULL OpticalSurfa 44 // 1998-11-07 - NULL OpticalSurface pointer before use
45 // comparison not sh 45 // comparison not sharp for: std::abs(cost1) < 1.0
46 // remove sin1, sin2 46 // remove sin1, sin2 in lines 556,567
47 // (thanks to Stefan 47 // (thanks to Stefano Magni)
48 // 1999-10-10 - Accommodate chang 48 // 1999-10-10 - Accommodate changes done in DoAbsorption by
49 // changing logic in 49 // changing logic in DielectricMetal
50 // 2001-10-18 - avoid Linux (gcc- 50 // 2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
51 // might be used uni 51 // might be used uninitialized in this function
52 // moved E2_perp, E2 52 // moved E2_perp, E2_parl and E2_total out of 'if'
53 // 2003-11-27 - Modified line 168 53 // 2003-11-27 - Modified line 168-9 to reflect changes made to
54 // G4OpticalSurface 54 // G4OpticalSurface class ( by Fan Lei)
55 // 2004-02-02 - Set theStatus = U 55 // 2004-02-02 - Set theStatus = Undefined at start of DoIt
56 // 2005-07-28 - add G4ProcessType 56 // 2005-07-28 - add G4ProcessType to constructor
57 // 2006-11-04 - add capability of 57 // 2006-11-04 - add capability of calculating the reflectivity
58 // off a metal surfa << 58 // off a metal surface by way of a complex index
59 // of refraction - T << 59 // of refraction - Thanks to Sehwook Lee and John
60 // Hauptman (Dept. o 60 // Hauptman (Dept. of Physics - Iowa State Univ.)
61 // 2009-11-10 - add capability of 61 // 2009-11-10 - add capability of simulating surface reflections
62 // with Look-Up-Tabl 62 // with Look-Up-Tables (LUT) containing measured
63 // optical reflectan 63 // optical reflectance for a variety of surface
64 // treatments - Than 64 // treatments - Thanks to Martin Janecek and
65 // William Moses (La 65 // William Moses (Lawrence Berkeley National Lab.)
66 // 2013-06-01 - add the capabilit 66 // 2013-06-01 - add the capability of simulating the transmission
67 // of a dichronic fi 67 // of a dichronic filter
68 // 2017-02-24 - add capability of <<
69 // with Look-Up-Tabl <<
70 // 68 //
71 // Author: Peter Gumplinger 69 // Author: Peter Gumplinger
72 // adopted from work by Werner Keil - April 70 // adopted from work by Werner Keil - April 2/96
>> 71 // mail: gum@triumf.ca
73 // 72 //
74 ////////////////////////////////////////////// 73 ////////////////////////////////////////////////////////////////////////
75 74
76 #include "G4OpBoundaryProcess.hh" <<
77 <<
78 #include "G4ios.hh" 75 #include "G4ios.hh"
79 #include "G4GeometryTolerance.hh" << 76 #include "G4PhysicalConstants.hh"
80 #include "G4LogicalBorderSurface.hh" <<
81 #include "G4LogicalSkinSurface.hh" <<
82 #include "G4OpProcessSubType.hh" 77 #include "G4OpProcessSubType.hh"
83 #include "G4OpticalParameters.hh" << 78
>> 79 #include "G4OpBoundaryProcess.hh"
>> 80 #include "G4GeometryTolerance.hh"
>> 81
>> 82 #include "G4VSensitiveDetector.hh"
84 #include "G4ParallelWorldProcess.hh" 83 #include "G4ParallelWorldProcess.hh"
85 #include "G4PhysicalConstants.hh" << 84
86 #include "G4SystemOfUnits.hh" 85 #include "G4SystemOfUnits.hh"
87 #include "G4TransportationManager.hh" <<
88 #include "G4VSensitiveDetector.hh" <<
89 86
90 //....oooOO0OOooo........oooOO0OOooo........oo << 87 /////////////////////////
>> 88 // Class Implementation
>> 89 /////////////////////////
>> 90
>> 91 //////////////
>> 92 // Operators
>> 93 //////////////
>> 94
>> 95 // G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
>> 96 // {
>> 97 // }
>> 98
>> 99 /////////////////
>> 100 // Constructors
>> 101 /////////////////
>> 102
91 G4OpBoundaryProcess::G4OpBoundaryProcess(const 103 G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
92 G4Pro << 104 G4ProcessType type)
93 : G4VDiscreteProcess(processName, ptype) << 105 : G4VDiscreteProcess(processName, type)
94 { 106 {
95 Initialise(); << 107 if ( verboseLevel > 0) {
>> 108 G4cout << GetProcessName() << " is created " << G4endl;
>> 109 }
96 110
97 if(verboseLevel > 0) << 111 SetProcessSubType(fOpBoundary);
98 { <<
99 G4cout << GetProcessName() << " is created <<
100 } <<
101 SetProcessSubType(fOpBoundary); <<
102 112
103 fStatus = Undefined; << 113 theStatus = Undefined;
104 fModel = glisur; << 114 theModel = glisur;
105 fFinish = polished; << 115 theFinish = polished;
106 fReflectivity = 1.; << 116 theReflectivity = 1.;
107 fEfficiency = 0.; << 117 theEfficiency = 0.;
108 fTransmittance = 0.; << 118 theTransmittance = 0.;
109 fSurfaceRoughness = 0.; <<
110 fProb_sl = 0.; <<
111 fProb_ss = 0.; <<
112 fProb_bs = 0.; <<
113 <<
114 fRealRIndexMPV = nullptr; <<
115 fImagRIndexMPV = nullptr; <<
116 fMaterial1 = nullptr; <<
117 fMaterial2 = nullptr; <<
118 fOpticalSurface = nullptr; <<
119 fCarTolerance = G4GeometryTolerance::GetIn <<
120 <<
121 f_iTE = f_iTM = 0; <<
122 fPhotonMomentum = 0.; <<
123 fRindex1 = fRindex2 = 1.; <<
124 fSint1 = 0.; <<
125 fDichroicVector = nullptr; <<
126 } <<
127 119
128 //....oooOO0OOooo........oooOO0OOooo........oo << 120 theSurfaceRoughness = 0.;
129 G4OpBoundaryProcess::~G4OpBoundaryProcess() = <<
130 121
131 //....oooOO0OOooo........oooOO0OOooo........oo << 122 prob_sl = 0.;
132 void G4OpBoundaryProcess::PreparePhysicsTable( << 123 prob_ss = 0.;
133 { << 124 prob_bs = 0.;
134 Initialise(); <<
135 } <<
136 125
137 //....oooOO0OOooo........oooOO0OOooo........oo << 126 PropertyPointer = NULL;
138 void G4OpBoundaryProcess::Initialise() << 127 PropertyPointer1 = NULL;
139 { << 128 PropertyPointer2 = NULL;
140 G4OpticalParameters* params = G4OpticalParam <<
141 SetInvokeSD(params->GetBoundaryInvokeSD()); <<
142 SetVerboseLevel(params->GetBoundaryVerboseLe <<
143 } <<
144 129
145 //....oooOO0OOooo........oooOO0OOooo........oo << 130 Material1 = NULL;
146 G4VParticleChange* G4OpBoundaryProcess::PostSt << 131 Material2 = NULL;
147 <<
148 { <<
149 fStatus = Undefined; <<
150 aParticleChange.Initialize(aTrack); <<
151 aParticleChange.ProposeVelocity(aTrack.GetVe <<
152 <<
153 // Get hyperStep from G4ParallelWorldProces <<
154 // NOTE: PostSetpDoIt of this process to be <<
155 // G4ParallelWorldProcess! <<
156 const G4Step* pStep = &aStep; <<
157 const G4Step* hStep = G4ParallelWorldProcess <<
158 if(hStep != nullptr) <<
159 pStep = hStep; <<
160 <<
161 if(pStep->GetPostStepPoint()->GetStepStatus( <<
162 { <<
163 fMaterial1 = pStep->GetPreStepPoint()->Get <<
164 fMaterial2 = pStep->GetPostStepPoint()->Ge <<
165 } <<
166 else <<
167 { <<
168 fStatus = NotAtBoundary; <<
169 if(verboseLevel > 1) <<
170 BoundaryProcessVerbose(); <<
171 return G4VDiscreteProcess::PostStepDoIt(aT <<
172 } <<
173 132
174 G4VPhysicalVolume* thePrePV = pStep->GetPre << 133 OpticalSurface = NULL;
175 G4VPhysicalVolume* thePostPV = pStep->GetPos <<
176 134
177 if(verboseLevel > 1) << 135 kCarTolerance = G4GeometryTolerance::GetInstance()
178 { << 136 ->GetSurfaceTolerance();
179 G4cout << " Photon at Boundary! " << G4end <<
180 if(thePrePV != nullptr) <<
181 G4cout << " thePrePV: " << thePrePV->Ge <<
182 if(thePostPV != nullptr) <<
183 G4cout << " thePostPV: " << thePostPV->G <<
184 } <<
185 137
186 G4double stepLength = aTrack.GetStepLength() << 138 iTE = iTM = 0;
187 if(stepLength <= fCarTolerance) << 139 thePhotonMomentum = 0.;
188 { << 140 Rindex1 = Rindex2 = 1.;
189 fStatus = StepTooSmall; << 141 cost1 = cost2 = sint1 = sint2 = 0.;
190 if(verboseLevel > 1) <<
191 BoundaryProcessVerbose(); <<
192 <<
193 G4MaterialPropertyVector* groupvel = nullp <<
194 G4MaterialPropertiesTable* aMPT = fMateria <<
195 if(aMPT != nullptr) <<
196 { <<
197 groupvel = aMPT->GetProperty(kGROUPVEL); <<
198 } <<
199 <<
200 if(groupvel != nullptr) <<
201 { <<
202 aParticleChange.ProposeVelocity( <<
203 groupvel->Value(fPhotonMomentum, idx_g <<
204 } <<
205 return G4VDiscreteProcess::PostStepDoIt(aT <<
206 } <<
207 else if (stepLength <= 10.*fCarTolerance && <<
208 { // see bug 2510 <<
209 ++fNumSmallStepWarnings; <<
210 if(verboseLevel > 0) <<
211 { <<
212 G4ExceptionDescription ed; <<
213 ed << "G4OpBoundaryProcess: " <<
214 << "Opticalphoton step length: " << s <<
215 << "This is larger than the threshold <<
216 "to set status StepTooSmall." << G <<
217 << "Boundary scattering may be incorr <<
218 if(fNumSmallStepWarnings == 10) <<
219 { <<
220 ed << G4endl << "*** Step size warning <<
221 } <<
222 G4Exception("G4OpBoundaryProcess", "OpBo <<
223 } <<
224 } <<
225 142
226 const G4DynamicParticle* aParticle = aTrack. << 143 idx = idy = 0;
>> 144 DichroicVector = NULL;
227 145
228 fPhotonMomentum = aParticle->GetTotalMoment << 146 fInvokeSD = true;
229 fOldMomentum = aParticle->GetMomentumDir << 147 }
230 fOldPolarization = aParticle->GetPolarizatio <<
231 <<
232 if(verboseLevel > 1) <<
233 { <<
234 G4cout << " Old Momentum Direction: " << f <<
235 << " Old Polarization: " << f <<
236 } <<
237 148
238 G4ThreeVector theGlobalPoint = pStep->GetPos << 149 // G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
239 G4bool valid; << 150 // {
>> 151 // }
>> 152
>> 153 ////////////////
>> 154 // Destructors
>> 155 ////////////////
>> 156
>> 157 G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
>> 158
>> 159 ////////////
>> 160 // Methods
>> 161 ////////////
240 162
241 // ID of Navigator which limits step << 163 // PostStepDoIt
242 G4int hNavId = G4ParallelWorldProcess::GetHy << 164 // ------------
243 auto iNav = G4TransportationManager::GetT << 165 //
244 ->GetActiveNavigatorsIterator( <<
245 fGlobalNormal = (iNav[hNavId])->GetGlobalExi <<
246 <<
247 if(valid) <<
248 { <<
249 fGlobalNormal = -fGlobalNormal; <<
250 } <<
251 else <<
252 { <<
253 G4ExceptionDescription ed; <<
254 ed << " G4OpBoundaryProcess/PostStepDoIt() <<
255 << " The Navigator reports that it retu <<
256 G4Exception( <<
257 "G4OpBoundaryProcess::PostStepDoIt", "Op <<
258 "Invalid Surface Normal - Geometry must <<
259 } <<
260 166
261 if(fOldMomentum * fGlobalNormal > 0.0) << 167 G4VParticleChange*
262 { << 168 G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
263 #ifdef G4OPTICAL_DEBUG << 169 {
264 G4ExceptionDescription ed; << 170 theStatus = Undefined;
265 ed << " G4OpBoundaryProcess/PostStepDoIt() <<
266 "wrong direction. " <<
267 << G4endl <<
268 << " The momentum of the photon arriv <<
269 << " must exit the volume cross in th <<
270 << " So it MUST have dot < 0 with the <<
271 "volume (globalNormal)." <<
272 << G4endl << " >> The dot product of <<
273 << fOldMomentum * fGlobalNormal << G4en <<
274 << " Old Momentum (during step) <<
275 << " Global Normal (Exiting New Vol <<
276 << G4endl; <<
277 G4Exception("G4OpBoundaryProcess::PostStep <<
278 EventMustBeAborted, // Or Jus <<
279 // repeat <<
280 ed, <<
281 "Invalid Surface Normal - Geom <<
282 "normal pointing in the right <<
283 #else <<
284 fGlobalNormal = -fGlobalNormal; <<
285 #endif <<
286 } <<
287 171
288 G4MaterialPropertyVector* rIndexMPV = nullpt << 172 aParticleChange.Initialize(aTrack);
289 G4MaterialPropertiesTable* MPT = fMaterial1- << 173 aParticleChange.ProposeVelocity(aTrack.GetVelocity());
290 if(MPT != nullptr) <<
291 { <<
292 rIndexMPV = MPT->GetProperty(kRINDEX); <<
293 } <<
294 if(rIndexMPV != nullptr) <<
295 { <<
296 fRindex1 = rIndexMPV->Value(fPhotonMomentu <<
297 } <<
298 else <<
299 { <<
300 fStatus = NoRINDEX; <<
301 if(verboseLevel > 1) <<
302 BoundaryProcessVerbose(); <<
303 aParticleChange.ProposeLocalEnergyDeposit( <<
304 aParticleChange.ProposeTrackStatus(fStopAn <<
305 return G4VDiscreteProcess::PostStepDoIt(aT <<
306 } <<
307 174
308 fReflectivity = 1.; << 175 // Get hyperStep from G4ParallelWorldProcess
309 fEfficiency = 0.; << 176 // NOTE: PostSetpDoIt of this process should be
310 fTransmittance = 0.; << 177 // invoked after G4ParallelWorldProcess!
311 fSurfaceRoughness = 0.; <<
312 fModel = glisur; <<
313 fFinish = polished; <<
314 G4SurfaceType type = dielectric_dielectric; <<
315 <<
316 rIndexMPV = nullptr; <<
317 fOpticalSurface = nullptr; <<
318 <<
319 G4LogicalSurface* surface = <<
320 G4LogicalBorderSurface::GetSurface(thePreP <<
321 if(surface == nullptr) <<
322 { <<
323 if(thePostPV->GetMotherLogical() == thePre <<
324 { <<
325 surface = G4LogicalSkinSurface::GetSurfa <<
326 if(surface == nullptr) <<
327 { <<
328 surface = <<
329 G4LogicalSkinSurface::GetSurface(the <<
330 } <<
331 } <<
332 else <<
333 { <<
334 surface = G4LogicalSkinSurface::GetSurfa <<
335 if(surface == nullptr) <<
336 { <<
337 surface = <<
338 G4LogicalSkinSurface::GetSurface(the <<
339 } <<
340 } <<
341 } <<
342 178
343 if(surface != nullptr) << 179 const G4Step* pStep = &aStep;
344 { <<
345 fOpticalSurface = <<
346 dynamic_cast<G4OpticalSurface*>(surface- <<
347 } <<
348 if(fOpticalSurface != nullptr) <<
349 { <<
350 type = fOpticalSurface->GetType(); <<
351 fModel = fOpticalSurface->GetModel(); <<
352 fFinish = fOpticalSurface->GetFinish(); <<
353 <<
354 G4MaterialPropertiesTable* sMPT = <<
355 fOpticalSurface->GetMaterialPropertiesTa <<
356 if(sMPT != nullptr) <<
357 { <<
358 if(fFinish == polishedbackpainted || fFi <<
359 { <<
360 rIndexMPV = sMPT->GetProperty(kRINDEX) <<
361 if(rIndexMPV != nullptr) <<
362 { <<
363 fRindex2 = rIndexMPV->Value(fPhotonM <<
364 } <<
365 else <<
366 { <<
367 fStatus = NoRINDEX; <<
368 if(verboseLevel > 1) <<
369 BoundaryProcessVerbose(); <<
370 aParticleChange.ProposeLocalEnergyDe <<
371 aParticleChange.ProposeTrackStatus(f <<
372 return G4VDiscreteProcess::PostStepD <<
373 } <<
374 } <<
375 180
376 fRealRIndexMPV = sMPT->GetProperty(kREAL << 181 const G4Step* hStep = G4ParallelWorldProcess::GetHyperStep();
377 fImagRIndexMPV = sMPT->GetProperty(kIMAG << 182
378 f_iTE = f_iTM = 1; << 183 if (hStep) pStep = hStep;
379 << 184
380 G4MaterialPropertyVector* pp; << 185 G4bool isOnBoundary =
381 if((pp = sMPT->GetProperty(kREFLECTIVITY << 186 (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
382 { << 187
383 fReflectivity = pp->Value(fPhotonMomen << 188 if (isOnBoundary) {
384 } << 189 Material1 = pStep->GetPreStepPoint()->GetMaterial();
385 else if(fRealRIndexMPV && fImagRIndexMPV << 190 Material2 = pStep->GetPostStepPoint()->GetMaterial();
386 { << 191 } else {
387 CalculateReflectivity(); << 192 theStatus = NotAtBoundary;
388 } << 193 if ( verboseLevel > 0) BoundaryProcessVerbose();
389 << 194 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
390 if((pp = sMPT->GetProperty(kEFFICIENCY)) <<
391 { <<
392 fEfficiency = pp->Value(fPhotonMomentu <<
393 } <<
394 if((pp = sMPT->GetProperty(kTRANSMITTANC <<
395 { <<
396 fTransmittance = pp->Value(fPhotonMome <<
397 } <<
398 if(sMPT->ConstPropertyExists(kSURFACEROU <<
399 { <<
400 fSurfaceRoughness = sMPT->GetConstProp <<
401 } <<
402 <<
403 if(fModel == unified) <<
404 { <<
405 fProb_sl = (pp = sMPT->GetProperty(kSP <<
406 ? pp->Value(fPhotonMoment <<
407 : 0.; <<
408 fProb_ss = (pp = sMPT->GetProperty(kSP <<
409 ? pp->Value(fPhotonMoment <<
410 : 0.; <<
411 fProb_bs = (pp = sMPT->GetProperty(kBA <<
412 ? pp->Value(fPhotonMoment <<
413 : 0.; <<
414 } <<
415 } // end of if(sMPT) <<
416 else if(fFinish == polishedbackpainted || <<
417 { <<
418 aParticleChange.ProposeLocalEnergyDeposi <<
419 aParticleChange.ProposeTrackStatus(fStop <<
420 return G4VDiscreteProcess::PostStepDoIt( <<
421 } <<
422 } // end of if(fOpticalSurface) <<
423 <<
424 // DIELECTRIC-DIELECTRIC <<
425 if(type == dielectric_dielectric) <<
426 { <<
427 if(fFinish == polished || fFinish == groun <<
428 { <<
429 if(fMaterial1 == fMaterial2) <<
430 { <<
431 fStatus = SameMaterial; <<
432 if(verboseLevel > 1) <<
433 BoundaryProcessVerbose(); <<
434 return G4VDiscreteProcess::PostStepDoI <<
435 } <<
436 MPT = fMaterial2->GetMaterialPrope <<
437 rIndexMPV = nullptr; <<
438 if(MPT != nullptr) <<
439 { <<
440 rIndexMPV = MPT->GetProperty(kRINDEX); <<
441 } <<
442 if(rIndexMPV != nullptr) <<
443 { <<
444 fRindex2 = rIndexMPV->Value(fPhotonMom <<
445 } <<
446 else <<
447 { <<
448 fStatus = NoRINDEX; <<
449 if(verboseLevel > 1) <<
450 BoundaryProcessVerbose(); <<
451 aParticleChange.ProposeLocalEnergyDepo <<
452 aParticleChange.ProposeTrackStatus(fSt <<
453 return G4VDiscreteProcess::PostStepDoI <<
454 } <<
455 } <<
456 if(fFinish == polishedbackpainted || fFini <<
457 { <<
458 DielectricDielectric(); <<
459 } <<
460 else <<
461 { <<
462 G4double rand = G4UniformRand(); <<
463 if(rand > fReflectivity + fTransmittance <<
464 { <<
465 DoAbsorption(); <<
466 } <<
467 else if(rand > fReflectivity) <<
468 { <<
469 fStatus = Transmission; <<
470 fNewMomentum = fOldMomentum; <<
471 fNewPolarization = fOldPolarization; <<
472 } <<
473 else <<
474 { <<
475 if(fFinish == polishedfrontpainted) <<
476 { <<
477 DoReflection(); <<
478 } <<
479 else if(fFinish == groundfrontpainted) <<
480 { <<
481 fStatus = LambertianReflection; <<
482 DoReflection(); <<
483 } <<
484 else <<
485 { <<
486 DielectricDielectric(); <<
487 } 195 }
488 } << 196
489 } << 197 G4VPhysicalVolume* thePrePV =
490 } << 198 pStep->GetPreStepPoint() ->GetPhysicalVolume();
491 else if(type == dielectric_metal) << 199 G4VPhysicalVolume* thePostPV =
492 { << 200 pStep->GetPostStepPoint()->GetPhysicalVolume();
493 DielectricMetal(); << 201
494 } << 202 if ( verboseLevel > 0 ) {
495 else if(type == dielectric_LUT) << 203 G4cout << " Photon at Boundary! " << G4endl;
496 { << 204 if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl;
497 DielectricLUT(); << 205 if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
498 } <<
499 else if(type == dielectric_LUTDAVIS) <<
500 { <<
501 DielectricLUTDAVIS(); <<
502 } <<
503 else if(type == dielectric_dichroic) <<
504 { <<
505 DielectricDichroic(); <<
506 } <<
507 else if(type == coated) <<
508 { <<
509 CoatedDielectricDielectric(); <<
510 } <<
511 else <<
512 { <<
513 if(fNumBdryTypeWarnings <= 10) <<
514 { <<
515 ++fNumBdryTypeWarnings; <<
516 if(verboseLevel > 0) <<
517 { <<
518 G4ExceptionDescription ed; <<
519 ed << " PostStepDoIt(): Illegal bounda <<
520 if(fNumBdryTypeWarnings == 10) <<
521 { <<
522 ed << "** Boundary type warnings sto <<
523 } 206 }
524 G4Exception("G4OpBoundaryProcess", "Op <<
525 } <<
526 } <<
527 return G4VDiscreteProcess::PostStepDoIt(aT <<
528 } <<
529 207
530 fNewMomentum = fNewMomentum.unit(); << 208 if (aTrack.GetStepLength()<=kCarTolerance/2){
531 fNewPolarization = fNewPolarization.unit(); << 209 theStatus = StepTooSmall;
>> 210 if ( verboseLevel > 0) BoundaryProcessVerbose();
>> 211 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 212 }
532 213
533 if(verboseLevel > 1) << 214 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
534 { <<
535 G4cout << " New Momentum Direction: " << f <<
536 << " New Polarization: " << f <<
537 BoundaryProcessVerbose(); <<
538 } <<
539 215
540 aParticleChange.ProposeMomentumDirection(fNe << 216 thePhotonMomentum = aParticle->GetTotalMomentum();
541 aParticleChange.ProposePolarization(fNewPola << 217 OldMomentum = aParticle->GetMomentumDirection();
>> 218 OldPolarization = aParticle->GetPolarization();
542 219
543 if(fStatus == FresnelRefraction || fStatus = << 220 if ( verboseLevel > 0 ) {
544 { << 221 G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
545 // not all surface types check that fMater << 222 G4cout << " Old Polarization: " << OldPolarization << G4endl;
546 G4MaterialPropertiesTable* aMPT = fMateria << 223 }
547 G4MaterialPropertyVector* groupvel = nullp <<
548 if(aMPT != nullptr) <<
549 { <<
550 groupvel = aMPT->GetProperty(kGROUPVEL); <<
551 } <<
552 if(groupvel != nullptr) <<
553 { <<
554 aParticleChange.ProposeVelocity( <<
555 groupvel->Value(fPhotonMomentum, idx_g <<
556 } <<
557 } <<
558 224
559 if(fStatus == Detection && fInvokeSD) << 225 G4ThreeVector theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
560 InvokeSD(pStep); <<
561 return G4VDiscreteProcess::PostStepDoIt(aTra <<
562 } <<
563 226
564 //....oooOO0OOooo........oooOO0OOooo........oo << 227 G4bool valid;
565 void G4OpBoundaryProcess::BoundaryProcessVerbo << 228 // Use the new method for Exit Normal in global coordinates,
566 { << 229 // which provides the normal more reliably.
567 G4cout << " *** "; <<
568 if(fStatus == Undefined) <<
569 G4cout << "Undefined"; <<
570 else if(fStatus == Transmission) <<
571 G4cout << "Transmission"; <<
572 else if(fStatus == FresnelRefraction) <<
573 G4cout << "FresnelRefraction"; <<
574 else if(fStatus == FresnelReflection) <<
575 G4cout << "FresnelReflection"; <<
576 else if(fStatus == TotalInternalReflection) <<
577 G4cout << "TotalInternalReflection"; <<
578 else if(fStatus == LambertianReflection) <<
579 G4cout << "LambertianReflection"; <<
580 else if(fStatus == LobeReflection) <<
581 G4cout << "LobeReflection"; <<
582 else if(fStatus == SpikeReflection) <<
583 G4cout << "SpikeReflection"; <<
584 else if(fStatus == BackScattering) <<
585 G4cout << "BackScattering"; <<
586 else if(fStatus == PolishedLumirrorAirReflec <<
587 G4cout << "PolishedLumirrorAirReflection"; <<
588 else if(fStatus == PolishedLumirrorGlueRefle <<
589 G4cout << "PolishedLumirrorGlueReflection" <<
590 else if(fStatus == PolishedAirReflection) <<
591 G4cout << "PolishedAirReflection"; <<
592 else if(fStatus == PolishedTeflonAirReflecti <<
593 G4cout << "PolishedTeflonAirReflection"; <<
594 else if(fStatus == PolishedTiOAirReflection) <<
595 G4cout << "PolishedTiOAirReflection"; <<
596 else if(fStatus == PolishedTyvekAirReflectio <<
597 G4cout << "PolishedTyvekAirReflection"; <<
598 else if(fStatus == PolishedVM2000AirReflecti <<
599 G4cout << "PolishedVM2000AirReflection"; <<
600 else if(fStatus == PolishedVM2000GlueReflect <<
601 G4cout << "PolishedVM2000GlueReflection"; <<
602 else if(fStatus == EtchedLumirrorAirReflecti <<
603 G4cout << "EtchedLumirrorAirReflection"; <<
604 else if(fStatus == EtchedLumirrorGlueReflect <<
605 G4cout << "EtchedLumirrorGlueReflection"; <<
606 else if(fStatus == EtchedAirReflection) <<
607 G4cout << "EtchedAirReflection"; <<
608 else if(fStatus == EtchedTeflonAirReflection <<
609 G4cout << "EtchedTeflonAirReflection"; <<
610 else if(fStatus == EtchedTiOAirReflection) <<
611 G4cout << "EtchedTiOAirReflection"; <<
612 else if(fStatus == EtchedTyvekAirReflection) <<
613 G4cout << "EtchedTyvekAirReflection"; <<
614 else if(fStatus == EtchedVM2000AirReflection <<
615 G4cout << "EtchedVM2000AirReflection"; <<
616 else if(fStatus == EtchedVM2000GlueReflectio <<
617 G4cout << "EtchedVM2000GlueReflection"; <<
618 else if(fStatus == GroundLumirrorAirReflecti <<
619 G4cout << "GroundLumirrorAirReflection"; <<
620 else if(fStatus == GroundLumirrorGlueReflect <<
621 G4cout << "GroundLumirrorGlueReflection"; <<
622 else if(fStatus == GroundAirReflection) <<
623 G4cout << "GroundAirReflection"; <<
624 else if(fStatus == GroundTeflonAirReflection <<
625 G4cout << "GroundTeflonAirReflection"; <<
626 else if(fStatus == GroundTiOAirReflection) <<
627 G4cout << "GroundTiOAirReflection"; <<
628 else if(fStatus == GroundTyvekAirReflection) <<
629 G4cout << "GroundTyvekAirReflection"; <<
630 else if(fStatus == GroundVM2000AirReflection <<
631 G4cout << "GroundVM2000AirReflection"; <<
632 else if(fStatus == GroundVM2000GlueReflectio <<
633 G4cout << "GroundVM2000GlueReflection"; <<
634 else if(fStatus == Absorption) <<
635 G4cout << "Absorption"; <<
636 else if(fStatus == Detection) <<
637 G4cout << "Detection"; <<
638 else if(fStatus == NotAtBoundary) <<
639 G4cout << "NotAtBoundary"; <<
640 else if(fStatus == SameMaterial) <<
641 G4cout << "SameMaterial"; <<
642 else if(fStatus == StepTooSmall) <<
643 G4cout << "StepTooSmall"; <<
644 else if(fStatus == NoRINDEX) <<
645 G4cout << "NoRINDEX"; <<
646 else if(fStatus == Dichroic) <<
647 G4cout << "Dichroic Transmission"; <<
648 else if(fStatus == CoatedDielectricReflectio <<
649 G4cout << "Coated Dielectric Reflection"; <<
650 else if(fStatus == CoatedDielectricRefractio <<
651 G4cout << "Coated Dielectric Refraction"; <<
652 else if(fStatus == CoatedDielectricFrustrate <<
653 G4cout << "Coated Dielectric Frustrated Tr <<
654 230
655 G4cout << " ***" << G4endl; << 231 // ID of Navigator which limits step
656 } <<
657 232
658 //....oooOO0OOooo........oooOO0OOooo........oo << 233 G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID();
659 G4ThreeVector G4OpBoundaryProcess::GetFacetNor << 234 std::vector<G4Navigator*>::iterator iNav =
660 const G4ThreeVector& momentum, const G4Three << 235 G4TransportationManager::GetTransportationManager()->
661 { << 236 GetActiveNavigatorsIterator();
662 G4ThreeVector facetNormal; << 237 theGlobalNormal =
663 if(fModel == unified || fModel == LUT || fMo << 238 (iNav[hNavId])->GetGlobalExitNormal(theGlobalPoint,&valid);
664 { <<
665 /* This function codes alpha to a random v <<
666 distribution p(alpha) = g(alpha; 0, sigma_ <<
667 for alpha > 0 and alpha < 90, where g(alph <<
668 gaussian distribution with mean 0 and stan <<
669 <<
670 G4double sigma_alpha = 0.0; <<
671 if(fOpticalSurface) <<
672 sigma_alpha = fOpticalSurface->GetSigmaA <<
673 if(sigma_alpha == 0.0) <<
674 { <<
675 return normal; <<
676 } <<
677 <<
678 G4double f_max = std::min(1.0, 4. * sigma_ <<
679 G4double alpha, phi, sinAlpha; <<
680 <<
681 do <<
682 { // Loop checking, 13-Aug-2015, Peter Gu <<
683 do <<
684 { // Loop checking, 13-Aug-2015, Peter <<
685 alpha = G4RandGauss::shoot(0.0, sig <<
686 sinAlpha = std::sin(alpha); <<
687 } while(G4UniformRand() * f_max > sinAlp <<
688 <<
689 phi = G4UniformRand() * twopi; <<
690 facetNormal.set(sinAlpha * std::cos(phi) <<
691 std::cos(alpha)); <<
692 facetNormal.rotateUz(normal); <<
693 } while(momentum * facetNormal >= 0.0); <<
694 } <<
695 else <<
696 { <<
697 G4double polish = 1.0; <<
698 if(fOpticalSurface) <<
699 polish = fOpticalSurface->GetPolish(); <<
700 if(polish < 1.0) <<
701 { <<
702 do <<
703 { // Loop checking, 13-Aug-2015, Peter <<
704 G4ThreeVector smear; <<
705 do <<
706 { // Loop checking, 13-Aug-2015, Pete <<
707 smear.setX(2. * G4UniformRand() - 1. <<
708 smear.setY(2. * G4UniformRand() - 1. <<
709 smear.setZ(2. * G4UniformRand() - 1. <<
710 } while(smear.mag2() > 1.0); <<
711 facetNormal = normal + (1. - polish) * <<
712 } while(momentum * facetNormal >= 0.0); <<
713 facetNormal = facetNormal.unit(); <<
714 } <<
715 else <<
716 { <<
717 facetNormal = normal; <<
718 } <<
719 } <<
720 return facetNormal; <<
721 } <<
722 239
723 //....oooOO0OOooo........oooOO0OOooo........oo << 240 if (valid) {
724 void G4OpBoundaryProcess::DielectricMetal() << 241 theGlobalNormal = -theGlobalNormal;
725 { <<
726 G4int n = 0; <<
727 G4double rand; <<
728 G4ThreeVector A_trans; <<
729 <<
730 do <<
731 { <<
732 ++n; <<
733 rand = G4UniformRand(); <<
734 if(rand > fReflectivity && n == 1) <<
735 { <<
736 if(rand > fReflectivity + fTransmittance <<
737 { <<
738 DoAbsorption(); <<
739 } <<
740 else <<
741 { <<
742 fStatus = Transmission; <<
743 fNewMomentum = fOldMomentum; <<
744 fNewPolarization = fOldPolarization; <<
745 } <<
746 break; <<
747 } <<
748 else <<
749 { <<
750 if(fRealRIndexMPV && fImagRIndexMPV) <<
751 { <<
752 if(n > 1) <<
753 { <<
754 CalculateReflectivity(); <<
755 if(!G4BooleanRand(fReflectivity)) <<
756 { <<
757 DoAbsorption(); <<
758 break; <<
759 } <<
760 } 242 }
761 } << 243 else
762 if(fModel == glisur || fFinish == polish <<
763 { <<
764 DoReflection(); <<
765 } <<
766 else <<
767 { <<
768 if(n == 1) <<
769 ChooseReflection(); <<
770 if(fStatus == LambertianReflection) <<
771 { 244 {
772 DoReflection(); << 245 G4ExceptionDescription ed;
>> 246 ed << " G4OpBoundaryProcess/PostStepDoIt(): "
>> 247 << " The Navigator reports that it returned an invalid normal"
>> 248 << G4endl;
>> 249 G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun01",
>> 250 EventMustBeAborted,ed,
>> 251 "Invalid Surface Normal - Geometry must return valid surface normal");
773 } 252 }
774 else if(fStatus == BackScattering) << 253
775 { << 254 if (OldMomentum * theGlobalNormal > 0.0) {
776 fNewMomentum = -fOldMomentum; << 255 #ifdef G4OPTICAL_DEBUG
777 fNewPolarization = -fOldPolarization << 256 G4ExceptionDescription ed;
>> 257 ed << " G4OpBoundaryProcess/PostStepDoIt(): "
>> 258 << " theGlobalNormal points in a wrong direction. "
>> 259 << G4endl;
>> 260 ed << " The momentum of the photon arriving at interface (oldMomentum)"
>> 261 << " must exit the volume cross in the step. " << G4endl;
>> 262 ed << " So it MUST have dot < 0 with the normal that Exits the new volume (globalNormal)." << G4endl;
>> 263 ed << " >> The dot product of oldMomentum and global Normal is " << OldMomentum*theGlobalNormal << G4endl;
>> 264 ed << " Old Momentum (during step) = " << OldMomentum << G4endl;
>> 265 ed << " Global Normal (Exiting New Vol) = " << theGlobalNormal << G4endl;
>> 266 ed << G4endl;
>> 267 G4Exception("G4OpBoundaryProcess::PostStepDoIt", "OpBoun02",
>> 268 EventMustBeAborted, // Or JustWarning to see if it happens repeatedbly on one ray
>> 269 ed,
>> 270 "Invalid Surface Normal - Geometry must return valid surface normal pointing in the right direction");
>> 271 #else
>> 272 theGlobalNormal = -theGlobalNormal;
>> 273 #endif
778 } 274 }
779 else <<
780 { <<
781 if(fStatus == LobeReflection) <<
782 { <<
783 if(!fRealRIndexMPV || !fImagRIndex <<
784 { <<
785 fFacetNormal = GetFacetNormal(fO <<
786 } <<
787 // else <<
788 // case of complex rindex needs t <<
789 } <<
790 fNewMomentum = <<
791 fOldMomentum - 2. * fOldMomentum * <<
792 275
793 if(f_iTE > 0 && f_iTM > 0) << 276 G4MaterialPropertiesTable* aMaterialPropertiesTable;
794 { << 277 G4MaterialPropertyVector* Rindex;
795 fNewPolarization = << 278
796 -fOldPolarization + << 279 aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
797 (2. * fOldPolarization * fFacetN << 280 if (aMaterialPropertiesTable) {
798 } << 281 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
799 else if(f_iTE > 0) << 282 }
800 { << 283 else {
801 A_trans = (fSint1 > 0.0) ? fOldMom << 284 theStatus = NoRINDEX;
802 : fOldPol << 285 if ( verboseLevel > 0) BoundaryProcessVerbose();
803 fNewPolarization = -A_trans; << 286 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
>> 287 aParticleChange.ProposeTrackStatus(fStopAndKill);
>> 288 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 289 }
>> 290
>> 291 if (Rindex) {
>> 292 Rindex1 = Rindex->Value(thePhotonMomentum);
>> 293 }
>> 294 else {
>> 295 theStatus = NoRINDEX;
>> 296 if ( verboseLevel > 0) BoundaryProcessVerbose();
>> 297 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
>> 298 aParticleChange.ProposeTrackStatus(fStopAndKill);
>> 299 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 300 }
>> 301
>> 302 theReflectivity = 1.;
>> 303 theEfficiency = 0.;
>> 304 theTransmittance = 0.;
>> 305
>> 306 theSurfaceRoughness = 0.;
>> 307
>> 308 theModel = glisur;
>> 309 theFinish = polished;
>> 310
>> 311 G4SurfaceType type = dielectric_dielectric;
>> 312
>> 313 Rindex = NULL;
>> 314 OpticalSurface = NULL;
>> 315
>> 316 G4LogicalSurface* Surface = NULL;
>> 317
>> 318 Surface = G4LogicalBorderSurface::GetSurface(thePrePV, thePostPV);
>> 319
>> 320 if (Surface == NULL){
>> 321 G4bool enteredDaughter= (thePostPV->GetMotherLogical() ==
>> 322 thePrePV ->GetLogicalVolume());
>> 323 if(enteredDaughter){
>> 324 Surface =
>> 325 G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
>> 326 if(Surface == NULL)
>> 327 Surface =
>> 328 G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
>> 329 }
>> 330 else {
>> 331 Surface =
>> 332 G4LogicalSkinSurface::GetSurface(thePrePV->GetLogicalVolume());
>> 333 if(Surface == NULL)
>> 334 Surface =
>> 335 G4LogicalSkinSurface::GetSurface(thePostPV->GetLogicalVolume());
>> 336 }
>> 337 }
>> 338
>> 339 if (Surface) OpticalSurface =
>> 340 dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
>> 341
>> 342 if (OpticalSurface) {
>> 343
>> 344 type = OpticalSurface->GetType();
>> 345 theModel = OpticalSurface->GetModel();
>> 346 theFinish = OpticalSurface->GetFinish();
>> 347
>> 348 aMaterialPropertiesTable = OpticalSurface->
>> 349 GetMaterialPropertiesTable();
>> 350
>> 351 if (aMaterialPropertiesTable) {
>> 352
>> 353 if (theFinish == polishedbackpainted ||
>> 354 theFinish == groundbackpainted ) {
>> 355 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
>> 356 if (Rindex) {
>> 357 Rindex2 = Rindex->Value(thePhotonMomentum);
>> 358 }
>> 359 else {
>> 360 theStatus = NoRINDEX;
>> 361 if ( verboseLevel > 0) BoundaryProcessVerbose();
>> 362 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
>> 363 aParticleChange.ProposeTrackStatus(fStopAndKill);
>> 364 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 365 }
>> 366 }
>> 367
>> 368 PropertyPointer =
>> 369 aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
>> 370 PropertyPointer1 =
>> 371 aMaterialPropertiesTable->GetProperty("REALRINDEX");
>> 372 PropertyPointer2 =
>> 373 aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
>> 374
>> 375 iTE = 1;
>> 376 iTM = 1;
>> 377
>> 378 if (PropertyPointer) {
>> 379
>> 380 theReflectivity =
>> 381 PropertyPointer->Value(thePhotonMomentum);
>> 382
>> 383 } else if (PropertyPointer1 && PropertyPointer2) {
>> 384
>> 385 CalculateReflectivity();
>> 386
>> 387 }
>> 388
>> 389 PropertyPointer =
>> 390 aMaterialPropertiesTable->GetProperty("EFFICIENCY");
>> 391 if (PropertyPointer) {
>> 392 theEfficiency =
>> 393 PropertyPointer->Value(thePhotonMomentum);
>> 394 }
>> 395
>> 396 PropertyPointer =
>> 397 aMaterialPropertiesTable->GetProperty("TRANSMITTANCE");
>> 398 if (PropertyPointer) {
>> 399 theTransmittance =
>> 400 PropertyPointer->Value(thePhotonMomentum);
>> 401 }
>> 402
>> 403 if (aMaterialPropertiesTable->
>> 404 ConstPropertyExists("SURFACEROUGHNESS"))
>> 405 theSurfaceRoughness = aMaterialPropertiesTable->
>> 406 GetConstProperty("SURFACEROUGHNESS");
>> 407
>> 408 if ( theModel == unified ) {
>> 409 PropertyPointer =
>> 410 aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
>> 411 if (PropertyPointer) {
>> 412 prob_sl =
>> 413 PropertyPointer->Value(thePhotonMomentum);
>> 414 } else {
>> 415 prob_sl = 0.0;
>> 416 }
>> 417
>> 418 PropertyPointer =
>> 419 aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
>> 420 if (PropertyPointer) {
>> 421 prob_ss =
>> 422 PropertyPointer->Value(thePhotonMomentum);
>> 423 } else {
>> 424 prob_ss = 0.0;
>> 425 }
>> 426
>> 427 PropertyPointer =
>> 428 aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
>> 429 if (PropertyPointer) {
>> 430 prob_bs =
>> 431 PropertyPointer->Value(thePhotonMomentum);
>> 432 } else {
>> 433 prob_bs = 0.0;
>> 434 }
>> 435 }
>> 436 }
>> 437 else if (theFinish == polishedbackpainted ||
>> 438 theFinish == groundbackpainted ) {
>> 439 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
>> 440 aParticleChange.ProposeTrackStatus(fStopAndKill);
>> 441 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 442 }
>> 443 }
>> 444
>> 445 if (type == dielectric_dielectric ) {
>> 446 if (theFinish == polished || theFinish == ground ) {
>> 447
>> 448 if (Material1 == Material2){
>> 449 theStatus = SameMaterial;
>> 450 if ( verboseLevel > 0) BoundaryProcessVerbose();
>> 451 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 452 }
>> 453 aMaterialPropertiesTable =
>> 454 Material2->GetMaterialPropertiesTable();
>> 455 if (aMaterialPropertiesTable)
>> 456 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
>> 457 if (Rindex) {
>> 458 Rindex2 = Rindex->Value(thePhotonMomentum);
>> 459 }
>> 460 else {
>> 461 theStatus = NoRINDEX;
>> 462 if ( verboseLevel > 0) BoundaryProcessVerbose();
>> 463 aParticleChange.ProposeLocalEnergyDeposit(thePhotonMomentum);
>> 464 aParticleChange.ProposeTrackStatus(fStopAndKill);
>> 465 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 466 }
>> 467 }
>> 468 }
>> 469
>> 470 if (type == dielectric_metal) {
>> 471
>> 472 DielectricMetal();
>> 473
>> 474 }
>> 475 else if (type == dielectric_LUT) {
>> 476
>> 477 DielectricLUT();
>> 478
>> 479 }
>> 480 else if (type == dielectric_dichroic) {
>> 481
>> 482 DielectricDichroic();
>> 483
>> 484 }
>> 485 else if (type == dielectric_dielectric) {
>> 486
>> 487 if ( theFinish == polishedbackpainted ||
>> 488 theFinish == groundbackpainted ) {
>> 489 DielectricDielectric();
804 } 490 }
805 else if(f_iTM > 0) << 491 else {
806 { << 492 G4double rand = G4UniformRand();
807 fNewPolarization = << 493 if ( rand > theReflectivity ) {
808 -fNewMomentum.cross(A_trans).uni << 494 if (rand > theReflectivity + theTransmittance) {
>> 495 DoAbsorption();
>> 496 } else {
>> 497 theStatus = Transmission;
>> 498 NewMomentum = OldMomentum;
>> 499 NewPolarization = OldPolarization;
>> 500 }
>> 501 }
>> 502 else {
>> 503 if ( theFinish == polishedfrontpainted ) {
>> 504 DoReflection();
>> 505 }
>> 506 else if ( theFinish == groundfrontpainted ) {
>> 507 theStatus = LambertianReflection;
>> 508 DoReflection();
>> 509 }
>> 510 else {
>> 511 DielectricDielectric();
>> 512 }
>> 513 }
809 } 514 }
810 } 515 }
811 } << 516 else {
812 fOldMomentum = fNewMomentum; << 517
813 fOldPolarization = fNewPolarization; << 518 G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
814 } << 519 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
815 // Loop checking, 13-Aug-2015, Peter Gumpl << 520
816 } while(fNewMomentum * fGlobalNormal < 0.0); << 521 }
>> 522
>> 523 NewMomentum = NewMomentum.unit();
>> 524 NewPolarization = NewPolarization.unit();
>> 525
>> 526 if ( verboseLevel > 0) {
>> 527 G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
>> 528 G4cout << " New Polarization: " << NewPolarization << G4endl;
>> 529 BoundaryProcessVerbose();
>> 530 }
>> 531
>> 532 aParticleChange.ProposeMomentumDirection(NewMomentum);
>> 533 aParticleChange.ProposePolarization(NewPolarization);
>> 534
>> 535 if ( theStatus == FresnelRefraction || theStatus == Transmission ) {
>> 536 G4MaterialPropertyVector* groupvel =
>> 537 Material2->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
>> 538 G4double finalVelocity = groupvel->Value(thePhotonMomentum);
>> 539 aParticleChange.ProposeVelocity(finalVelocity);
>> 540 }
>> 541
>> 542 if ( theStatus == Detection && fInvokeSD ) InvokeSD(pStep);
>> 543
>> 544 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
817 } 545 }
818 546
819 //....oooOO0OOooo........oooOO0OOooo........oo << 547 void G4OpBoundaryProcess::BoundaryProcessVerbose() const
820 void G4OpBoundaryProcess::DielectricLUT() <<
821 { 548 {
822 G4int thetaIndex, phiIndex; << 549 if ( theStatus == Undefined )
823 G4double angularDistVal, thetaRad, phiRad; << 550 G4cout << " *** Undefined *** " << G4endl;
824 G4ThreeVector perpVectorTheta, perpVectorPhi << 551 if ( theStatus == Transmission )
825 << 552 G4cout << " *** Transmission *** " << G4endl;
826 fStatus = G4OpBoundaryProcessStatus( << 553 if ( theStatus == FresnelRefraction )
827 G4int(fFinish) + (G4int(NoRINDEX) - G4int( << 554 G4cout << " *** FresnelRefraction *** " << G4endl;
828 << 555 if ( theStatus == FresnelReflection )
829 G4int thetaIndexMax = fOpticalSurface->GetTh << 556 G4cout << " *** FresnelReflection *** " << G4endl;
830 G4int phiIndexMax = fOpticalSurface->GetPh << 557 if ( theStatus == TotalInternalReflection )
831 << 558 G4cout << " *** TotalInternalReflection *** " << G4endl;
832 G4double rand; << 559 if ( theStatus == LambertianReflection )
833 << 560 G4cout << " *** LambertianReflection *** " << G4endl;
834 do << 561 if ( theStatus == LobeReflection )
835 { << 562 G4cout << " *** LobeReflection *** " << G4endl;
836 rand = G4UniformRand(); << 563 if ( theStatus == SpikeReflection )
837 if(rand > fReflectivity) << 564 G4cout << " *** SpikeReflection *** " << G4endl;
838 { << 565 if ( theStatus == BackScattering )
839 if(rand > fReflectivity + fTransmittance << 566 G4cout << " *** BackScattering *** " << G4endl;
840 { << 567 if ( theStatus == PolishedLumirrorAirReflection )
841 DoAbsorption(); << 568 G4cout << " *** PolishedLumirrorAirReflection *** " << G4endl;
842 } << 569 if ( theStatus == PolishedLumirrorGlueReflection )
843 else << 570 G4cout << " *** PolishedLumirrorGlueReflection *** " << G4endl;
844 { << 571 if ( theStatus == PolishedAirReflection )
845 fStatus = Transmission; << 572 G4cout << " *** PolishedAirReflection *** " << G4endl;
846 fNewMomentum = fOldMomentum; << 573 if ( theStatus == PolishedTeflonAirReflection )
847 fNewPolarization = fOldPolarization; << 574 G4cout << " *** PolishedTeflonAirReflection *** " << G4endl;
848 } << 575 if ( theStatus == PolishedTiOAirReflection )
849 break; << 576 G4cout << " *** PolishedTiOAirReflection *** " << G4endl;
850 } << 577 if ( theStatus == PolishedTyvekAirReflection )
851 else << 578 G4cout << " *** PolishedTyvekAirReflection *** " << G4endl;
852 { << 579 if ( theStatus == PolishedVM2000AirReflection )
853 // Calculate Angle between Normal and Ph << 580 G4cout << " *** PolishedVM2000AirReflection *** " << G4endl;
854 G4double anglePhotonToNormal = fOldMomen << 581 if ( theStatus == PolishedVM2000GlueReflection )
855 // Round to closest integer: LBNL model << 582 G4cout << " *** PolishedVM2000GlueReflection *** " << G4endl;
856 G4int angleIncident = (G4int)std::lrint( << 583 if ( theStatus == EtchedLumirrorAirReflection )
857 << 584 G4cout << " *** EtchedLumirrorAirReflection *** " << G4endl;
858 // Take random angles THETA and PHI, << 585 if ( theStatus == EtchedLumirrorGlueReflection )
859 // and see if below Probability - if not << 586 G4cout << " *** EtchedLumirrorGlueReflection *** " << G4endl;
860 do << 587 if ( theStatus == EtchedAirReflection )
861 { << 588 G4cout << " *** EtchedAirReflection *** " << G4endl;
862 thetaIndex = (G4int)G4RandFlat::shootI << 589 if ( theStatus == EtchedTeflonAirReflection )
863 phiIndex = (G4int)G4RandFlat::shootI << 590 G4cout << " *** EtchedTeflonAirReflection *** " << G4endl;
864 // Find probability with the new indec << 591 if ( theStatus == EtchedTiOAirReflection )
865 angularDistVal = fOpticalSurface->GetA << 592 G4cout << " *** EtchedTiOAirReflection *** " << G4endl;
866 angleIncident, thetaIndex, phiIndex) << 593 if ( theStatus == EtchedTyvekAirReflection )
867 // Loop checking, 13-Aug-2015, Peter G << 594 G4cout << " *** EtchedTyvekAirReflection *** " << G4endl;
868 } while(!G4BooleanRand(angularDistVal)); << 595 if ( theStatus == EtchedVM2000AirReflection )
869 << 596 G4cout << " *** EtchedVM2000AirReflection *** " << G4endl;
870 thetaRad = G4double(-90 + 4 * thetaIndex << 597 if ( theStatus == EtchedVM2000GlueReflection )
871 phiRad = G4double(-90 + 5 * phiIndex) << 598 G4cout << " *** EtchedVM2000GlueReflection *** " << G4endl;
872 // Rotate Photon Momentum in Theta, then << 599 if ( theStatus == GroundLumirrorAirReflection )
873 fNewMomentum = -fOldMomentum; << 600 G4cout << " *** GroundLumirrorAirReflection *** " << G4endl;
874 << 601 if ( theStatus == GroundLumirrorGlueReflection )
875 perpVectorTheta = fNewMomentum.cross(fGl << 602 G4cout << " *** GroundLumirrorGlueReflection *** " << G4endl;
876 if(perpVectorTheta.mag() < fCarTolerance << 603 if ( theStatus == GroundAirReflection )
877 { << 604 G4cout << " *** GroundAirReflection *** " << G4endl;
878 perpVectorTheta = fNewMomentum.orthogo << 605 if ( theStatus == GroundTeflonAirReflection )
879 } << 606 G4cout << " *** GroundTeflonAirReflection *** " << G4endl;
880 fNewMomentum = << 607 if ( theStatus == GroundTiOAirReflection )
881 fNewMomentum.rotate(anglePhotonToNorma << 608 G4cout << " *** GroundTiOAirReflection *** " << G4endl;
882 perpVectorPhi = perpVectorTheta.cross(fN << 609 if ( theStatus == GroundTyvekAirReflection )
883 fNewMomentum = fNewMomentum.rotate(-phi << 610 G4cout << " *** GroundTyvekAirReflection *** " << G4endl;
884 << 611 if ( theStatus == GroundVM2000AirReflection )
885 // Rotate Polarization too: << 612 G4cout << " *** GroundVM2000AirReflection *** " << G4endl;
886 fFacetNormal = (fNewMomentum - fOldM << 613 if ( theStatus == GroundVM2000GlueReflection )
887 fNewPolarization = -fOldPolarization + << 614 G4cout << " *** GroundVM2000GlueReflection *** " << G4endl;
888 (2. * fOldPolarizatio << 615 if ( theStatus == Absorption )
889 } << 616 G4cout << " *** Absorption *** " << G4endl;
890 // Loop checking, 13-Aug-2015, Peter Gumpl << 617 if ( theStatus == Detection )
891 } while(fNewMomentum * fGlobalNormal <= 0.0) << 618 G4cout << " *** Detection *** " << G4endl;
>> 619 if ( theStatus == NotAtBoundary )
>> 620 G4cout << " *** NotAtBoundary *** " << G4endl;
>> 621 if ( theStatus == SameMaterial )
>> 622 G4cout << " *** SameMaterial *** " << G4endl;
>> 623 if ( theStatus == StepTooSmall )
>> 624 G4cout << " *** StepTooSmall *** " << G4endl;
>> 625 if ( theStatus == NoRINDEX )
>> 626 G4cout << " *** NoRINDEX *** " << G4endl;
>> 627 if ( theStatus == Dichroic )
>> 628 G4cout << " *** Dichroic Transmission *** " << G4endl;
892 } 629 }
893 630
894 //....oooOO0OOooo........oooOO0OOooo........oo << 631 G4ThreeVector
895 void G4OpBoundaryProcess::DielectricLUTDAVIS() << 632 G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
>> 633 const G4ThreeVector& Normal ) const
896 { 634 {
897 G4int angindex, random, angleIncident; << 635 G4ThreeVector FacetNormal;
898 G4double reflectivityValue, elevation, azimu <<
899 G4double anglePhotonToNormal; <<
900 <<
901 G4int lutbin = fOpticalSurface->GetLUTbins( <<
902 G4double rand = G4UniformRand(); <<
903 <<
904 G4double sinEl; <<
905 G4ThreeVector u, vNorm, w; <<
906 <<
907 do <<
908 { <<
909 anglePhotonToNormal = fOldMomentum.angle(- <<
910 <<
911 // Davis model has 90 reflection bins: rou <<
912 // don't allow angleIncident to be 90 for <<
913 angleIncident = std::min( <<
914 static_cast<G4int>(std::floor(anglePhoto <<
915 reflectivityValue = fOpticalSurface->GetRe <<
916 <<
917 if(rand > reflectivityValue) <<
918 { <<
919 if(fEfficiency > 0.) <<
920 { <<
921 DoAbsorption(); <<
922 break; <<
923 } <<
924 else <<
925 { <<
926 fStatus = Transmission; <<
927 636
928 if(angleIncident <= 0.01) << 637 if (theModel == unified || theModel == LUT) {
929 { <<
930 fNewMomentum = fOldMomentum; <<
931 break; <<
932 } <<
933 638
934 do << 639 /* This function code alpha to a random value taken from the
935 { << 640 distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
936 random = (G4int)G4RandFlat::shootInt << 641 for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
937 angindex = << 642 is a gaussian distribution with mean 0 and standard deviation
938 (((random * 2) - 1)) + angleIncide << 643 sigma_alpha. */
939 << 644
940 azimuth = << 645 G4double alpha;
941 fOpticalSurface->GetAngularDistrib << 646
942 elevation = fOpticalSurface->GetAngu << 647 G4double sigma_alpha = 0.0;
943 } while(elevation == 0. && azimuth == << 648 if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
944 << 649
945 sinEl = std::sin(elevation); << 650 if (sigma_alpha == 0.0) return FacetNormal = Normal;
946 vNorm = (fGlobalNormal.cross(fOldMomen << 651
947 u = vNorm.cross(fGlobalNormal) * ( << 652 G4double f_max = std::min(1.0,4.*sigma_alpha);
948 vNorm *= (sinEl * std::sin(azimuth)); << 653
949 // fGlobalNormal shouldn't be modified << 654 G4double phi, SinAlpha, CosAlpha, SinPhi, CosPhi, unit_x, unit_y, unit_z;
950 w = (fGlobalNormal *= std:: << 655 G4ThreeVector tmpNormal;
951 fNewMomentum = u + vNorm + w; << 656
952 << 657 do {
953 // Rotate Polarization too: << 658 do {
954 fFacetNormal = (fNewMomentum - fOl << 659 alpha = G4RandGauss::shoot(0.0,sigma_alpha);
955 fNewPolarization = -fOldPolarization + << 660 // Loop checking, 13-Aug-2015, Peter Gumplinger
956 << 661 } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
957 } << 662
958 } << 663 phi = G4UniformRand()*twopi;
959 else << 664
960 { << 665 SinAlpha = std::sin(alpha);
961 fStatus = LobeReflection; << 666 CosAlpha = std::cos(alpha);
962 << 667 SinPhi = std::sin(phi);
963 if(angleIncident == 0) << 668 CosPhi = std::cos(phi);
964 { << 669
965 fNewMomentum = -fOldMomentum; << 670 unit_x = SinAlpha * CosPhi;
966 break; << 671 unit_y = SinAlpha * SinPhi;
967 } << 672 unit_z = CosAlpha;
968 << 673
969 do << 674 FacetNormal.setX(unit_x);
970 { << 675 FacetNormal.setY(unit_y);
971 random = (G4int)G4RandFlat::shootInt << 676 FacetNormal.setZ(unit_z);
972 angindex = (((random * 2) - 1)) + (ang << 677
973 << 678 tmpNormal = Normal;
974 azimuth = fOpticalSurface->GetAngularD << 679
975 elevation = fOpticalSurface->GetAngula << 680 FacetNormal.rotateUz(tmpNormal);
976 } while(elevation == 0. && azimuth == 0. << 681 // Loop checking, 13-Aug-2015, Peter Gumplinger
977 << 682 } while (Momentum * FacetNormal >= 0.0);
978 sinEl = std::sin(elevation); << 683 }
979 vNorm = (fGlobalNormal.cross(fOldMomentu << 684 else {
980 u = vNorm.cross(fGlobalNormal) * (si << 685
981 vNorm *= (sinEl * std::sin(azimuth)); << 686 G4double polish = 1.0;
982 // fGlobalNormal shouldn't be modified h << 687 if (OpticalSurface) polish = OpticalSurface->GetPolish();
983 w = (fGlobalNormal *= std::cos(elevation << 688
984 << 689 if (polish < 1.0) {
985 fNewMomentum = u + vNorm + w; << 690 do {
986 << 691 G4ThreeVector smear;
987 // Rotate Polarization too: (needs revis << 692 do {
988 fNewPolarization = fOldPolarization; << 693 smear.setX(2.*G4UniformRand()-1.0);
989 } << 694 smear.setY(2.*G4UniformRand()-1.0);
990 } while(fNewMomentum * fGlobalNormal <= 0.0) << 695 smear.setZ(2.*G4UniformRand()-1.0);
>> 696 // Loop checking, 13-Aug-2015, Peter Gumplinger
>> 697 } while (smear.mag()>1.0);
>> 698 smear = (1.-polish) * smear;
>> 699 FacetNormal = Normal + smear;
>> 700 // Loop checking, 13-Aug-2015, Peter Gumplinger
>> 701 } while (Momentum * FacetNormal >= 0.0);
>> 702 FacetNormal = FacetNormal.unit();
>> 703 }
>> 704 else {
>> 705 FacetNormal = Normal;
>> 706 }
>> 707 }
>> 708 return FacetNormal;
991 } 709 }
992 710
993 //....oooOO0OOooo........oooOO0OOooo........oo << 711 void G4OpBoundaryProcess::DielectricMetal()
994 void G4OpBoundaryProcess::DielectricDichroic() <<
995 { 712 {
996 // Calculate Angle between Normal and Photon << 713 G4int n = 0;
997 G4double anglePhotonToNormal = fOldMomentum. << 714 G4double rand, PdotN, EdotN;
>> 715 G4ThreeVector A_trans, A_paral;
>> 716
>> 717 do {
>> 718
>> 719 n++;
>> 720
>> 721 rand = G4UniformRand();
>> 722 if ( rand > theReflectivity && n == 1 ) {
>> 723 if (rand > theReflectivity + theTransmittance) {
>> 724 DoAbsorption();
>> 725 } else {
>> 726 theStatus = Transmission;
>> 727 NewMomentum = OldMomentum;
>> 728 NewPolarization = OldPolarization;
>> 729 }
>> 730 break;
>> 731 }
>> 732 else {
>> 733
>> 734 if (PropertyPointer1 && PropertyPointer2) {
>> 735 if ( n > 1 ) {
>> 736 CalculateReflectivity();
>> 737 if ( !G4BooleanRand(theReflectivity) ) {
>> 738 DoAbsorption();
>> 739 break;
>> 740 }
>> 741 }
>> 742 }
>> 743
>> 744 if ( theModel == glisur || theFinish == polished ) {
>> 745
>> 746 DoReflection();
>> 747
>> 748 } else {
>> 749
>> 750 if ( n == 1 ) ChooseReflection();
>> 751
>> 752 if ( theStatus == LambertianReflection ) {
>> 753 DoReflection();
>> 754 }
>> 755 else if ( theStatus == BackScattering ) {
>> 756 NewMomentum = -OldMomentum;
>> 757 NewPolarization = -OldPolarization;
>> 758 }
>> 759 else {
>> 760
>> 761 if(theStatus==LobeReflection){
>> 762 if ( PropertyPointer1 && PropertyPointer2 ){
>> 763 } else {
>> 764 theFacetNormal =
>> 765 GetFacetNormal(OldMomentum,theGlobalNormal);
>> 766 }
>> 767 }
>> 768
>> 769 PdotN = OldMomentum * theFacetNormal;
>> 770 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
>> 771 EdotN = OldPolarization * theFacetNormal;
>> 772
>> 773 if (sint1 > 0.0 ) {
>> 774 A_trans = OldMomentum.cross(theFacetNormal);
>> 775 A_trans = A_trans.unit();
>> 776 } else {
>> 777 A_trans = OldPolarization;
>> 778 }
>> 779 A_paral = NewMomentum.cross(A_trans);
>> 780 A_paral = A_paral.unit();
>> 781
>> 782 if(iTE>0&&iTM>0) {
>> 783 NewPolarization =
>> 784 -OldPolarization + (2.*EdotN)*theFacetNormal;
>> 785 } else if (iTE>0) {
>> 786 NewPolarization = -A_trans;
>> 787 } else if (iTM>0) {
>> 788 NewPolarization = -A_paral;
>> 789 }
998 790
999 // Round it to closest integer << 791 }
1000 G4double angleIncident = std::floor(180. / <<
1001 792
1002 if(!fDichroicVector) << 793 }
1003 { <<
1004 if(fOpticalSurface) <<
1005 fDichroicVector = fOpticalSurface->GetD <<
1006 } <<
1007 794
1008 if(fDichroicVector) << 795 OldMomentum = NewMomentum;
1009 { << 796 OldPolarization = NewPolarization;
1010 G4double wavelength = h_Planck * c_light << 797
1011 fTransmittance = fDichroicVector->Va << 798 }
1012 i <<
1013 perCent; <<
1014 // G4cout << "wavelength: " << std::flo <<
1015 // << "nm" << <<
1016 // G4cout << "Incident angle: " << angl <<
1017 // G4cout << "Transmittance: " <<
1018 // << std::floor(fTransmittance/ <<
1019 } <<
1020 else <<
1021 { <<
1022 G4ExceptionDescription ed; <<
1023 ed << " G4OpBoundaryProcess/DielectricDic <<
1024 << " The dichroic surface has no G4Phy <<
1025 G4Exception("G4OpBoundaryProcess::Dielect <<
1026 FatalException, ed, <<
1027 "A dichroic surface must have <<
1028 } <<
1029 799
1030 if(!G4BooleanRand(fTransmittance)) <<
1031 { // Not transmitted, so reflect <<
1032 if(fModel == glisur || fFinish == polishe <<
1033 { <<
1034 DoReflection(); <<
1035 } <<
1036 else <<
1037 { <<
1038 ChooseReflection(); <<
1039 if(fStatus == LambertianReflection) <<
1040 { <<
1041 DoReflection(); <<
1042 } <<
1043 else if(fStatus == BackScattering) <<
1044 { <<
1045 fNewMomentum = -fOldMomentum; <<
1046 fNewPolarization = -fOldPolarization; <<
1047 } <<
1048 else <<
1049 { <<
1050 G4double PdotN, EdotN; <<
1051 do <<
1052 { <<
1053 if(fStatus == LobeReflection) <<
1054 { <<
1055 fFacetNormal = GetFacetNormal(fOl <<
1056 } <<
1057 PdotN = fOldMomentum * fFace <<
1058 fNewMomentum = fOldMomentum - (2. * <<
1059 // Loop checking, 13-Aug-2015, Pete 800 // Loop checking, 13-Aug-2015, Peter Gumplinger
1060 } while(fNewMomentum * fGlobalNormal << 801 } while (NewMomentum * theGlobalNormal < 0.0);
>> 802 }
1061 803
1062 EdotN = fOldPolarization * << 804 void G4OpBoundaryProcess::DielectricLUT()
1063 fNewPolarization = -fOldPolarization << 805 {
1064 } << 806 G4int thetaIndex, phiIndex;
1065 } << 807 G4double AngularDistributionValue, thetaRad, phiRad, EdotN;
1066 } << 808 G4ThreeVector PerpendicularVectorTheta, PerpendicularVectorPhi;
1067 else << 809
1068 { << 810 theStatus = G4OpBoundaryProcessStatus(G4int(theFinish) +
1069 fStatus = Dichroic; << 811 (G4int(NoRINDEX)-G4int(groundbackpainted)));
1070 fNewMomentum = fOldMomentum; << 812
1071 fNewPolarization = fOldPolarization; << 813 G4int thetaIndexMax = OpticalSurface->GetThetaIndexMax();
1072 } << 814 G4int phiIndexMax = OpticalSurface->GetPhiIndexMax();
>> 815
>> 816 G4double rand;
>> 817
>> 818 do {
>> 819 rand = G4UniformRand();
>> 820 if ( rand > theReflectivity ) {
>> 821 if (rand > theReflectivity + theTransmittance) {
>> 822 DoAbsorption();
>> 823 } else {
>> 824 theStatus = Transmission;
>> 825 NewMomentum = OldMomentum;
>> 826 NewPolarization = OldPolarization;
>> 827 }
>> 828 break;
>> 829 }
>> 830 else {
>> 831 // Calculate Angle between Normal and Photon Momentum
>> 832 G4double anglePhotonToNormal =
>> 833 OldMomentum.angle(-theGlobalNormal);
>> 834 // Round it to closest integer
>> 835 G4int angleIncident = G4int(std::floor(180/pi*anglePhotonToNormal+0.5));
>> 836
>> 837 // Take random angles THETA and PHI,
>> 838 // and see if below Probability - if not - Redo
>> 839 do {
>> 840 thetaIndex = G4RandFlat::shootInt(thetaIndexMax-1);
>> 841 phiIndex = G4RandFlat::shootInt(phiIndexMax-1);
>> 842 // Find probability with the new indeces from LUT
>> 843 AngularDistributionValue = OpticalSurface ->
>> 844 GetAngularDistributionValue(angleIncident,
>> 845 thetaIndex,
>> 846 phiIndex);
>> 847 // Loop checking, 13-Aug-2015, Peter Gumplinger
>> 848 } while ( !G4BooleanRand(AngularDistributionValue) );
>> 849
>> 850 thetaRad = (-90 + 4*thetaIndex)*pi/180;
>> 851 phiRad = (-90 + 5*phiIndex)*pi/180;
>> 852 // Rotate Photon Momentum in Theta, then in Phi
>> 853 NewMomentum = -OldMomentum;
>> 854
>> 855 PerpendicularVectorTheta = NewMomentum.cross(theGlobalNormal);
>> 856 if (PerpendicularVectorTheta.mag() < kCarTolerance )
>> 857 PerpendicularVectorTheta = NewMomentum.orthogonal();
>> 858 NewMomentum =
>> 859 NewMomentum.rotate(anglePhotonToNormal-thetaRad,
>> 860 PerpendicularVectorTheta);
>> 861 PerpendicularVectorPhi =
>> 862 PerpendicularVectorTheta.cross(NewMomentum);
>> 863 NewMomentum = NewMomentum.rotate(-phiRad,PerpendicularVectorPhi);
>> 864
>> 865 // Rotate Polarization too:
>> 866 theFacetNormal = (NewMomentum - OldMomentum).unit();
>> 867 EdotN = OldPolarization * theFacetNormal;
>> 868 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
>> 869 }
>> 870 // Loop checking, 13-Aug-2015, Peter Gumplinger
>> 871 } while (NewMomentum * theGlobalNormal <= 0.0);
1073 } 872 }
1074 873
1075 //....oooOO0OOooo........oooOO0OOooo........o << 874 void G4OpBoundaryProcess::DielectricDichroic()
1076 void G4OpBoundaryProcess::DielectricDielectri <<
1077 { 875 {
1078 G4bool inside = false; << 876 // Calculate Angle between Normal and Photon Momentum
1079 G4bool swap = false; << 877 G4double anglePhotonToNormal = OldMomentum.angle(-theGlobalNormal);
1080 878
1081 if(fFinish == polished) << 879 // Round it to closest integer
1082 { << 880 G4double angleIncident = std::floor(180/pi*anglePhotonToNormal+0.5);
1083 fFacetNormal = fGlobalNormal; <<
1084 } <<
1085 else <<
1086 { <<
1087 fFacetNormal = GetFacetNormal(fOldMomentu <<
1088 } <<
1089 G4double cost1 = -fOldMomentum * fFacetNorm <<
1090 G4double cost2 = 0.; <<
1091 G4double sint2 = 0.; <<
1092 <<
1093 G4bool surfaceRoughnessCriterionPass = true <<
1094 if(fSurfaceRoughness != 0. && fRindex1 > fR <<
1095 { <<
1096 G4double wavelength = h_Pl <<
1097 G4double surfaceRoughnessCriterion = std: <<
1098 (4. * pi * fSurfaceRoughness * fRindex1 <<
1099 surfaceRoughnessCriterionPass = G4Boolean <<
1100 } <<
1101 881
1102 leap: << 882 if (!DichroicVector) {
>> 883 if (OpticalSurface) DichroicVector = OpticalSurface->GetDichroicVector();
>> 884 }
1103 885
1104 G4bool through = false; <<
1105 G4bool done = false; <<
1106 886
1107 G4ThreeVector A_trans, A_paral, E1pp, E1pl; << 887 if (DichroicVector) {
1108 G4double E1_perp, E1_parl; << 888 G4double wavelength = h_Planck*c_light/thePhotonMomentum;
1109 G4double s1, s2, E2_perp, E2_parl, E2_total << 889 theTransmittance =
1110 G4double E2_abs, C_parl, C_perp; << 890 DichroicVector->Value(wavelength/nm,angleIncident,idx,idy)*perCent;
1111 G4double alpha; << 891 // G4cout << "wavelength: " << std::floor(wavelength/nm)
1112 << 892 // << "nm" << G4endl;
1113 do << 893 // G4cout << "Incident angle: " << angleIncident << "deg" << G4endl;
1114 { << 894 // G4cout << "Transmittance: "
1115 if(through) << 895 // << std::floor(theTransmittance/perCent) << "%" << G4endl;
1116 { << 896 } else {
1117 swap = !swap; << 897 G4ExceptionDescription ed;
1118 through = false; << 898 ed << " G4OpBoundaryProcess/DielectricDichroic(): "
1119 fGlobalNormal = -fGlobalNormal; << 899 << " The dichroic surface has no G4Physics2DVector"
1120 G4SwapPtr(fMaterial1, fMaterial2); << 900 << G4endl;
1121 G4SwapObj(&fRindex1, &fRindex2); << 901 G4Exception("G4OpBoundaryProcess::DielectricDichroic", "OpBoun03",
1122 } << 902 FatalException,ed,
1123 << 903 "A dichroic surface must have an associated G4Physics2DVector");
1124 if(fFinish == polished) << 904 }
1125 { << 905
1126 fFacetNormal = fGlobalNormal; << 906 if ( !G4BooleanRand(theTransmittance) ) { // Not transmitted, so reflect
1127 } << 907
1128 else << 908 if ( theModel == glisur || theFinish == polished ) {
1129 { << 909 DoReflection();
1130 fFacetNormal = GetFacetNormal(fOldMomen << 910 } else {
1131 } << 911 ChooseReflection();
1132 << 912 if ( theStatus == LambertianReflection ) {
1133 cost1 = -fOldMomentum * fFacetNormal; << 913 DoReflection();
1134 if(std::abs(cost1) < 1.0 - fCarTolerance) << 914 } else if ( theStatus == BackScattering ) {
1135 { << 915 NewMomentum = -OldMomentum;
1136 fSint1 = std::sqrt(1. - cost1 * cost1); << 916 NewPolarization = -OldPolarization;
1137 sint2 = fSint1 * fRindex1 / fRindex2; << 917 } else {
1138 // this isn't a sine as we might expect << 918 G4double PdotN, EdotN;
1139 } << 919 do {
1140 else << 920 if (theStatus==LobeReflection)
1141 { << 921 theFacetNormal = GetFacetNormal(OldMomentum,theGlobalNormal);
1142 fSint1 = 0.0; << 922 PdotN = OldMomentum * theFacetNormal;
1143 sint2 = 0.0; << 923 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
1144 } << 924 // Loop checking, 13-Aug-2015, Peter Gumplinger
1145 << 925 } while (NewMomentum * theGlobalNormal <= 0.0);
1146 // TOTAL INTERNAL REFLECTION << 926 EdotN = OldPolarization * theFacetNormal;
1147 if(sint2 >= 1.0) << 927 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
1148 { << 928 }
1149 swap = false; << 929 }
1150 << 930
1151 fStatus = TotalInternalReflection; << 931 } else {
1152 if(!surfaceRoughnessCriterionPass) << 932
1153 fStatus = LambertianReflection; << 933 theStatus = Dichroic;
1154 if(fModel == unified && fFinish != poli << 934 NewMomentum = OldMomentum;
1155 ChooseReflection(); << 935 NewPolarization = OldPolarization;
1156 if(fStatus == LambertianReflection) <<
1157 { <<
1158 DoReflection(); <<
1159 } <<
1160 else if(fStatus == BackScattering) <<
1161 { <<
1162 fNewMomentum = -fOldMomentum; <<
1163 fNewPolarization = -fOldPolarization; <<
1164 } <<
1165 else <<
1166 { <<
1167 fNewMomentum = <<
1168 fOldMomentum - 2. * fOldMomentum * <<
1169 fNewPolarization = -fOldPolarization <<
1170 <<
1171 } <<
1172 } <<
1173 // NOT TIR <<
1174 else if(sint2 < 1.0) <<
1175 { <<
1176 // Calculate amplitude for transmission <<
1177 if(cost1 > 0.0) <<
1178 { <<
1179 cost2 = std::sqrt(1. - sint2 * sint2) <<
1180 } <<
1181 else <<
1182 { <<
1183 cost2 = -std::sqrt(1. - sint2 * sint2 <<
1184 } <<
1185 <<
1186 if(fSint1 > 0.0) <<
1187 { <<
1188 A_trans = (fOldMomentum.cross(fFacetN <<
1189 E1_perp = fOldPolarization * A_trans; <<
1190 E1pp = E1_perp * A_trans; <<
1191 E1pl = fOldPolarization - E1pp; <<
1192 E1_parl = E1pl.mag(); <<
1193 } <<
1194 else <<
1195 { <<
1196 A_trans = fOldPolarization; <<
1197 // Here we Follow Jackson's conventio <<
1198 // component = 1 in case of a ray per <<
1199 E1_perp = 0.0; <<
1200 E1_parl = 1.0; <<
1201 } <<
1202 <<
1203 s1 = fRindex1 * cost1; <<
1204 E2_perp = 2. * s1 * E1_perp / (fRindex <<
1205 E2_parl = 2. * s1 * E1_parl / (fRindex <<
1206 E2_total = E2_perp * E2_perp + E2_parl <<
1207 s2 = fRindex2 * cost2 * E2_total; <<
1208 <<
1209 // D.Sawkey, 24 May 24 <<
1210 // Transmittance has already been taken <<
1211 // For e.g. specular surfaces, the rati <<
1212 // reflection should be given by the ma <<
1213 // TRANSMITTANCE <<
1214 //if(fTransmittance > 0.) <<
1215 // transCoeff = fTransmittance; <<
1216 //else if(cost1 != 0.0) <<
1217 if(cost1 != 0.0) <<
1218 transCoeff = s2 / s1; <<
1219 else <<
1220 transCoeff = 0.0; <<
1221 <<
1222 // NOT TIR: REFLECTION <<
1223 if(!G4BooleanRand(transCoeff)) <<
1224 { <<
1225 swap = false; <<
1226 fStatus = FresnelReflection; <<
1227 <<
1228 if(!surfaceRoughnessCriterionPass) <<
1229 fStatus = LambertianReflection; <<
1230 if(fModel == unified && fFinish != po <<
1231 ChooseReflection(); <<
1232 if(fStatus == LambertianReflection) <<
1233 { <<
1234 DoReflection(); <<
1235 } <<
1236 else if(fStatus == BackScattering) <<
1237 { <<
1238 fNewMomentum = -fOldMomentum; <<
1239 fNewPolarization = -fOldPolarizatio <<
1240 } <<
1241 else <<
1242 { <<
1243 fNewMomentum = <<
1244 fOldMomentum - 2. * fOldMomentum <<
1245 if(fSint1 > 0.0) <<
1246 { // incident ray oblique <<
1247 E2_parl = fRindex2 * E2_parl / f <<
1248 E2_perp = E2_perp - E1_perp; <<
1249 E2_total = E2_perp * E2_perp + E2 <<
1250 A_paral = (fNewMomentum.cross(A_ <<
1251 E2_abs = std::sqrt(E2_total); <<
1252 C_parl = E2_parl / E2_abs; <<
1253 C_perp = E2_perp / E2_abs; <<
1254 936
1255 fNewPolarization = C_parl * A_par <<
1256 } <<
1257 else <<
1258 { // incident ray perpendicular <<
1259 if(fRindex2 > fRindex1) <<
1260 { <<
1261 fNewPolarization = -fOldPolariz <<
1262 } <<
1263 else <<
1264 { <<
1265 fNewPolarization = fOldPolariza <<
1266 } <<
1267 } <<
1268 } 937 }
1269 } << 938 }
1270 // NOT TIR: TRANSMISSION <<
1271 else <<
1272 { <<
1273 inside = !inside; <<
1274 through = true; <<
1275 fStatus = FresnelRefraction; <<
1276 <<
1277 if(fSint1 > 0.0) <<
1278 { // incident ray oblique <<
1279 alpha = cost1 - cost2 * (fRi <<
1280 fNewMomentum = (fOldMomentum + alph <<
1281 A_paral = (fNewMomentum.cross( <<
1282 E2_abs = std::sqrt(E2_total); <<
1283 C_parl = E2_parl / E2_abs; <<
1284 C_perp = E2_perp / E2_abs; <<
1285 939
1286 fNewPolarization = C_parl * A_paral << 940 void G4OpBoundaryProcess::DielectricDielectric()
1287 } << 941 {
1288 else << 942 G4bool Inside = false;
1289 { // incident ray perpendicular << 943 G4bool Swap = false;
1290 fNewMomentum = fOldMomentum; << 944
1291 fNewPolarization = fOldPolarization << 945 G4bool SurfaceRoughnessCriterionPass = 1;
>> 946 if (theSurfaceRoughness != 0. && Rindex1 > Rindex2) {
>> 947 G4double wavelength = h_Planck*c_light/thePhotonMomentum;
>> 948 G4double SurfaceRoughnessCriterion =
>> 949 std::exp(-std::pow((4*pi*theSurfaceRoughness*Rindex1*cost1/wavelength),2));
>> 950 SurfaceRoughnessCriterionPass =
>> 951 G4BooleanRand(SurfaceRoughnessCriterion);
1292 } 952 }
1293 } <<
1294 } <<
1295 953
1296 fOldMomentum = fNewMomentum.unit(); << 954 leap:
1297 fOldPolarization = fNewPolarization.unit( <<
1298 955
1299 if(fStatus == FresnelRefraction) << 956 G4bool Through = false;
1300 { << 957 G4bool Done = false;
1301 done = (fNewMomentum * fGlobalNormal <= <<
1302 } <<
1303 else <<
1304 { <<
1305 done = (fNewMomentum * fGlobalNormal >= <<
1306 } <<
1307 // Loop checking, 13-Aug-2015, Peter Gump <<
1308 } while(!done); <<
1309 958
1310 if(inside && !swap) << 959 G4double PdotN, EdotN;
1311 { <<
1312 if(fFinish == polishedbackpainted || fFin <<
1313 { <<
1314 G4double rand = G4UniformRand(); <<
1315 if(rand > fReflectivity + fTransmittanc <<
1316 { <<
1317 DoAbsorption(); <<
1318 } <<
1319 else if(rand > fReflectivity) <<
1320 { <<
1321 fStatus = Transmission; <<
1322 fNewMomentum = fOldMomentum; <<
1323 fNewPolarization = fOldPolarization; <<
1324 } <<
1325 else <<
1326 { <<
1327 if(fStatus != FresnelRefraction) <<
1328 { <<
1329 fGlobalNormal = -fGlobalNormal; <<
1330 } <<
1331 else <<
1332 { <<
1333 swap = !swap; <<
1334 G4SwapPtr(fMaterial1, fMaterial2); <<
1335 G4SwapObj(&fRindex1, &fRindex2); <<
1336 } <<
1337 if(fFinish == groundbackpainted) <<
1338 fStatus = LambertianReflection; <<
1339 960
1340 DoReflection(); << 961 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
>> 962 G4double E1_perp, E1_parl;
>> 963 G4double s1, s2, E2_perp, E2_parl, E2_total, TransCoeff;
>> 964 G4double E2_abs, C_parl, C_perp;
>> 965 G4double alpha;
>> 966
>> 967 do {
>> 968
>> 969 if (Through) {
>> 970 Swap = !Swap;
>> 971 Through = false;
>> 972 theGlobalNormal = -theGlobalNormal;
>> 973 G4SwapPtr(Material1,Material2);
>> 974 G4SwapObj(&Rindex1,&Rindex2);
>> 975 }
>> 976
>> 977 if ( theFinish == polished ) {
>> 978 theFacetNormal = theGlobalNormal;
>> 979 }
>> 980 else {
>> 981 theFacetNormal =
>> 982 GetFacetNormal(OldMomentum,theGlobalNormal);
>> 983 }
>> 984
>> 985 PdotN = OldMomentum * theFacetNormal;
>> 986 EdotN = OldPolarization * theFacetNormal;
>> 987
>> 988 cost1 = - PdotN;
>> 989 if (std::abs(cost1) < 1.0-kCarTolerance){
>> 990 sint1 = std::sqrt(1.-cost1*cost1);
>> 991 sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
>> 992 }
>> 993 else {
>> 994 sint1 = 0.0;
>> 995 sint2 = 0.0;
>> 996 }
>> 997
>> 998 if (sint2 >= 1.0) {
>> 999
>> 1000 // Simulate total internal reflection
>> 1001
>> 1002 if (Swap) Swap = !Swap;
>> 1003
>> 1004 theStatus = TotalInternalReflection;
>> 1005
>> 1006 if ( !SurfaceRoughnessCriterionPass ) theStatus =
>> 1007 LambertianReflection;
>> 1008
>> 1009 if ( theModel == unified && theFinish != polished )
>> 1010 ChooseReflection();
>> 1011
>> 1012 if ( theStatus == LambertianReflection ) {
>> 1013 DoReflection();
>> 1014 }
>> 1015 else if ( theStatus == BackScattering ) {
>> 1016 NewMomentum = -OldMomentum;
>> 1017 NewPolarization = -OldPolarization;
>> 1018 }
>> 1019 else {
>> 1020
>> 1021 PdotN = OldMomentum * theFacetNormal;
>> 1022 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
>> 1023 EdotN = OldPolarization * theFacetNormal;
>> 1024 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
>> 1025
>> 1026 }
>> 1027 }
>> 1028 else if (sint2 < 1.0) {
>> 1029
>> 1030 // Calculate amplitude for transmission (Q = P x N)
>> 1031
>> 1032 if (cost1 > 0.0) {
>> 1033 cost2 = std::sqrt(1.-sint2*sint2);
>> 1034 }
>> 1035 else {
>> 1036 cost2 = -std::sqrt(1.-sint2*sint2);
>> 1037 }
>> 1038
>> 1039 if (sint1 > 0.0) {
>> 1040 A_trans = OldMomentum.cross(theFacetNormal);
>> 1041 A_trans = A_trans.unit();
>> 1042 E1_perp = OldPolarization * A_trans;
>> 1043 E1pp = E1_perp * A_trans;
>> 1044 E1pl = OldPolarization - E1pp;
>> 1045 E1_parl = E1pl.mag();
>> 1046 }
>> 1047 else {
>> 1048 A_trans = OldPolarization;
>> 1049 // Here we Follow Jackson's conventions and we set the
>> 1050 // parallel component = 1 in case of a ray perpendicular
>> 1051 // to the surface
>> 1052 E1_perp = 0.0;
>> 1053 E1_parl = 1.0;
>> 1054 }
>> 1055
>> 1056 s1 = Rindex1*cost1;
>> 1057 E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
>> 1058 E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
>> 1059 E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
>> 1060 s2 = Rindex2*cost2*E2_total;
>> 1061
>> 1062 if (theTransmittance > 0) TransCoeff = theTransmittance;
>> 1063 else if (cost1 != 0.0) TransCoeff = s2/s1;
>> 1064 else TransCoeff = 0.0;
>> 1065
>> 1066 if ( !G4BooleanRand(TransCoeff) ) {
>> 1067
>> 1068 // Simulate reflection
>> 1069
>> 1070 if (Swap) Swap = !Swap;
>> 1071
>> 1072 theStatus = FresnelReflection;
>> 1073
>> 1074 if ( !SurfaceRoughnessCriterionPass ) theStatus =
>> 1075 LambertianReflection;
>> 1076
>> 1077 if ( theModel == unified && theFinish != polished )
>> 1078 ChooseReflection();
>> 1079
>> 1080 if ( theStatus == LambertianReflection ) {
>> 1081 DoReflection();
>> 1082 }
>> 1083 else if ( theStatus == BackScattering ) {
>> 1084 NewMomentum = -OldMomentum;
>> 1085 NewPolarization = -OldPolarization;
>> 1086 }
>> 1087 else {
>> 1088
>> 1089 PdotN = OldMomentum * theFacetNormal;
>> 1090 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
>> 1091
>> 1092 if (sint1 > 0.0) { // incident ray oblique
>> 1093
>> 1094 E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
>> 1095 E2_perp = E2_perp - E1_perp;
>> 1096 E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
>> 1097 A_paral = NewMomentum.cross(A_trans);
>> 1098 A_paral = A_paral.unit();
>> 1099 E2_abs = std::sqrt(E2_total);
>> 1100 C_parl = E2_parl/E2_abs;
>> 1101 C_perp = E2_perp/E2_abs;
>> 1102
>> 1103 NewPolarization = C_parl*A_paral + C_perp*A_trans;
>> 1104
>> 1105 }
>> 1106
>> 1107 else { // incident ray perpendicular
>> 1108
>> 1109 if (Rindex2 > Rindex1) {
>> 1110 NewPolarization = - OldPolarization;
>> 1111 }
>> 1112 else {
>> 1113 NewPolarization = OldPolarization;
>> 1114 }
>> 1115
>> 1116 }
>> 1117 }
>> 1118 }
>> 1119 else { // photon gets transmitted
>> 1120
>> 1121 // Simulate transmission/refraction
>> 1122
>> 1123 Inside = !Inside;
>> 1124 Through = true;
>> 1125 theStatus = FresnelRefraction;
>> 1126
>> 1127 if (sint1 > 0.0) { // incident ray oblique
>> 1128
>> 1129 alpha = cost1 - cost2*(Rindex2/Rindex1);
>> 1130 NewMomentum = OldMomentum + alpha*theFacetNormal;
>> 1131 NewMomentum = NewMomentum.unit();
>> 1132 // PdotN = -cost2;
>> 1133 A_paral = NewMomentum.cross(A_trans);
>> 1134 A_paral = A_paral.unit();
>> 1135 E2_abs = std::sqrt(E2_total);
>> 1136 C_parl = E2_parl/E2_abs;
>> 1137 C_perp = E2_perp/E2_abs;
>> 1138
>> 1139 NewPolarization = C_parl*A_paral + C_perp*A_trans;
>> 1140
>> 1141 }
>> 1142 else { // incident ray perpendicular
>> 1143
>> 1144 NewMomentum = OldMomentum;
>> 1145 NewPolarization = OldPolarization;
>> 1146
>> 1147 }
>> 1148 }
>> 1149 }
>> 1150
>> 1151 OldMomentum = NewMomentum.unit();
>> 1152 OldPolarization = NewPolarization.unit();
>> 1153
>> 1154 if (theStatus == FresnelRefraction) {
>> 1155 Done = (NewMomentum * theGlobalNormal <= 0.0);
>> 1156 }
>> 1157 else {
>> 1158 Done = (NewMomentum * theGlobalNormal >= -kCarTolerance);
>> 1159 }
1341 1160
1342 fGlobalNormal = -fGlobalNormal; << 1161 // Loop checking, 13-Aug-2015, Peter Gumplinger
1343 fOldMomentum = fNewMomentum; << 1162 } while (!Done);
1344 1163
1345 goto leap; << 1164 if (Inside && !Swap) {
1346 } << 1165 if( theFinish == polishedbackpainted ||
1347 } << 1166 theFinish == groundbackpainted ) {
1348 } << 1167
>> 1168 G4double rand = G4UniformRand();
>> 1169 if ( rand > theReflectivity ) {
>> 1170 if (rand > theReflectivity + theTransmittance) {
>> 1171 DoAbsorption();
>> 1172 } else {
>> 1173 theStatus = Transmission;
>> 1174 NewMomentum = OldMomentum;
>> 1175 NewPolarization = OldPolarization;
>> 1176 }
>> 1177 }
>> 1178 else {
>> 1179 if (theStatus != FresnelRefraction ) {
>> 1180 theGlobalNormal = -theGlobalNormal;
>> 1181 }
>> 1182 else {
>> 1183 Swap = !Swap;
>> 1184 G4SwapPtr(Material1,Material2);
>> 1185 G4SwapObj(&Rindex1,&Rindex2);
>> 1186 }
>> 1187 if ( theFinish == groundbackpainted )
>> 1188 theStatus = LambertianReflection;
>> 1189
>> 1190 DoReflection();
>> 1191
>> 1192 theGlobalNormal = -theGlobalNormal;
>> 1193 OldMomentum = NewMomentum;
>> 1194
>> 1195 goto leap;
>> 1196 }
>> 1197 }
>> 1198 }
1349 } 1199 }
1350 1200
1351 //....oooOO0OOooo........oooOO0OOooo........o << 1201 // GetMeanFreePath
1352 G4double G4OpBoundaryProcess::GetMeanFreePath << 1202 // ---------------
>> 1203 //
>> 1204 G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
>> 1205 G4double ,
1353 1206 G4ForceCondition* condition)
1354 { 1207 {
1355 *condition = Forced; << 1208 *condition = Forced;
1356 return DBL_MAX; << 1209
>> 1210 return DBL_MAX;
1357 } 1211 }
1358 1212
1359 //....oooOO0OOooo........oooOO0OOooo........o << 1213 G4double G4OpBoundaryProcess::GetIncidentAngle()
1360 G4double G4OpBoundaryProcess::GetIncidentAngl <<
1361 { 1214 {
1362 return pi - std::acos(fOldMomentum * fFacet << 1215 G4double PdotN = OldMomentum * theFacetNormal;
1363 (fOldMomentum.mag() * << 1216 G4double magP= OldMomentum.mag();
>> 1217 G4double magN= theFacetNormal.mag();
>> 1218 G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
>> 1219
>> 1220 return incidentangle;
1364 } 1221 }
1365 1222
1366 //....oooOO0OOooo........oooOO0OOooo........o <<
1367 G4double G4OpBoundaryProcess::GetReflectivity 1223 G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
1368 1224 G4double E1_parl,
1369 1225 G4double incidentangle,
1370 << 1226 G4double RealRindex,
1371 << 1227 G4double ImaginaryRindex)
1372 { 1228 {
1373 G4complex reflectivity, reflectivity_TE, re << 1229 G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
1374 G4complex N1(fRindex1, 0.), N2(realRindex, << 1230 G4complex N1(Rindex1, 0), N2(RealRindex, ImaginaryRindex);
1375 G4complex cosPhi; << 1231 G4complex CosPhi;
1376 1232
1377 G4complex u(1., 0.); // unit number 1 << 1233 G4complex u(1,0); //unit number 1
1378 1234
1379 G4complex numeratorTE; // E1_perp=1 E1_par << 1235 G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
1380 G4complex numeratorTM; // E1_parl=1 E1_per << 1236 G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
1381 G4complex denominatorTE, denominatorTM; 1237 G4complex denominatorTE, denominatorTM;
1382 G4complex rTM, rTE; 1238 G4complex rTM, rTE;
1383 1239
1384 G4MaterialPropertiesTable* MPT = fMaterial1 << 1240 G4MaterialPropertiesTable* aMaterialPropertiesTable =
1385 G4MaterialPropertyVector* ppR = MPT->GetPr << 1241 Material1->GetMaterialPropertiesTable();
1386 G4MaterialPropertyVector* ppI = MPT->GetPr << 1242 G4MaterialPropertyVector* aPropertyPointerR =
1387 if(ppR && ppI) << 1243 aMaterialPropertiesTable->GetProperty("REALRINDEX");
1388 { << 1244 G4MaterialPropertyVector* aPropertyPointerI =
1389 G4double rRindex = ppR->Value(fPhotonMome << 1245 aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
1390 G4double iRindex = ppI->Value(fPhotonMome << 1246 if (aPropertyPointerR && aPropertyPointerI) {
1391 N1 = G4complex(rRindex, iRi << 1247 G4double RRindex = aPropertyPointerR->Value(thePhotonMomentum);
>> 1248 G4double IRindex = aPropertyPointerI->Value(thePhotonMomentum);
>> 1249 N1 = G4complex(RRindex,IRindex);
1392 } 1250 }
1393 1251
1394 // Following two equations, rTM and rTE, ar 1252 // Following two equations, rTM and rTE, are from: "Introduction To Modern
1395 // Optics" written by Fowles 1253 // Optics" written by Fowles
1396 cosPhi = std::sqrt(u - ((std::sin(incidenta <<
1397 (N1 * N1) / (N2 * N <<
1398 1254
1399 numeratorTE = N1 * std::cos(incidentangle << 1255 CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))*(N1*N1)/(N2*N2)));
1400 denominatorTE = N1 * std::cos(incidentangle <<
1401 rTE = numeratorTE / denominatorTE <<
1402 <<
1403 numeratorTM = N2 * std::cos(incidentangle <<
1404 denominatorTM = N2 * std::cos(incidentangle <<
1405 rTM = numeratorTM / denominatorTM <<
1406 1256
1407 // This is my (PG) calculaton for reflectiv << 1257 numeratorTE = N1*std::cos(incidentangle) - N2*CosPhi;
>> 1258 denominatorTE = N1*std::cos(incidentangle) + N2*CosPhi;
>> 1259 rTE = numeratorTE/denominatorTE;
>> 1260
>> 1261 numeratorTM = N2*std::cos(incidentangle) - N1*CosPhi;
>> 1262 denominatorTM = N2*std::cos(incidentangle) + N1*CosPhi;
>> 1263 rTM = numeratorTM/denominatorTM;
>> 1264
>> 1265 // This is my calculaton for reflectivity on a metalic surface
1408 // depending on the fraction of TE and TM p 1266 // depending on the fraction of TE and TM polarization
1409 // when TE polarization, E1_parl=0 and E1_p 1267 // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
1410 // when TM polarization, E1_parl=1 and E1_p 1268 // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
1411 1269
1412 reflectivity_TE = (rTE * conj(rTE)) * (E1_p << 1270 Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
1413 (E1_perp * E1_perp + E1_p << 1271 / (E1_perp*E1_perp + E1_parl*E1_parl);
1414 reflectivity_TM = (rTM * conj(rTM)) * (E1_p << 1272 Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
1415 (E1_perp * E1_perp + E1_p << 1273 / (E1_perp*E1_perp + E1_parl*E1_parl);
1416 reflectivity = reflectivity_TE + reflectivi << 1274 Reflectivity = Reflectivity_TE + Reflectivity_TM;
1417 << 1275
1418 do << 1276 do {
1419 { << 1277 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))
1420 if(G4UniformRand() * real(reflectivity) > << 1278 {iTE = -1;}else{iTE = 1;}
1421 { << 1279 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))
1422 f_iTE = -1; << 1280 {iTM = -1;}else{iTM = 1;}
1423 } <<
1424 else <<
1425 { <<
1426 f_iTE = 1; <<
1427 } <<
1428 if(G4UniformRand() * real(reflectivity) > <<
1429 { <<
1430 f_iTM = -1; <<
1431 } <<
1432 else <<
1433 { <<
1434 f_iTM = 1; <<
1435 } <<
1436 // Loop checking, 13-Aug-2015, Peter Gump 1281 // Loop checking, 13-Aug-2015, Peter Gumplinger
1437 } while(f_iTE < 0 && f_iTM < 0); << 1282 } while(iTE<0&&iTM<0);
>> 1283
>> 1284 return real(Reflectivity);
1438 1285
1439 return real(reflectivity); <<
1440 } 1286 }
1441 1287
1442 //....oooOO0OOooo........oooOO0OOooo........o <<
1443 void G4OpBoundaryProcess::CalculateReflectivi 1288 void G4OpBoundaryProcess::CalculateReflectivity()
1444 { 1289 {
1445 G4double realRindex = fRealRIndexMPV->Value << 1290 G4double RealRindex =
1446 G4double imaginaryRindex = << 1291 PropertyPointer1->Value(thePhotonMomentum);
1447 fImagRIndexMPV->Value(fPhotonMomentum, id << 1292 G4double ImaginaryRindex =
>> 1293 PropertyPointer2->Value(thePhotonMomentum);
1448 1294
1449 // calculate FacetNormal 1295 // calculate FacetNormal
1450 if(fFinish == ground) << 1296 if ( theFinish == ground ) {
1451 { << 1297 theFacetNormal =
1452 fFacetNormal = GetFacetNormal(fOldMomentu << 1298 GetFacetNormal(OldMomentum, theGlobalNormal);
1453 } << 1299 } else {
1454 else << 1300 theFacetNormal = theGlobalNormal;
1455 { <<
1456 fFacetNormal = fGlobalNormal; <<
1457 } 1301 }
1458 1302
1459 G4double cost1 = -fOldMomentum * fFacetNorm << 1303 G4double PdotN = OldMomentum * theFacetNormal;
1460 if(std::abs(cost1) < 1.0 - fCarTolerance) << 1304 cost1 = -PdotN;
1461 { << 1305
1462 fSint1 = std::sqrt(1. - cost1 * cost1); << 1306 if (std::abs(cost1) < 1.0 - kCarTolerance) {
1463 } << 1307 sint1 = std::sqrt(1. - cost1*cost1);
1464 else << 1308 } else {
1465 { << 1309 sint1 = 0.0;
1466 fSint1 = 0.0; <<
1467 } 1310 }
1468 1311
1469 G4ThreeVector A_trans, A_paral, E1pp, E1pl; 1312 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
1470 G4double E1_perp, E1_parl; 1313 G4double E1_perp, E1_parl;
1471 1314
1472 if(fSint1 > 0.0) << 1315 if (sint1 > 0.0 ) {
1473 { << 1316 A_trans = OldMomentum.cross(theFacetNormal);
1474 A_trans = (fOldMomentum.cross(fFacetNorma << 1317 A_trans = A_trans.unit();
1475 E1_perp = fOldPolarization * A_trans; << 1318 E1_perp = OldPolarization * A_trans;
1476 E1pp = E1_perp * A_trans; << 1319 E1pp = E1_perp * A_trans;
1477 E1pl = fOldPolarization - E1pp; << 1320 E1pl = OldPolarization - E1pp;
1478 E1_parl = E1pl.mag(); << 1321 E1_parl = E1pl.mag();
1479 } << 1322 }
1480 else << 1323 else {
1481 { << 1324 A_trans = OldPolarization;
1482 A_trans = fOldPolarization; << 1325 // Here we Follow Jackson's conventions and we set the
1483 // Here we Follow Jackson's conventions a << 1326 // parallel component = 1 in case of a ray perpendicular
1484 // component = 1 in case of a ray perpend << 1327 // to the surface
1485 E1_perp = 0.0; << 1328 E1_perp = 0.0;
1486 E1_parl = 1.0; << 1329 E1_parl = 1.0;
1487 } 1330 }
1488 1331
>> 1332 //calculate incident angle
1489 G4double incidentangle = GetIncidentAngle() 1333 G4double incidentangle = GetIncidentAngle();
1490 1334
1491 // calculate the reflectivity depending on << 1335 //calculate the reflectivity depending on incident angle,
1492 // polarization and complex refractive << 1336 //polarization and complex refractive
1493 fReflectivity = GetReflectivity(E1_perp, E1 << 1337
1494 imaginaryRi << 1338 theReflectivity =
>> 1339 GetReflectivity(E1_perp, E1_parl, incidentangle,
>> 1340 RealRindex, ImaginaryRindex);
1495 } 1341 }
1496 1342
1497 //....oooOO0OOooo........oooOO0OOooo........o <<
1498 G4bool G4OpBoundaryProcess::InvokeSD(const G4 1343 G4bool G4OpBoundaryProcess::InvokeSD(const G4Step* pStep)
1499 { 1344 {
1500 G4Step aStep = *pStep; 1345 G4Step aStep = *pStep;
1501 aStep.AddTotalEnergyDeposit(fPhotonMomentum <<
1502 <<
1503 G4VSensitiveDetector* sd = aStep.GetPostSte <<
1504 if(sd != nullptr) <<
1505 return sd->Hit(&aStep); <<
1506 else <<
1507 return false; <<
1508 } <<
1509 <<
1510 //....oooOO0OOooo........oooOO0OOooo........o <<
1511 inline void G4OpBoundaryProcess::SetInvokeSD( <<
1512 { <<
1513 fInvokeSD = flag; <<
1514 G4OpticalParameters::Instance()->SetBoundar <<
1515 } <<
1516 <<
1517 //....oooOO0OOooo........oooOO0OOooo........o <<
1518 void G4OpBoundaryProcess::SetVerboseLevel(G4i <<
1519 { <<
1520 verboseLevel = verbose; <<
1521 G4OpticalParameters::Instance()->SetBoundar <<
1522 } <<
1523 <<
1524 //....oooOO0OOooo........oooOO0OOooo........o <<
1525 void G4OpBoundaryProcess::CoatedDielectricDie <<
1526 { <<
1527 G4MaterialPropertyVector* pp = nullptr; <<
1528 <<
1529 G4MaterialPropertiesTable* MPT = fMaterial2 <<
1530 if((pp = MPT->GetProperty(kRINDEX))) <<
1531 { <<
1532 fRindex2 = pp->Value(fPhotonMomentum, idx <<
1533 } <<
1534 <<
1535 MPT = fOpticalSurface->GetMaterialPropertie <<
1536 if((pp = MPT->GetProperty(kCOATEDRINDEX))) <<
1537 { <<
1538 fCoatedRindex = pp->Value(fPhotonMomentum <<
1539 } <<
1540 if(MPT->ConstPropertyExists(kCOATEDTHICKNES <<
1541 { <<
1542 fCoatedThickness = MPT->GetConstProperty( <<
1543 } <<
1544 if(MPT->ConstPropertyExists(kCOATEDFRUSTRAT <<
1545 { <<
1546 fCoatedFrustratedTransmission = <<
1547 (G4bool)MPT->GetConstProperty(kCOATEDFR <<
1548 } <<
1549 <<
1550 G4double sintTL; <<
1551 G4double wavelength = h_Planck * c_light / <<
1552 G4double PdotN; <<
1553 G4double E1_perp, E1_parl; <<
1554 G4double s1, E2_perp, E2_parl, E2_total, tr <<
1555 G4double E2_abs, C_parl, C_perp; <<
1556 G4double alpha; <<
1557 G4ThreeVector A_trans, A_paral, E1pp, E1pl; <<
1558 //G4bool Inside = false; <<
1559 //G4bool Swap = false; <<
1560 G4bool through = false; <<
1561 G4bool done = false; <<
1562 <<
1563 do { <<
1564 if (through) <<
1565 { <<
1566 //Swap = !Swap; <<
1567 through = false; <<
1568 fGlobalNormal = -fGlobalNormal; <<
1569 G4SwapPtr(fMaterial1, fMaterial2); <<
1570 G4SwapObj(&fRindex1, &fRindex2); <<
1571 } <<
1572 <<
1573 if(fFinish == polished) <<
1574 { <<
1575 fFacetNormal = fGlobalNormal; <<
1576 } <<
1577 else <<
1578 { <<
1579 fFacetNormal = GetFacetNormal(fOldMomen <<
1580 } <<
1581 <<
1582 PdotN = fOldMomentum * fFacetNormal; <<
1583 G4double cost1 = -PdotN; <<
1584 G4double sint2, cost2 = 0.; <<
1585 <<
1586 if (std::abs(cost1) < 1.0 - fCarTolerance <<
1587 { <<
1588 fSint1 = std::sqrt(1. - cost1 * cost1); <<
1589 sint2 = fSint1 * fRindex1 / fRindex2; <<
1590 sintTL = fSint1 * fRindex1 / fCoatedRin <<
1591 } else <<
1592 { <<
1593 fSint1 = 0.0; <<
1594 sint2 = 0.0; <<
1595 sintTL = 0.0; <<
1596 } <<
1597 <<
1598 if (fSint1 > 0.0) <<
1599 { <<
1600 A_trans = fOldMomentum.cross(fFacetNorm <<
1601 A_trans = A_trans.unit(); <<
1602 E1_perp = fOldPolarization * A_trans; <<
1603 E1pp = E1_perp * A_trans; <<
1604 E1pl = fOldPolarization - E1pp; <<
1605 E1_parl = E1pl.mag(); <<
1606 } <<
1607 else <<
1608 { <<
1609 A_trans = fOldPolarization; <<
1610 E1_perp = 0.0; <<
1611 E1_parl = 1.0; <<
1612 } <<
1613 <<
1614 s1 = fRindex1 * cost1; <<
1615 <<
1616 if (cost1 > 0.0) <<
1617 { <<
1618 cost2 = std::sqrt(1. - sint2 * sint2); <<
1619 } <<
1620 else <<
1621 { <<
1622 cost2 = -std::sqrt(1. - sint2 * sint2); <<
1623 } <<
1624 <<
1625 transCoeff = 0.0; <<
1626 <<
1627 if (sintTL >= 1.0) <<
1628 { // --> Angle > Angle Limit <<
1629 //Swap = false; <<
1630 } <<
1631 E2_perp = 2. * s1 * E1_perp / (fRindex1 * <<
1632 E2_parl = 2. * s1 * E1_parl / (fRindex2 * <<
1633 E2_total = E2_perp * E2_perp + E2_parl * <<
1634 <<
1635 transCoeff = 1. - GetReflectivityThroughT <<
1636 sintTL, E1_perp, E1_p <<
1637 if (!G4BooleanRand(transCoeff)) <<
1638 { <<
1639 if(verboseLevel > 2) <<
1640 G4cout << "Reflection from " << fMate <<
1641 << fMaterial2->GetName() << G4 <<
1642 <<
1643 //Swap = false; <<
1644 <<
1645 if (sintTL >= 1.0) <<
1646 { <<
1647 fStatus = TotalInternalReflection; <<
1648 } <<
1649 else <<
1650 { <<
1651 fStatus = CoatedDielectricReflection; <<
1652 } <<
1653 <<
1654 PdotN = fOldMomentum * fFacetNormal; <<
1655 fNewMomentum = fOldMomentum - (2. * Pdo <<
1656 <<
1657 if (fSint1 > 0.0) { // incident ray o <<
1658 <<
1659 E2_parl = fRindex2 * E2_parl / fRinde <<
1660 E2_perp = E2_perp - E1_perp; <<
1661 E2_total = E2_perp * E2_perp + E2_par <<
1662 A_paral = fNewMomentum.cross(A_trans) <<
1663 A_paral = A_paral.unit(); <<
1664 E2_abs = std::sqrt(E2_total); <<
1665 C_parl = E2_parl / E2_abs; <<
1666 C_perp = E2_perp / E2_abs; <<
1667 <<
1668 fNewPolarization = C_parl * A_paral + <<
1669 <<
1670 } <<
1671 else <<
1672 { // incident ray perpend <<
1673 if (fRindex2 > fRindex1) <<
1674 { <<
1675 fNewPolarization = -fOldPolarizatio <<
1676 } <<
1677 else <<
1678 { <<
1679 fNewPolarization = fOldPolarization <<
1680 } <<
1681 } <<
1682 <<
1683 } else { // photon gets transmitted <<
1684 if (verboseLevel > 2) <<
1685 G4cout << "Transmission from " << fMa <<
1686 << fMaterial2->GetName() << G4 <<
1687 <<
1688 //Inside = !Inside; <<
1689 through = true; <<
1690 <<
1691 if (fEfficiency > 0.) <<
1692 { <<
1693 DoAbsorption(); <<
1694 return; <<
1695 } <<
1696 else <<
1697 { <<
1698 if (sintTL >= 1.0) <<
1699 { <<
1700 fStatus = CoatedDielectricFrustrate <<
1701 } <<
1702 else <<
1703 { <<
1704 fStatus = CoatedDielectricRefractio <<
1705 } <<
1706 <<
1707 if (fSint1 > 0.0) { // incident <<
1708 <<
1709 alpha = cost1 - cost2 * (fRindex2 / <<
1710 fNewMomentum = fOldMomentum + alpha <<
1711 fNewMomentum = fNewMomentum.unit(); <<
1712 A_paral = fNewMomentum.cross(A_tran <<
1713 A_paral = A_paral.unit(); <<
1714 E2_abs = std::sqrt(E2_total); <<
1715 C_parl = E2_parl / E2_abs; <<
1716 C_perp = E2_perp / E2_abs; <<
1717 <<
1718 fNewPolarization = C_parl * A_paral <<
1719 <<
1720 } <<
1721 else <<
1722 { // incident ray pe <<
1723 fNewMomentum = fOldMomentum; <<
1724 fNewPolarization = fOldPolarization <<
1725 } <<
1726 } <<
1727 } <<
1728 1346
1729 fOldMomentum = fNewMomentum.unit(); << 1347 aStep.AddTotalEnergyDeposit(thePhotonMomentum);
1730 fOldPolarization = fNewPolarization.unit( <<
1731 if ((fStatus == CoatedDielectricFrustrate <<
1732 (fStatus == CoatedDielectricRefractio <<
1733 { <<
1734 done = (fNewMomentum * fGlobalNormal <= <<
1735 } <<
1736 else <<
1737 { <<
1738 done = (fNewMomentum * fGlobalNormal >= <<
1739 } <<
1740 1348
1741 } while (!done); << 1349 G4VSensitiveDetector* sd = aStep.GetPostStepPoint()->GetSensitiveDetector();
1742 } << 1350 if (sd) return sd->Hit(&aStep);
1743 << 1351 else return false;
1744 //....oooOO0OOooo........oooOO0OOooo........o <<
1745 G4double G4OpBoundaryProcess::GetReflectivity <<
1746 G4double E1_perp, <<
1747 G4double E1_parl, <<
1748 G4double wavelength, G4dou <<
1749 G4complex Reflectivity, Reflectivity_TE, Re <<
1750 G4double gammaTL, costTL; <<
1751 <<
1752 G4complex i(0, 1); <<
1753 G4complex rTM, rTE; <<
1754 G4complex r1toTL, rTLto2; <<
1755 G4double k0 = 2 * pi / wavelength; <<
1756 <<
1757 // Angle > Angle limit <<
1758 if (sinTL >= 1.0) { <<
1759 if (fCoatedFrustratedTransmission) { //Fr <<
1760 <<
1761 if (cost1 > 0.0) <<
1762 { <<
1763 gammaTL = std::sqrt(fRindex1 * fRinde <<
1764 fCoatedRindex * fCoatedRin <<
1765 } <<
1766 else <<
1767 { <<
1768 gammaTL = -std::sqrt(fRindex1 * fRind <<
1769 fCoatedRindex * fCoatedRin <<
1770 } <<
1771 <<
1772 // TE <<
1773 r1toTL = (fRindex1 * cost1 - i * gammaT <<
1774 rTLto2 = (i * gammaTL - fRindex2 * cost <<
1775 if (cost1 != 0.0) <<
1776 { <<
1777 rTE = (r1toTL + rTLto2 * std::exp(-2 <<
1778 (1.0 + r1toTL * rTLto2 * std <<
1779 } <<
1780 // TM <<
1781 r1toTL = (fRindex1 * i * gammaTL - fCoa <<
1782 (fRindex1 * i * gammaTL + f <<
1783 rTLto2 = (fCoatedRindex * fCoatedRindex <<
1784 (fCoatedRindex * fCoatedRin <<
1785 if (cost1 != 0.0) <<
1786 { <<
1787 rTM = (r1toTL + rTLto2 * std::exp(-2 <<
1788 (1.0 + r1toTL * rTLto2 * std <<
1789 } <<
1790 } <<
1791 else <<
1792 { //Total reflection <<
1793 return(1.); <<
1794 } <<
1795 } <<
1796 <<
1797 // Angle <= Angle limit <<
1798 else //if (sinTL < 1.0) <<
1799 { <<
1800 if (cost1 > 0.0) <<
1801 { <<
1802 costTL = std::sqrt(1. - sinTL * sinTL); <<
1803 } <<
1804 else <<
1805 { <<
1806 costTL = -std::sqrt(1. - sinTL * sinTL) <<
1807 } <<
1808 // TE <<
1809 r1toTL = (fRindex1 * cost1 - fCoatedRinde <<
1810 rTLto2 = (fCoatedRindex * costTL - fRinde <<
1811 if (cost1 != 0.0) <<
1812 { <<
1813 rTE = (r1toTL + rTLto2 * std::exp(2.0 * <<
1814 (1.0 + r1toTL * rTLto2 * std::exp <<
1815 } <<
1816 // TM <<
1817 r1toTL = (fRindex1 * costTL - fCoatedRind <<
1818 rTLto2 = (fCoatedRindex * cost2 - fRindex <<
1819 if (cost1 != 0.0) <<
1820 { <<
1821 rTM = (r1toTL + rTLto2 * std::exp(2.0 * <<
1822 (1.0 + r1toTL * rTLto2 * std::exp <<
1823 } <<
1824 } <<
1825 <<
1826 Reflectivity_TE = (rTE * conj(rTE)) * (E1_p <<
1827 Reflectivity_TM = (rTM * conj(rTM)) * (E1_p <<
1828 Reflectivity = Reflectivity_TE + Reflectivi <<
1829 <<
1830 return real(Reflectivity); <<
1831 } 1352 }
1832 1353