source: trunk/source/processes/optical/src/G4OpBoundaryProcess.cc@ 831

Last change on this file since 831 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 30.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26////////////////////////////////////////////////////////////////////////
27// Optical Photon Boundary Process Class Implementation
28////////////////////////////////////////////////////////////////////////
29//
30// File: G4OpBoundaryProcess.cc
31// Description: Discrete Process -- reflection/refraction at
32// optical interfaces
33// Version: 1.1
34// Created: 1997-06-18
35// Modified: 1998-05-25 - Correct parallel component of polarization
36// (thanks to: Stefano Magni + Giovanni Pieri)
37// 1998-05-28 - NULL Rindex pointer before reuse
38// (thanks to: Stefano Magni)
39// 1998-06-11 - delete *sint1 in oblique reflection
40// (thanks to: Giovanni Pieri)
41// 1998-06-19 - move from GetLocalExitNormal() to the new
42// method: GetLocalExitNormal(&valid) to get
43// the surface normal in all cases
44// 1998-11-07 - NULL OpticalSurface pointer before use
45// comparison not sharp for: std::abs(cost1) < 1.0
46// remove sin1, sin2 in lines 556,567
47// (thanks to Stefano Magni)
48// 1999-10-10 - Accommodate changes done in DoAbsorption by
49// changing logic in DielectricMetal
50// 2001-10-18 - avoid Linux (gcc-2.95.2) warning about variables
51// might be used uninitialized in this function
52// moved E2_perp, E2_parl and E2_total out of 'if'
53// 2003-11-27 - Modified line 168-9 to reflect changes made to
54// G4OpticalSurface class ( by Fan Lei)
55// 2004-02-02 - Set theStatus = Undefined at start of DoIt
56// 2005-07-28 - add G4ProcessType to constructor
57// 2006-11-04 - add capability of calculating the reflectivity
58// off a metal surface by way of a complex index
59// of refraction - Thanks to Sehwook Lee and John
60// Hauptman (Dept. of Physics - Iowa State Univ.)
61//
62// Author: Peter Gumplinger
63// adopted from work by Werner Keil - April 2/96
64// mail: gum@triumf.ca
65//
66////////////////////////////////////////////////////////////////////////
67
68#include "G4ios.hh"
69#include "G4OpBoundaryProcess.hh"
70#include "G4GeometryTolerance.hh"
71
72/////////////////////////
73// Class Implementation
74/////////////////////////
75
76 //////////////
77 // Operators
78 //////////////
79
80// G4OpBoundaryProcess::operator=(const G4OpBoundaryProcess &right)
81// {
82// }
83
84 /////////////////
85 // Constructors
86 /////////////////
87
88G4OpBoundaryProcess::G4OpBoundaryProcess(const G4String& processName,
89 G4ProcessType type)
90 : G4VDiscreteProcess(processName, type)
91{
92 if ( verboseLevel > 0) {
93 G4cout << GetProcessName() << " is created " << G4endl;
94 }
95
96 theStatus = Undefined;
97 theModel = glisur;
98 theFinish = polished;
99 theReflectivity = 1.;
100 theEfficiency = 0.;
101
102 prob_sl = 0.;
103 prob_ss = 0.;
104 prob_bs = 0.;
105
106 kCarTolerance = G4GeometryTolerance::GetInstance()
107 ->GetSurfaceTolerance();
108
109}
110
111// G4OpBoundaryProcess::G4OpBoundaryProcess(const G4OpBoundaryProcess &right)
112// {
113// }
114
115 ////////////////
116 // Destructors
117 ////////////////
118
119G4OpBoundaryProcess::~G4OpBoundaryProcess(){}
120
121 ////////////
122 // Methods
123 ////////////
124
125// PostStepDoIt
126// ------------
127//
128G4VParticleChange*
129G4OpBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
130{
131 theStatus = Undefined;
132
133 aParticleChange.Initialize(aTrack);
134
135 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
136 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
137
138 if (pPostStepPoint->GetStepStatus() != fGeomBoundary){
139 theStatus = NotAtBoundary;
140 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
141 }
142 if (aTrack.GetStepLength()<=kCarTolerance/2){
143 theStatus = StepTooSmall;
144 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
145 }
146
147 Material1 = pPreStepPoint -> GetMaterial();
148 Material2 = pPostStepPoint -> GetMaterial();
149
150 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
151
152 thePhotonMomentum = aParticle->GetTotalMomentum();
153 OldMomentum = aParticle->GetMomentumDirection();
154 OldPolarization = aParticle->GetPolarization();
155
156 G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition();
157
158 G4Navigator* theNavigator =
159 G4TransportationManager::GetTransportationManager()->
160 GetNavigatorForTracking();
161
162 G4ThreeVector theLocalPoint = theNavigator->
163 GetGlobalToLocalTransform().
164 TransformPoint(theGlobalPoint);
165
166 G4ThreeVector theLocalNormal; // Normal points back into volume
167
168 G4bool valid;
169 theLocalNormal = theNavigator->GetLocalExitNormal(&valid);
170
171 if (valid) {
172 theLocalNormal = -theLocalNormal;
173 }
174 else {
175 G4cerr << " G4OpBoundaryProcess/PostStepDoIt(): "
176 << " The Navigator reports that it returned an invalid normal"
177 << G4endl;
178 }
179
180 theGlobalNormal = theNavigator->GetLocalToGlobalTransform().
181 TransformAxis(theLocalNormal);
182
183 if (OldMomentum * theGlobalNormal > 0.0) {
184#ifdef G4DEBUG_OPTICAL
185 G4cerr << " G4OpBoundaryProcess/PostStepDoIt(): "
186 << " theGlobalNormal points the wrong direction "
187 << G4endl;
188#endif
189 theGlobalNormal = -theGlobalNormal;
190 }
191
192 G4MaterialPropertiesTable* aMaterialPropertiesTable;
193 G4MaterialPropertyVector* Rindex;
194
195 aMaterialPropertiesTable = Material1->GetMaterialPropertiesTable();
196 if (aMaterialPropertiesTable) {
197 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
198 }
199 else {
200 theStatus = NoRINDEX;
201 aParticleChange.ProposeTrackStatus(fStopAndKill);
202 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
203 }
204
205 if (Rindex) {
206 Rindex1 = Rindex->GetProperty(thePhotonMomentum);
207 }
208 else {
209 theStatus = NoRINDEX;
210 aParticleChange.ProposeTrackStatus(fStopAndKill);
211 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
212 }
213
214 theModel = glisur;
215 theFinish = polished;
216
217 G4SurfaceType type = dielectric_dielectric;
218
219 Rindex = NULL;
220 OpticalSurface = NULL;
221
222 G4LogicalSurface* Surface = G4LogicalBorderSurface::GetSurface
223 (pPreStepPoint ->GetPhysicalVolume(),
224 pPostStepPoint->GetPhysicalVolume());
225
226 if (Surface == NULL){
227 G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume()
228 ->GetMotherLogical() ==
229 pPreStepPoint->GetPhysicalVolume()
230 ->GetLogicalVolume());
231 if(enteredDaughter){
232 Surface = G4LogicalSkinSurface::GetSurface
233 (pPostStepPoint->GetPhysicalVolume()->
234 GetLogicalVolume());
235 if(Surface == NULL)
236 Surface = G4LogicalSkinSurface::GetSurface
237 (pPreStepPoint->GetPhysicalVolume()->
238 GetLogicalVolume());
239 }
240 else {
241 Surface = G4LogicalSkinSurface::GetSurface
242 (pPreStepPoint->GetPhysicalVolume()->
243 GetLogicalVolume());
244 if(Surface == NULL)
245 Surface = G4LogicalSkinSurface::GetSurface
246 (pPostStepPoint->GetPhysicalVolume()->
247 GetLogicalVolume());
248 }
249 }
250
251 // if (Surface) OpticalSurface = dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
252 if (Surface) OpticalSurface = (G4OpticalSurface*) Surface->GetSurfaceProperty();
253
254 if (OpticalSurface) {
255
256 type = OpticalSurface->GetType();
257 theModel = OpticalSurface->GetModel();
258 theFinish = OpticalSurface->GetFinish();
259
260 aMaterialPropertiesTable = OpticalSurface->
261 GetMaterialPropertiesTable();
262
263 if (aMaterialPropertiesTable) {
264
265 if (theFinish == polishedbackpainted ||
266 theFinish == groundbackpainted ) {
267 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
268 if (Rindex) {
269 Rindex2 = Rindex->GetProperty(thePhotonMomentum);
270 }
271 else {
272 theStatus = NoRINDEX;
273 aParticleChange.ProposeTrackStatus(fStopAndKill);
274 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
275 }
276 }
277
278 G4MaterialPropertyVector* PropertyPointer;
279 G4MaterialPropertyVector* PropertyPointer1;
280 G4MaterialPropertyVector* PropertyPointer2;
281
282 PropertyPointer =
283 aMaterialPropertiesTable->GetProperty("REFLECTIVITY");
284 PropertyPointer1 =
285 aMaterialPropertiesTable->GetProperty("REALRINDEX");
286 PropertyPointer2 =
287 aMaterialPropertiesTable->GetProperty("IMAGINARYRINDEX");
288
289 iTE = 1;
290 iTM = 1;
291
292 if (PropertyPointer) {
293
294 theReflectivity =
295 PropertyPointer->GetProperty(thePhotonMomentum);
296
297 } else if (PropertyPointer1 && PropertyPointer2) {
298
299 G4double RealRindex =
300 PropertyPointer1->GetProperty(thePhotonMomentum);
301 G4double ImaginaryRindex =
302 PropertyPointer2->GetProperty(thePhotonMomentum);
303
304 // calculate FacetNormal
305 if ( theFinish == ground ) {
306 theFacetNormal =
307 GetFacetNormal(OldMomentum, theGlobalNormal);
308 } else {
309 theFacetNormal = theGlobalNormal;
310 }
311
312 G4double PdotN = OldMomentum * theFacetNormal;
313 cost1 = -PdotN;
314
315 if (std::abs(cost1) < 1.0 - kCarTolerance) {
316 sint1 = std::sqrt(1. - cost1*cost1);
317 } else {
318 sint1 = 0.0;
319 }
320
321 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
322 G4double E1_perp, E1_parl;
323
324 if (sint1 > 0.0 ) {
325 A_trans = OldMomentum.cross(theFacetNormal);
326 A_trans = A_trans.unit();
327 E1_perp = OldPolarization * A_trans;
328 E1pp = E1_perp * A_trans;
329 E1pl = OldPolarization - E1pp;
330 E1_parl = E1pl.mag();
331 }
332 else {
333 A_trans = OldPolarization;
334 // Here we Follow Jackson's conventions and we set the
335 // parallel component = 1 in case of a ray perpendicular
336 // to the surface
337 E1_perp = 0.0;
338 E1_parl = 1.0;
339 }
340
341 //calculate incident angle
342 G4double incidentangle = GetIncidentAngle();
343
344 //calculate the reflectivity depending on incident angle,
345 //polarization and complex refractive
346
347 theReflectivity =
348 GetReflectivity(E1_perp, E1_parl, incidentangle,
349 RealRindex, ImaginaryRindex);
350
351 } else {
352 theReflectivity = 1.0;
353 }
354
355 PropertyPointer =
356 aMaterialPropertiesTable->GetProperty("EFFICIENCY");
357 if (PropertyPointer) {
358 theEfficiency =
359 PropertyPointer->GetProperty(thePhotonMomentum);
360 } else {
361 theEfficiency = 0.0;
362 }
363
364 if ( theModel == unified ) {
365 PropertyPointer =
366 aMaterialPropertiesTable->GetProperty("SPECULARLOBECONSTANT");
367 if (PropertyPointer) {
368 prob_sl =
369 PropertyPointer->GetProperty(thePhotonMomentum);
370 } else {
371 prob_sl = 0.0;
372 }
373
374 PropertyPointer =
375 aMaterialPropertiesTable->GetProperty("SPECULARSPIKECONSTANT");
376 if (PropertyPointer) {
377 prob_ss =
378 PropertyPointer->GetProperty(thePhotonMomentum);
379 } else {
380 prob_ss = 0.0;
381 }
382
383 PropertyPointer =
384 aMaterialPropertiesTable->GetProperty("BACKSCATTERCONSTANT");
385 if (PropertyPointer) {
386 prob_bs =
387 PropertyPointer->GetProperty(thePhotonMomentum);
388 } else {
389 prob_bs = 0.0;
390 }
391 }
392 }
393 else if (theFinish == polishedbackpainted ||
394 theFinish == groundbackpainted ) {
395 aParticleChange.ProposeTrackStatus(fStopAndKill);
396 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
397 }
398 }
399
400 if (type == dielectric_dielectric ) {
401 if (theFinish == polished || theFinish == ground ) {
402
403 if (Material1 == Material2){
404 theStatus = SameMaterial;
405 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
406 }
407 aMaterialPropertiesTable =
408 Material2->GetMaterialPropertiesTable();
409 if (aMaterialPropertiesTable)
410 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
411 if (Rindex) {
412 Rindex2 = Rindex->GetProperty(thePhotonMomentum);
413 }
414 else {
415 theStatus = NoRINDEX;
416 aParticleChange.ProposeTrackStatus(fStopAndKill);
417 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
418 }
419 }
420 }
421
422 if ( verboseLevel > 0 ) {
423 G4cout << " Photon at Boundary! " << G4endl;
424 G4cout << " Old Momentum Direction: " << OldMomentum << G4endl;
425 G4cout << " Old Polarization: " << OldPolarization << G4endl;
426 }
427
428 if (type == dielectric_metal) {
429
430 DielectricMetal();
431
432 }
433 else if (type == dielectric_dielectric) {
434
435 if ( theFinish == polishedfrontpainted ||
436 theFinish == groundfrontpainted ) {
437
438 if( !G4BooleanRand(theReflectivity) ) {
439 DoAbsorption();
440 }
441 else {
442 if ( theFinish == groundfrontpainted )
443 theStatus = LambertianReflection;
444 DoReflection();
445 }
446 }
447 else {
448 DielectricDielectric();
449 }
450 }
451 else {
452
453 G4cerr << " Error: G4BoundaryProcess: illegal boundary type " << G4endl;
454 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
455
456 }
457
458 NewMomentum = NewMomentum.unit();
459 NewPolarization = NewPolarization.unit();
460
461 if ( verboseLevel > 0) {
462 G4cout << " New Momentum Direction: " << NewMomentum << G4endl;
463 G4cout << " New Polarization: " << NewPolarization << G4endl;
464 if ( theStatus == Undefined )
465 G4cout << " *** Undefined *** " << G4endl;
466 if ( theStatus == FresnelRefraction )
467 G4cout << " *** FresnelRefraction *** " << G4endl;
468 if ( theStatus == FresnelReflection )
469 G4cout << " *** FresnelReflection *** " << G4endl;
470 if ( theStatus == TotalInternalReflection )
471 G4cout << " *** TotalInternalReflection *** " << G4endl;
472 if ( theStatus == LambertianReflection )
473 G4cout << " *** LambertianReflection *** " << G4endl;
474 if ( theStatus == LobeReflection )
475 G4cout << " *** LobeReflection *** " << G4endl;
476 if ( theStatus == SpikeReflection )
477 G4cout << " *** SpikeReflection *** " << G4endl;
478 if ( theStatus == BackScattering )
479 G4cout << " *** BackScattering *** " << G4endl;
480 if ( theStatus == Absorption )
481 G4cout << " *** Absorption *** " << G4endl;
482 if ( theStatus == Detection )
483 G4cout << " *** Detection *** " << G4endl;
484 if ( theStatus == NotAtBoundary )
485 G4cout << " *** NotAtBoundary *** " << G4endl;
486 if ( theStatus == SameMaterial )
487 G4cout << " *** SameMaterial *** " << G4endl;
488 if ( theStatus == StepTooSmall )
489 G4cout << " *** StepTooSmall *** " << G4endl;
490 if ( theStatus == NoRINDEX )
491 G4cout << " *** NoRINDEX *** " << G4endl;
492 }
493
494 aParticleChange.ProposeMomentumDirection(NewMomentum);
495 aParticleChange.ProposePolarization(NewPolarization);
496
497 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
498}
499
500G4ThreeVector
501G4OpBoundaryProcess::GetFacetNormal(const G4ThreeVector& Momentum,
502 const G4ThreeVector& Normal ) const
503{
504 G4ThreeVector FacetNormal;
505
506 if (theModel == unified) {
507
508 /* This function code alpha to a random value taken from the
509 distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
510 for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
511 is a gaussian distribution with mean 0 and standard deviation
512 sigma_alpha. */
513
514 G4double alpha;
515
516 G4double sigma_alpha = 0.0;
517 if (OpticalSurface) sigma_alpha = OpticalSurface->GetSigmaAlpha();
518
519 G4double f_max = std::min(1.0,4.*sigma_alpha);
520
521 do {
522 do {
523 alpha = G4RandGauss::shoot(0.0,sigma_alpha);
524 } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= halfpi );
525
526 G4double phi = G4UniformRand()*twopi;
527
528 G4double SinAlpha = std::sin(alpha);
529 G4double CosAlpha = std::cos(alpha);
530 G4double SinPhi = std::sin(phi);
531 G4double CosPhi = std::cos(phi);
532
533 G4double unit_x = SinAlpha * CosPhi;
534 G4double unit_y = SinAlpha * SinPhi;
535 G4double unit_z = CosAlpha;
536
537 FacetNormal.setX(unit_x);
538 FacetNormal.setY(unit_y);
539 FacetNormal.setZ(unit_z);
540
541 G4ThreeVector tmpNormal = Normal;
542
543 FacetNormal.rotateUz(tmpNormal);
544 } while (Momentum * FacetNormal >= 0.0);
545 }
546 else {
547
548 G4double polish = 1.0;
549 if (OpticalSurface) polish = OpticalSurface->GetPolish();
550
551 if (polish < 1.0) {
552 do {
553 G4ThreeVector smear;
554 do {
555 smear.setX(2.*G4UniformRand()-1.0);
556 smear.setY(2.*G4UniformRand()-1.0);
557 smear.setZ(2.*G4UniformRand()-1.0);
558 } while (smear.mag()>1.0);
559 smear = (1.-polish) * smear;
560 FacetNormal = Normal + smear;
561 } while (Momentum * FacetNormal >= 0.0);
562 FacetNormal = FacetNormal.unit();
563 }
564 else {
565 FacetNormal = Normal;
566 }
567 }
568 return FacetNormal;
569}
570
571void G4OpBoundaryProcess::DielectricMetal()
572{
573 G4int n = 0;
574
575 do {
576
577 n++;
578
579 if( !G4BooleanRand(theReflectivity) && n == 1 ) {
580
581 DoAbsorption();
582 break;
583
584 }
585 else {
586
587 if ( theModel == glisur || theFinish == polished ) {
588
589 DoReflection();
590
591 } else {
592
593 if ( n == 1 ) ChooseReflection();
594
595 if ( theStatus == LambertianReflection ) {
596 DoReflection();
597 }
598 else if ( theStatus == BackScattering ) {
599 NewMomentum = -OldMomentum;
600 NewPolarization = -OldPolarization;
601 }
602 else {
603
604 if(theStatus==LobeReflection)theFacetNormal =
605 GetFacetNormal(OldMomentum,theGlobalNormal);
606
607 G4double PdotN = OldMomentum * theFacetNormal;
608 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
609 G4double EdotN = OldPolarization * theFacetNormal;
610
611 G4ThreeVector A_trans, A_paral;
612
613 if (sint1 > 0.0 ) {
614 A_trans = OldMomentum.cross(theFacetNormal);
615 A_trans = A_trans.unit();
616 } else {
617 A_trans = OldPolarization;
618 }
619 A_paral = NewMomentum.cross(A_trans);
620 A_paral = A_paral.unit();
621
622 if(iTE>0&&iTM>0) {
623 NewPolarization =
624 -OldPolarization + (2.*EdotN)*theFacetNormal;
625 } else if (iTE>0) {
626 NewPolarization = -A_trans;
627 } else if (iTM>0) {
628 NewPolarization = -A_paral;
629 }
630
631 }
632
633 }
634
635 OldMomentum = NewMomentum;
636 OldPolarization = NewPolarization;
637
638 }
639
640 } while (NewMomentum * theGlobalNormal < 0.0);
641}
642
643void G4OpBoundaryProcess::DielectricDielectric()
644{
645 G4bool Inside = false;
646 G4bool Swap = false;
647
648 leap:
649
650 G4bool Through = false;
651 G4bool Done = false;
652
653 do {
654
655 if (Through) {
656 Swap = !Swap;
657 Through = false;
658 theGlobalNormal = -theGlobalNormal;
659 G4Swap(Material1,Material2);
660 G4Swap(&Rindex1,&Rindex2);
661 }
662
663 if ( theFinish == ground || theFinish == groundbackpainted ) {
664 theFacetNormal =
665 GetFacetNormal(OldMomentum,theGlobalNormal);
666 }
667 else {
668 theFacetNormal = theGlobalNormal;
669 }
670
671 G4double PdotN = OldMomentum * theFacetNormal;
672 G4double EdotN = OldPolarization * theFacetNormal;
673
674 cost1 = - PdotN;
675 if (std::abs(cost1) < 1.0-kCarTolerance){
676 sint1 = std::sqrt(1.-cost1*cost1);
677 sint2 = sint1*Rindex1/Rindex2; // *** Snell's Law ***
678 }
679 else {
680 sint1 = 0.0;
681 sint2 = 0.0;
682 }
683
684 if (sint2 >= 1.0) {
685
686 // Simulate total internal reflection
687
688 if (Swap) Swap = !Swap;
689
690 theStatus = TotalInternalReflection;
691
692 if ( theModel == unified && theFinish != polished )
693 ChooseReflection();
694
695 if ( theStatus == LambertianReflection ) {
696 DoReflection();
697 }
698 else if ( theStatus == BackScattering ) {
699 NewMomentum = -OldMomentum;
700 NewPolarization = -OldPolarization;
701 }
702 else {
703
704 PdotN = OldMomentum * theFacetNormal;
705 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
706 EdotN = OldPolarization * theFacetNormal;
707 NewPolarization = -OldPolarization + (2.*EdotN)*theFacetNormal;
708
709 }
710 }
711 else if (sint2 < 1.0) {
712
713 // Calculate amplitude for transmission (Q = P x N)
714
715 if (cost1 > 0.0) {
716 cost2 = std::sqrt(1.-sint2*sint2);
717 }
718 else {
719 cost2 = -std::sqrt(1.-sint2*sint2);
720 }
721
722 G4ThreeVector A_trans, A_paral, E1pp, E1pl;
723 G4double E1_perp, E1_parl;
724
725 if (sint1 > 0.0) {
726 A_trans = OldMomentum.cross(theFacetNormal);
727 A_trans = A_trans.unit();
728 E1_perp = OldPolarization * A_trans;
729 E1pp = E1_perp * A_trans;
730 E1pl = OldPolarization - E1pp;
731 E1_parl = E1pl.mag();
732 }
733 else {
734 A_trans = OldPolarization;
735 // Here we Follow Jackson's conventions and we set the
736 // parallel component = 1 in case of a ray perpendicular
737 // to the surface
738 E1_perp = 0.0;
739 E1_parl = 1.0;
740 }
741
742 G4double s1 = Rindex1*cost1;
743 G4double E2_perp = 2.*s1*E1_perp/(Rindex1*cost1+Rindex2*cost2);
744 G4double E2_parl = 2.*s1*E1_parl/(Rindex2*cost1+Rindex1*cost2);
745 G4double E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
746 G4double s2 = Rindex2*cost2*E2_total;
747
748 G4double TransCoeff;
749
750 if (cost1 != 0.0) {
751 TransCoeff = s2/s1;
752 }
753 else {
754 TransCoeff = 0.0;
755 }
756
757 G4double E2_abs, C_parl, C_perp;
758
759 if ( !G4BooleanRand(TransCoeff) ) {
760
761 // Simulate reflection
762
763 if (Swap) Swap = !Swap;
764
765 theStatus = FresnelReflection;
766
767 if ( theModel == unified && theFinish != polished )
768 ChooseReflection();
769
770 if ( theStatus == LambertianReflection ) {
771 DoReflection();
772 }
773 else if ( theStatus == BackScattering ) {
774 NewMomentum = -OldMomentum;
775 NewPolarization = -OldPolarization;
776 }
777 else {
778
779 PdotN = OldMomentum * theFacetNormal;
780 NewMomentum = OldMomentum - (2.*PdotN)*theFacetNormal;
781
782 if (sint1 > 0.0) { // incident ray oblique
783
784 E2_parl = Rindex2*E2_parl/Rindex1 - E1_parl;
785 E2_perp = E2_perp - E1_perp;
786 E2_total = E2_perp*E2_perp + E2_parl*E2_parl;
787 A_paral = NewMomentum.cross(A_trans);
788 A_paral = A_paral.unit();
789 E2_abs = std::sqrt(E2_total);
790 C_parl = E2_parl/E2_abs;
791 C_perp = E2_perp/E2_abs;
792
793 NewPolarization = C_parl*A_paral + C_perp*A_trans;
794
795 }
796
797 else { // incident ray perpendicular
798
799 if (Rindex2 > Rindex1) {
800 NewPolarization = - OldPolarization;
801 }
802 else {
803 NewPolarization = OldPolarization;
804 }
805
806 }
807 }
808 }
809 else { // photon gets transmitted
810
811 // Simulate transmission/refraction
812
813 Inside = !Inside;
814 Through = true;
815 theStatus = FresnelRefraction;
816
817 if (sint1 > 0.0) { // incident ray oblique
818
819 G4double alpha = cost1 - cost2*(Rindex2/Rindex1);
820 NewMomentum = OldMomentum + alpha*theFacetNormal;
821 NewMomentum = NewMomentum.unit();
822 PdotN = -cost2;
823 A_paral = NewMomentum.cross(A_trans);
824 A_paral = A_paral.unit();
825 E2_abs = std::sqrt(E2_total);
826 C_parl = E2_parl/E2_abs;
827 C_perp = E2_perp/E2_abs;
828
829 NewPolarization = C_parl*A_paral + C_perp*A_trans;
830
831 }
832 else { // incident ray perpendicular
833
834 NewMomentum = OldMomentum;
835 NewPolarization = OldPolarization;
836
837 }
838 }
839 }
840
841 OldMomentum = NewMomentum.unit();
842 OldPolarization = NewPolarization.unit();
843
844 if (theStatus == FresnelRefraction) {
845 Done = (NewMomentum * theGlobalNormal <= 0.0);
846 }
847 else {
848 Done = (NewMomentum * theGlobalNormal >= 0.0);
849 }
850
851 } while (!Done);
852
853 if (Inside && !Swap) {
854 if( theFinish == polishedbackpainted ||
855 theFinish == groundbackpainted ) {
856
857 if( !G4BooleanRand(theReflectivity) ) {
858 DoAbsorption();
859 }
860 else {
861 if (theStatus != FresnelRefraction ) {
862 theGlobalNormal = -theGlobalNormal;
863 }
864 else {
865 Swap = !Swap;
866 G4Swap(Material1,Material2);
867 G4Swap(&Rindex1,&Rindex2);
868 }
869 if ( theFinish == groundbackpainted )
870 theStatus = LambertianReflection;
871
872 DoReflection();
873
874 theGlobalNormal = -theGlobalNormal;
875 OldMomentum = NewMomentum;
876
877 goto leap;
878 }
879 }
880 }
881}
882
883// GetMeanFreePath
884// ---------------
885//
886G4double G4OpBoundaryProcess::GetMeanFreePath(const G4Track& ,
887 G4double ,
888 G4ForceCondition* condition)
889{
890 *condition = Forced;
891
892 return DBL_MAX;
893}
894
895G4double G4OpBoundaryProcess::GetIncidentAngle()
896{
897 G4double PdotN = OldMomentum * theFacetNormal;
898 G4double magP= OldMomentum.mag();
899 G4double magN= theFacetNormal.mag();
900 G4double incidentangle = pi - std::acos(PdotN/(magP*magN));
901
902 return incidentangle;
903}
904
905G4double G4OpBoundaryProcess::GetReflectivity(G4double E1_perp,
906 G4double E1_parl,
907 G4double incidentangle,
908 G4double RealRindex,
909 G4double ImaginaryRindex)
910{
911
912 G4complex Reflectivity, Reflectivity_TE, Reflectivity_TM;
913 G4complex N(RealRindex, ImaginaryRindex);
914 G4complex CosPhi;
915
916 G4complex u(1,0); //unit number 1
917
918 G4complex numeratorTE; // E1_perp=1 E1_parl=0 -> TE polarization
919 G4complex numeratorTM; // E1_parl=1 E1_perp=0 -> TM polarization
920 G4complex denominatorTE, denominatorTM;
921 G4complex rTM, rTE;
922
923 // Following two equations, rTM and rTE, are from: "Introduction To Modern
924 // Optics" written by Fowles
925
926 CosPhi=std::sqrt(u-((std::sin(incidentangle)*std::sin(incidentangle))/(N*N)));
927
928 numeratorTE = std::cos(incidentangle) - N*CosPhi;
929 denominatorTE = std::cos(incidentangle) + N*CosPhi;
930 rTE = numeratorTE/denominatorTE;
931
932 numeratorTM = N*std::cos(incidentangle) - CosPhi;
933 denominatorTM = N*std::cos(incidentangle) + CosPhi;
934 rTM = numeratorTM/denominatorTM;
935
936 // This is my calculaton for reflectivity on a metalic surface
937 // depending on the fraction of TE and TM polarization
938 // when TE polarization, E1_parl=0 and E1_perp=1, R=abs(rTE)^2 and
939 // when TM polarization, E1_parl=1 and E1_perp=0, R=abs(rTM)^2
940
941 Reflectivity_TE = (rTE*conj(rTE))*(E1_perp*E1_perp)
942 / (E1_perp*E1_perp + E1_parl*E1_parl);
943 Reflectivity_TM = (rTM*conj(rTM))*(E1_parl*E1_parl)
944 / (E1_perp*E1_perp + E1_parl*E1_parl);
945 Reflectivity = Reflectivity_TE + Reflectivity_TM;
946
947 do {
948 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TE))iTE = -1;
949 if(G4UniformRand()*real(Reflectivity) > real(Reflectivity_TM))iTM = -1;
950 } while(iTE<0&&iTM<0);
951
952 return real(Reflectivity);
953
954}
Note: See TracBrowser for help on using the repository browser.