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

Last change on this file since 1109 was 1055, checked in by garnier, 17 years ago

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