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

Last change on this file since 1132 was 1055, checked in by garnier, 15 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.