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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

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