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

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

update to geant4.9.2
File size: 31.1 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        }
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
226        G4LogicalSurface* Surface = NULL;
227
228        Surface = G4LogicalBorderSurface::GetSurface
229                  (pPreStepPoint ->GetPhysicalVolume(),
230                   pPostStepPoint->GetPhysicalVolume());
231
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
257        if (Surface) OpticalSurface = 
258           dynamic_cast <G4OpticalSurface*> (Surface->GetSurfaceProperty());
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
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
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
593             // Comment out DoAbsorption if you wish to have
594             //       Transmission instead of Absorption
595
596             DoAbsorption();
597             theStatus = Absorption;
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;
675              G4SwapPtr(Material1,Material2);
676              G4SwapObj(&Rindex1,&Rindex2);
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;
882                   G4SwapPtr(Material1,Material2);
883                   G4SwapObj(&Rindex1,&Rindex2);
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.