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

Last change on this file since 1038 was 1007, checked in by garnier, 17 years ago

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