FigureDayltgCoeffsAtPointsForSunPosition Subroutine

SUBROUTINE FigureDayltgCoeffsAtPointsForSunPosition(ZoneNum, iRefPoint, iXelement, NWX, iYelement, NWY, WinEl,  &
                            iWin, iWin2, iHour,iSunPos,SkyObstructionMult,  RWIN2, Ray, PHRAY, LSHCAL,  &
                            InShelfSurf, COSB, ObTrans,TVISB,DOMEGA, ICtrl, ShType, BlNum, THRAY, WNORM2, ExtWinType,  &
                            IConst, AZVIEW, RREF2, LoopWin, IHitIntObs,IHitExtObs, CalledFrom,   &
                            TVISIntWin, TVISIntWinDisk,                                               &
                            MapNum, MapWindowSolidAngAtRefPtWtd)

          ! SUBROUTINE INFORMATION:
          !       AUTHOR         B. Griffith
          !       DATE WRITTEN   November 2012, refactor from legacy code by Fred Winklemann
          !       MODIFIED       na
          !       RE-ENGINEERED  na

          ! PURPOSE OF THIS SUBROUTINE:
          ! collect code for calculations sun position aspects for daylighting coefficients

          ! METHODOLOGY EMPLOYED:
          ! switch as need to serve both reference points and map points based on calledFrom

          ! REFERENCES:
          ! na

          ! USE STATEMENTS:
  USE General,             ONLY: POLYF, BlindBeamBeamTrans
  USE DataSystemVariables, ONLY: DetailedSkyDiffuseAlgorithm
  USE SolarReflectionManager, ONLY: SolReflRecSurf

  IMPLICIT NONE ! Enforce explicit typing of all variables in this routine

          ! SUBROUTINE ARGUMENT DEFINITIONS:
  INTEGER, INTENT(IN)  :: ZoneNum
  INTEGER, INTENT(IN)  :: iRefPoint
  INTEGER, INTENT(IN)  :: iXelement
  INTEGER, Intent(IN)  :: NWX                          ! Number of window elements in x direction for dayltg calc
  INTEGER, Intent(IN)  :: iYelement
  INTEGER, Intent(IN)  :: NWY                          ! Number of window elements in y direction for dayltg calc
  INTEGER, Intent(IN)  :: WinEl                        ! Current window element counter
  INTEGER, INTENT(IN)  :: iWin
  INTEGER, INTENT(IN)  :: iWin2
  INTEGER, INTENT(IN)  :: iHour
  INTEGER, INTENT(INOUT)  :: iSunPos
  REAL(r64), INTENT(IN) :: SkyObstructionMult
  REAL(r64), DIMENSION(3),INTENT(IN) :: RWIN2  ! Center of a window element for TDD:DOME (if exists) in abs coord sys
  REAL(r64), DIMENSION(3),INTENT(IN) :: RAY    ! Unit vector along ray from reference point to window element
  REAL(r64),INTENT(IN) :: PHRAY                ! Altitude of ray from reference point to window element (radians)
  INTEGER, INTENT(IN)  :: LSHCAL               ! Interior shade calculation flag:  0=not yet calculated, 1=already calculated
  INTEGER, INTENT(IN)  :: InShelfSurf          ! Inside daylighting shelf surface number
  REAL(r64), INTENT(IN) :: COSB                ! Cosine of angle between window outward normal and ray from reference point
                                               ! to window element
  REAL(r64), INTENT(IN) :: ObTrans             ! Product of solar transmittances of exterior obstructions hit by ray
                                               ! from reference point through a window element
  REAL(r64), INTENT(IN) :: TVISB               ! Visible transmittance of window for COSB angle of incidence (times light well
                                               !   efficiency, if appropriate)
  REAL(r64), INTENT(IN)  :: DOMEGA             ! Solid angle subtended by window element wrt reference point (steradians)
  INTEGER, INTENT(IN)   :: ICtrl               ! Window control counter
  INTEGER, INTENT(IN)   :: ShType                       ! Window shading type
  INTEGER, INTENT(IN)   :: BlNum                        ! Window blind number
  REAL(r64), INTENT(IN) :: THRAY                        ! Azimuth of ray from reference point to window element (radians)
  REAL(r64), DIMENSION(3), INTENT(IN) :: WNORM2                    ! Unit vector normal to window
  INTEGER, INTENT(IN)   :: ExtWinType                   ! Exterior window type (InZoneExtWin, AdjZoneExtWin, NotInOrAdjZoneExtWin)
  INTEGER, INTENT(IN)   :: IConst                       ! Construction counter
  REAL(r64), INTENT(IN) :: AZVIEW                       ! Azimuth of view vector in absolute coord system for
                                            !  glare calculation (radians)
  REAL(r64), DIMENSION(3), INTENT(IN) :: RREF2  ! Location of virtual reference point in absolute coordinate system
  INTEGER, INTENT(IN)   :: LoopWin
  INTEGER, INTENT(IN)   :: IHitIntObs                   ! = 1 if interior obstruction hit, = 0 otherwise
  INTEGER, INTENT(IN)   :: IHitExtObs                   ! 1 if ray from ref pt to ext win hits an exterior obstruction
  INTEGER, INTENT(IN)   :: CalledFrom ! indicate  which type of routine called this routine
  REAL(r64), INTENT(INOUT) :: TVISIntWin    ! Visible transmittance of int win at COSBIntWin for light from ext win
  REAL(r64), INTENT(INOUT) :: TVISIntWinDisk    ! Visible transmittance of int win at COSBIntWin for sun
  INTEGER,   INTENT(IN),   OPTIONAL :: MapNum
  REAL(r64),  DIMENSION(:,:) , INTENT(IN), OPTIONAL :: MapWindowSolidAngAtRefPtWtd

          ! SUBROUTINE PARAMETER DEFINITIONS:
          ! na

          ! INTERFACE BLOCK SPECIFICATIONS:
          ! na

          ! DERIVED TYPE DEFINITIONS:
          ! na

          ! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
  REAL(r64), DIMENSION(3) :: RREF           ! Location of a reference point in absolute coordinate system
  REAL(r64), DIMENSION(4) :: XEDIRSK        ! Illuminance contribution from luminance element, sky-related
  REAL(r64) :: XEDIRSU                      ! Illuminance contribution from luminance element, sun-related
  REAL(r64), DIMENSION(4) :: XAVWLSK        ! Luminance of window element, sky-related
  REAL(r64), DIMENSION(3) :: RAYCOS         ! Unit vector from reference point to sun
  INTEGER   :: JB                           ! Slat angle counter
  REAL(r64) :: TransBmBmMult(MaxSlatAngs)   ! Beam-beam transmittance of isolated blind
  REAL(r64) :: TransBmBmMultRefl(MaxSlatAngs) ! As above but for beam reflected from exterior obstruction
  REAL(r64) :: ProfAng                      ! Solar profile angle on a window (radians)
  REAL(r64) :: POSFAC                       ! Position factor for a window element / ref point / view vector combination
  REAL(r64) :: XR                           ! Horizontal displacement ratio
  REAL(r64) :: YR                           ! Vertical displacement ratio
  INTEGER   :: IHIT                         ! Hit flag; =1 if ray from ref point thru window element hits
                                            !  an obstruction, =0 otherwise.

  REAL(r64) :: ObTransDisk                  ! Product of solar transmittances of exterior obstructions hit by ray
                                            ! from reference point to sun
  REAL(r64), DIMENSION(3) :: HP             ! Hit coordinates, if ray hits
  REAL(r64) :: LumAtHitPtFrSun              ! Luminance at hit point of obstruction by reflection of direct light from
                                            !  sun (cd/m2)
  INTEGER   :: ISky                         ! Sky type index: 1=clear, 2=clear turbid, 3=intermediate, 4=overcast

  REAL(r64) :: ELUM                         ! Sky or ground luminance (cd/m2)
  REAL(r64) :: DEDIR                        ! Illuminance contribution at reference point from window element (lux)
  REAL(r64) :: COSI                         ! Cosine of angle between direct sun and window outward normal
  INTEGER   :: IP                           ! IP=1 if ray passes thru window, =0 if not
  REAL(r64) :: TVISS                        ! Direct solar visible transmittance of window at given angle of incidence
                                            !  (times light well efficiency, if appropriate)
  REAL(r64) :: XAVWL                        ! XAVWL*TVISS is contribution of window luminance from solar disk (cd/m2)

  REAL(r64) :: SlatAng                      ! Blind slat angle (rad)
  INTEGER   :: Loop,Loop2                   ! DO loop indices
  INTEGER   :: NearestHitSurfNum            ! Surface number of nearest obstruction
  INTEGER   :: NearestHitSurfNumX           ! Surface number to use when obstruction is a shadowing surface
  REAL(r64), DIMENSION(3) :: NearestHitPt   ! Hit point of ray on nearest obstruction
  REAL(r64) :: SunObstructionMult           ! = 1.0 if sun hits a ground point; otherwise = 0.0
  REAL(r64) :: HorDis                       ! Distance between ground hit point and proj'n of center
                                            !  of window element onto ground (m)
  REAL(r64) :: Alfa,Beta                    ! Intermediate variables
  REAL(r64), DIMENSION(3) :: GroundHitPt    ! Coordinates of point that ray hits ground (m)
  INTEGER   :: IHitObs                      ! 1 if obstruction is hit; 0 otherwise
  REAL(r64), DIMENSION(3) :: ObsHitPt       ! Coordinates of hit point on an obstruction (m)
  INTEGER   :: ObsSurfNum                   ! Surface number of obstruction
  INTEGER   :: ObsConstrNum                 ! Construction number of obstruction
  REAL(r64) :: ObsVisRefl                   ! Visible reflectance of obstruction
  REAL(r64) :: SkyReflVisLum                ! Reflected sky luminance at hit point divided by

  INTEGER   :: RecSurfNum                   ! Receiving surface number
  INTEGER   :: ReflSurfNum,ReflSurfNumX     ! Reflecting surface number
  REAL(r64), DIMENSION(3) :: ReflNorm       ! Normal vector to reflecting surface
  REAL(r64) :: CosIncAngRefl                ! Cos of angle of incidence of beam on reflecting surface
  REAL(r64), DIMENSION(3) :: SunVecMir      ! Sun ray mirrored in reflecting surface
  REAL(r64) :: CosIncAngRec                 ! Cos of angle of incidence of reflected beam on receiving window
  INTEGER   :: IHitRefl                     ! 1 if ray hits reflecting surface; 0 otherwise
  REAL(r64), DIMENSION(3) :: HitPtRefl      ! Point that ray hits reflecting surface
  REAL(r64) :: ReflDistance                 ! Distance between ref pt and hit point on reflecting surf (m)
  INTEGER   :: IHitObsRefl                  ! > 0 if obstruction hit between ref pt and reflection point
  REAL(r64), DIMENSION(3) :: HitPtObs       ! Hit point on obstruction
  REAL(r64) :: ObsDistance                  ! Distance from ref pt to reflection point
  INTEGER   :: ReflSurfRecNum               ! Receiving surface number for a reflecting window
  REAL(r64) :: SpecReflectance              ! Specular reflectance of a reflecting surface
  REAL(r64) :: TVisRefl                     ! Bare window vis trans for reflected beam
                                            !  (times light well efficiency, if appropriate)
  INTEGER   :: ConstrNumRefl                ! Window construction number for a specularly reflecting shading surf
  REAL(r64) :: PHSUNrefl                    ! Altitude angle of reflected sun (radians)
  REAL(r64) :: THSUNrefl                    ! Azimuth anggle of reflected sun (radians)

  INTEGER   :: IHitIntWinDisk               ! 1 if ray from ref pt to sun passes thru an int window; 0 otherwise
  INTEGER   :: IHitIntObsDisk               ! 1 if ray from ref pt to sun hits an interior obstruction; 0 otherwise
  INTEGER   :: IHitExtObsDisk               ! 1 if ray from ref pt to sun hits an exterior obstruction; 0 otherwise

  INTEGER   :: IntWinDisk                   ! Surface loop index for finding int windows betw ref pt and sun
  REAL(r64), DIMENSION(3) :: HitPtIntWinDisk  ! Intersection point on an interior window for ray from ref pt to sun (m)
  INTEGER   :: IntWinDiskHitNum             ! Surface number of int window intersected by ray betw ref pt and sun
  REAL(r64) :: COSBIntWin                   ! Cos of angle between int win outward normal and ray betw ref pt and
                                            !  exterior window element or between ref pt and sun
  REAL(r64) :: TVisIntWinMult               ! Interior window vis trans multiplier for ext win in adjacent zone
  REAL(r64) :: TVisIntWinDiskMult           ! Interior window vis trans solar disk multiplier for ext win in adj zone
  REAL(r64) :: WindowSolidAngleDaylightPoint

  IF (SUNCOSHR(3,iHour) < SunIsUpValue) RETURN

  ISunPos = ISunPos + 1

 ! Altitude of sun (degrees)
  PHSUN = PHSUNHR(iHour)
  SPHSUN = SPHSUNHR(iHour)
  CPHSUN = CPHSUNHR(iHour)

  ! Azimuth of sun in absolute coord sys
  THSUN = THSUNHR(iHour)

            ! First time through, call routine to calculate inter-reflected illuminance
            ! at reference point and luminance of window with shade, screen or blind.

            ! Rob/TH - Not sure whether this call is necessary for interior zones with interior windows only.
!  new code would be -
! IF (LSHCAL == 1 .AND. ExtWinType /= AdjZoneExtWin) CALL DayltgInterReflectedIllum(ISunPos,IHR,ZoneNum,IWin2)
  IF (SurfaceWindow(IWin)%WindowModelType /= WindowBSDFModel) THEN
    IF (LSHCAL == 1) CALL DayltgInterReflectedIllum(iSunPos,iHour,ZoneNum,IWin2)
  ELSE
    IF (LSHCAL == 1) CALL DayltgInterReflectedIllumComplexFenestration(IWin2, WinEl, iHour, ZoneNum, iRefPoint, CalledFrom, MapNum)
    IF (COSB <= 0.d0) RETURN
    CALL DayltgDirectIllumComplexFenestration(IWin, WinEl, iHour, ZoneNum, iRefPoint, CalledFrom, MapNum)
    ! Call direct sun component only once since calculation is done for entire window
    IF (WinEl == (NWX * NWY)) THEN
      CALL DayltgDirectSunDiskComplexFenestration(IWin2, ZoneNum, iHour, loopwin, iRefPoint, WinEl, AZVIEW, CalledFrom, &
       MapNum, MapWindowSolidAngAtRefPtWtd)
    END IF
    RETURN
  END IF

            ! Daylighting shelf simplification:  The shelf completely blocks all view of the window,
            ! only interrelflected illumination is allowed (see DayltgInterReflectedIllum above).
            ! Everything else in this loop has to do with direct luminance from the window.
  IF (InShelfSurf > 0) RETURN

  IF (COSB <= 0.d0) RETURN

  XEDIRSK = 0.d0
  XEDIRSU = 0.d0
  XAVWLSK = 0.d0

  ! Add contribution of this window element to glare and to
  ! direct illuminance at reference point

  ! The I,J,K indices for sky and sun components of direct illuminance
  ! (EDIRSK, EDIRSU) and average window luminance (AVWLSK, AVWLSU) are:
  ! I=1 for clear sky, =2 Clear turbid, =3 Intermediate, =4 Overcast;
  ! J=1 for bare window, =2 for window with shade or fixed slat-angle blind;
  !  = 2,3,...,MaxSlatAngs+1 for window with variable slat-angle blind;
  ! K = sun position index.

  ! ----- CASE I -- BARE WINDOW (no shading device)

  ! Beam solar and sky solar reflected from nearest obstruction.
  ! In the following IHitIntObs == 0  ==> no interior obstructions hit, and
  !                  IHitExtObs == 1  ==> one or more exterior obstructions hit.
  IF(CalcSolRefl .AND. IHitIntObs == 0 .AND. IHitExtObs == 1) THEN
    ! One or more exterior obstructions was hit; get contribution of reflection
    ! from nearest obstruction.
    ! Find obstruction whose hit point is closest to this ray's window element
    CALL DayltgClosestObstruction(RWIN2,Ray,NearestHitSurfNum,NearestHitPt)
    IF(NearestHitSurfNum > 0) THEN

      ! Beam solar reflected from nearest obstruction

      CALL DayltgSurfaceLumFromSun(iHour,Ray,NearestHitSurfNum,NearestHitPt,LumAtHitPtFrSun)
      AVWLSU(1,iHour) = AVWLSU(1,iHour) + LumAtHitPtFrSun * TVISB
      IF (PHRAY >= 0.d0) EDIRSU(1,iHour) = EDIRSU(1,iHour) + LumAtHitPtFrSun * DOMEGA * RAY(3) * TVISB

      ! Sky solar reflected from nearest obstruction

      ObsConstrNum = Surface(NearestHitSurfNum)%Construction
      IF(ObsConstrNum > 0) THEN
        ! Exterior building surface is nearest hit
        IF(.NOT.Construct(ObsConstrNum)%TypeIsWindow) THEN
          ! Obstruction is not a window, i.e., is an opaque surface
          ObsVisRefl = 1.0d0 - Material(Construct(ObsConstrNum)%LayerPoint(1))%AbsorpVisible
        ELSE
          ! Obstruction is a window; assume it is bare
          IF(SurfaceWindow(NearestHitSurfNum)%StormWinFlag==1)  &
            ObsConstrNum = Surface(NearestHitSurfNum)%StormWinConstruction
          ObsVisRefl = Construct(ObsConstrNum)%ReflectVisDiffFront
        END IF
      ELSE
        ! Shadowing surface is nearest hit
        IF(Surface(NearestHitSurfNum)%Shelf > 0) THEN
          ! This is a daylighting shelf, for which reflection is separately calculated
          ObsVisRefl = 0.d0
        ELSE
          ObsVisRefl = Surface(NearestHitSurfNum)%ShadowSurfDiffuseVisRefl
          IF(Surface(NearestHitSurfNum)%ShadowSurfGlazingConstruct > 0) &
            ObsVisRefl = ObsVisRefl + Surface(NearestHitSurfNum)%ShadowSurfGlazingFrac *  &
            Construct(Surface(NearestHitSurfNum)%ShadowSurfGlazingConstruct)%ReflectVisDiffFront
        END IF
      END IF
      NearestHitSurfNumX = NearestHitSurfNum
      ! Each shadowing surface has a "mirror" duplicate surface facing in the opposite direction.
      ! The following gets the correct side of a shadowing surface for reflection.
      IF(Surface(NearestHitSurfNum)%ShadowingSurf) THEN
        IF(DOT_PRODUCT(RAY,Surface(NearestHitSurfNum)%OutNormVec) > 0.d0)  NearestHitSurfNumX = &
             NearestHitSurfNum + 1
      END IF
      IF (.not. DetailedSkyDiffuseAlgorithm .or. .not.  ShadingTransmittanceVaries .or.  &
          SolarDistribution == MinimalShadowing) THEN
        SkyReflVisLum = ObsVisRefl * Surface(NearestHitSurfNumX)%ViewFactorSky *  &
                 DifShdgRatioIsoSky(NearestHitSurfNumX) / Pi
      ELSE
        SkyReflVisLum = ObsVisRefl * Surface(NearestHitSurfNumX)%ViewFactorSky *  &
                 DifShdgRatioIsoSkyHRTS(NearestHitSurfNumX,iHour,1) / Pi
      ENDIF
      DO ISky = 1,4
        XAVWLSK(ISky) = GILSK(ISky,iHour) * SkyReflVisLum
        AVWLSK(ISky,1,iHour) = AVWLSK(ISky,1,iHour) + XAVWLSK(ISky) * TVISB
        IF (PHRAY >= 0.0d0) THEN
          XEDIRSK(ISky) = GILSK(ISky,iHour) * SkyReflVisLum * DOMEGA * RAY(3)
          EDIRSK(ISky,1,iHour) = EDIRSK(ISky,1,iHour) + XEDIRSK(ISky) * TVISB
        END IF
      END DO
    END IF
  END IF  ! End of check if solar reflection calculation is in effect

  IF(ObTrans > 1.d-6) THEN
    ! Ray did not hit an obstruction or the transmittance product of hit obstructions is non-zero.
    ! Contribution of sky or ground luminance in cd/m2
    IF (SurfaceWindow(IWin)%OriginalClass .EQ. SurfaceClass_TDD_Diffuser) THEN
      ! Make all transmitted light diffuse for a TDD with a bare diffuser
      DO ISKY = 1,4
        AVWLSK(ISKY,1,iHour) = AVWLSK(ISKY,1,iHour) + WLUMSK(ISKY,1,iHour)
        IF(ISky == 1) THEN
          AVWLSU(1,iHour) = AVWLSU(1,iHour) + WLUMSU(1,iHour)
          AVWLSUdisk(1,iHour) = AVWLSUdisk(1,iHour) + WLUMSUdisk(1,iHour)
        END IF
        IF(PHRAY > 0.d0) THEN
          EDIRSK(ISKY,1,iHour) = EDIRSK(ISKY,1,iHour) + WLUMSK(ISKY,1,iHour) * DOMEGA * RAY(3)
          IF(ISky == 1) EDIRSU(1,iHour) = EDIRSU(1,iHour) + WLUMSU(1,iHour) * DOMEGA * RAY(3)
        END IF
      END DO

    ELSE ! Bare window

      DO ISKY = 1,4
        IF (PHRAY > 0.d0) THEN
          ! Ray heads upward to sky
          ELUM = DayltgSkyLuminance(ISky,THRAY,PHRAY)
          XEDIRSK(ISKY) = ELUM * DOMEGA * RAY(3)
          DEDIR  = XEDIRSK(ISKY) * TVISB
          EDIRSK(ISKY,1,iHour) = EDIRSK(ISKY,1,iHour) + DEDIR * ObTrans
          AVWLSK(ISKY,1,iHour) = AVWLSK(ISKY,1,iHour) + ELUM * TVISB * ObTrans
          XAVWLSK(ISKY) = ELUM * ObTrans
        ELSE ! PHRAY <= 0.
          ! Ray heads downward to ground.
          ! Contribution from sky diffuse reflected from ground
          XAVWLSK(ISKY) = GILSK(ISKY,iHour) * (GndReflectanceForDayltg / PI) * ObTrans *  &
                                    SkyObstructionMult
          AVWLSK(ISKY,1,iHour) = AVWLSK(ISKY,1,iHour) + TVISB * XAVWLSK(ISKY)
          ! Contribution from beam solar reflected from ground (beam reaching ground point
          ! can be obstructed [SunObstructionMult < 1.0] if CalcSolRefl = .TRUE.)
          IF(ISky == 1) THEN
            SunObstructionMult = 1.d0
            IF(CalcSolRefl) THEN
              ! Coordinates of ground point hit by the ray
              Alfa = ACOS(-RAY(3))
              Beta = ATAN2(RAY(2),RAY(1))
              HorDis = (RWIN2(3)-GroundLevelZ)*TAN(Alfa)
              GroundHitPt(3) = GroundLevelZ
              GroundHitPt(1) = RWIN2(1) + HorDis*COS(Beta)
              GroundHitPt(2) = RWIN2(2) + HorDis*SIN(Beta)
              ! Sun reaches ground point if vector from this point to the sun is unobstructed
              IHitObs = 0
              DO ObsSurfNum = 1,TotSurfaces
                 IF(.NOT.Surface(ObsSurfNum)%ShadowSurfPossibleObstruction) CYCLE
                 CALL DayltgPierceSurface(ObsSurfNum,GroundHitPt,SunCosHr(1:3,iHour),IHitObs,ObsHitPt)
                 IF(IHitObs > 0) EXIT
              END DO
              IF(IHitObs > 0) SunObstructionMult = 0.d0
            END IF
            AVWLSU(1,iHour) = AVWLSU(1,iHour) +  &
                 TVISB * GILSU(iHour) * (GndReflectanceForDayltg / PI) * ObTrans * SunObstructionMult
          END IF  ! End of check if ISky = 1
        END IF  ! End of check if ray is going up or down
      END DO  ! End of loop over sky types
    END IF ! End of check if bare window or TDD:DIFFUSER
  END IF ! End of check if ObTrans > 1.E-6

  !
  ! Illuminance from beam solar (without interior reflection)
  !
  ! Just run this once on the last pass
  IF (iXelement == NWX .AND. iYelement == NWY) THEN ! Last pass

    ! Beam solar reaching reference point directly without exterior reflection

    ! Unit vector from ref. pt. to sun
    RAYCOS(1) = CPHSUN * COS(THSUN)
    RAYCOS(2) = CPHSUN * SIN(THSUN)
    RAYCOS(3) = SPHSUN

    ! Is sun on front side of exterior window?
    COSI = DOT_PRODUCT(WNORM2,RAYCOS)
    IF (COSI > 0.d0) THEN

      ! Does RAYCOS pass thru exterior window? HP is point that RAYCOS intersects window plane.
      CALL DayltgPierceSurface(IWin2,RREF2,RAYCOS,IP,HP)
      IHitIntObsDisk = 0
      IF (IP > 0) THEN
        IF(ExtWinType == InZoneExtWin) THEN
          ! Check for interior obstructions between reference point and HP.
          CALL DayltgHitInteriorObstruction(IWin2,RREF2,HP,IHitIntObsDisk)
        END IF
        ObTransDisk = 0.d0  ! Init value
        ! Init flag for vector from RP to sun passing through interior window
        IHitIntWinDisk = 0
        IF(ExtWinType == AdjZoneExtWin) THEN ! This block is for RPs in zones with interior windows
                                             ! adjacent to zones with exterior windows
          ! Does RAYCOS pass through interior window in zone containing RP?
          DO IntWinDisk = Zone(ZoneNum)%SurfaceFirst,Zone(ZoneNum)%SurfaceLast
            IF(Surface(IntWinDisk)%Class==SurfaceClass_Window .AND. Surface(IntWinDisk)%ExtBoundCond >= 1) THEN
              IF(Surface(Surface(IntWinDisk)%ExtBoundCond)%Zone == Surface(IWin2)%Zone) THEN
                CALL DayltgPierceSurface(IntWinDisk,RREF,RAYCOS,IHitIntWinDisk,HitPtIntWinDisk)
                IF(IHitIntWinDisk > 0) THEN
                  IntWinDiskHitNum = IntWinDisk
                  COSBIntWin = DOT_PRODUCT(Surface(IntWinDisk)%OutNormVec(1:3),RAYCOS)
                  IF(COSBIntWin <= 0.d0) THEN
                    IHitIntWinDisk = 0
                    IntWinDiskHitNum = 0
                    CYCLE
                  END IF
                  TVISIntWinDisk = POLYF(COSBIntWin,Construct(Surface(IntWinDisk)%Construction)%TransVisBeamCoef(1))
                  EXIT
                END IF
              END IF
            END IF
          END DO

          IF(IHitIntWinDisk == 0) THEN  ! Vector from RP to sun does not pass through interior window
            ObTransDisk = 0.d0
            IHIT = 1             !!fcw Is this needed?
          END IF

          ! Check for interior obstructions between ref point and interior window
          IHitIntObsDisk = 0
          IF(IHitIntWinDisk > 0) THEN
            CALL DayltgHitInteriorObstruction(IntWinDiskHitNum,RREF,HitPtIntWinDisk,IHitIntObsDisk)
            ! If no obstruction between RP and hit int win, check for obstruction
            ! between int win and ext win
            IF(IHitIntObsDisk == 0) THEN
              CALL DayltgHitBetWinObstruction(IntWinDiskHitNum,IWin2,HitPtIntWinDisk,HP,IHitIntObsDisk)
            ENDIF
          END IF
          IF(IHitIntObsDisk == 1) ObTransDisk = 0.d0
        END IF  ! case where RP is in zone with interior window adjacent to zone with exterior window

        IHitExtObsDisk = 0
        ! RJH 08-25-07 IHitIntWinDisk should not be reset to 0 here, and should be tested below.
        ! This is to correct logic flaw causing direct solar to reach adjacent zone refpt
        ! when vector to sun does not pass through interior window
        ! IHitIntWinDisk = 0
        IF(IHitIntObsDisk == 0) THEN  ! No interior obstruction was hit
          ! Net transmittance of exterior obstructions encountered by RAYCOS
          ! ObTransDisk = 1.0 will be returned if no exterior obstructions are hit.
          CALL DayltgHitObstruction(iHour,IWin2,RREF2,RAYCOS,ObTransDisk)
          IF(ObTransDisk < 1.d0) IHitExtObsDisk = 1
          ! RJH 08-26-07 However, if this is a case of interior window
          ! and vector to sun does not pass through interior window
          ! then reset ObTransDisk to 0.0 since it is the key test for adding
          ! contribution of sun to RP below.
          IF((ExtWinType == AdjZoneExtWin) .AND. (IHitIntWinDisk == 0)) THEN
            ObTransDisk = 0.d0
          END IF
        END IF

        ! PETER: need side wall mounted TDD to test this
        ! PETER: probably need to replace RREF2 with RWIN2
        ! PETER: need to check for interior obstructions too.

        IF (ObTransDisk > 1.d-6) THEN

          ! Sun reaches reference point;  increment illuminance.
          ! Direct normal illuminance is normalized to 1.0

          IF (SurfaceWindow(IWin)%OriginalClass .EQ. SurfaceClass_TDD_Diffuser) THEN
            ! No beam is transmitted.  Takes care of TDD with a bare diffuser and all types of blinds.
            TVISS = 0.d0
          ELSE
            ! Beam transmittance for bare window and all types of blinds
            TVISS = POLYF(COSI,Construct(IConst)%TransVisBeamCoef(1)) * SurfaceWindow(IWin)%GlazedFrac * &
                       SurfaceWindow(IWin)%LightWellEff
            IF(ExtWinType==AdjZoneExtWin.AND.IHitIntWinDisk==1)  &
              TVISS = TVISS * TVISIntWinDisk
          END IF

          EDIRSUdisk(1,iHour) = RAYCOS(3) * TVISS * ObTransDisk  ! Bare window

          TransBmBmMult = 0.d0
          IF (ShType == WSC_ST_ExteriorBlind .OR. ShType == WSC_ST_InteriorBlind .OR.   &
              ShType == WSC_ST_BetweenGlassBlind) THEN
            CALL ProfileAngle(IWin,RAYCOS,Blind(BlNum)%SlatOrientation,ProfAng)
            ! Contribution of beam passing through slats and reaching reference point
            DO JB = 1,MaxSlatAngs
              !IF (.NOT.SurfaceWindow(IWin)%MovableSlats .AND. JB > 1) EXIT
              IF (SurfaceWindow(IWin)%MovableSlats) THEN
                SlatAng = (JB - 1) * PI / (MaxSlatAngs - 1)
              ELSE
                SlatAng = Blind(BlNum)%SlatAngle * DegToRadians
              END IF
              TransBmBmMult(JB) = BlindBeamBeamTrans(ProfAng,SlatAng,Blind(BlNum)%SlatWidth, &
                                   Blind(BlNum)%SlatSeparation,Blind(BlNum)%SlatThickness)
              EDIRSUdisk(JB+1,iHour) = RAYCOS(3) * TVISS * TransBmBmMult(JB) * ObTransDisk

              ! do this only once for fixed slat blinds
              IF (.NOT.SurfaceWindow(IWin)%MovableSlats) EXIT
            END DO
          ELSE IF(ShType == WSC_ST_ExteriorScreen) THEN
!                          pass angle from sun to window normal here using PHSUN and THSUN from above and surface angles
!                          SunAltitudeToWindowNormalAngle = PHSUN - SurfaceWindow(IWin)%Phi
!                          SunAzimuthToWindowNormalAngle = THSUN - SurfaceWindow(IWin)%Theta
            CALL CalcScreenTransmittance(IWin, Phi=(PHSUN - SurfaceWindow(IWin)%Phi), &
                                             Theta=(THSUN - SurfaceWindow(IWin)%Theta))
            TransBmBmMult(1) = SurfaceScreens(SurfaceWindow(IWin)%ScreenNumber)%BmBmTrans
            EDIRSUdisk(2,iHour) = RAYCOS(3) * TVISS * TransBmBmMult(1) * ObTransDisk
          END IF

          ! Glare from solar disk

          ! Position factor for sun (note that AZVIEW is wrt y-axis and THSUN is wrt
          ! x-axis of absolute coordinate system.
          XR = TAN(ABS(PIOVR2 - AZVIEW - THSUN) + 0.001d0)
          YR = TAN(PHSUN + 0.001d0)
          POSFAC = DayltgGlarePositionFactor(XR,YR)

          SELECT CASE (CalledFrom)

          CASE ( CalledForRefPoint )
            WindowSolidAngleDaylightPoint = SurfaceWindow(IWin)%SolidAngAtRefPtWtd(iRefPoint)
          CASE ( CalledForMapPoint )
            WindowSolidAngleDaylightPoint = MapWindowSolidAngAtRefPtWtd(loopwin,iRefPoint)
          END SELECT


          IF (POSFAC /= 0.0d0.AND. WindowSolidAngleDaylightPoint > 0.000001d0) THEN
            ! Increment window luminance.  Luminance of solar disk (cd/m2)
            ! is 1.47*10^4*(direct normal solar illuminance) for direct normal solar
            ! illuminance in lux (lumens/m2). For purposes of calculating daylight factors
            ! direct normal solar illuminance = 1.0.
            ! Solid angle subtended by sun is 0.000068 steradians

            XAVWL = 14700.d0 * SQRT(0.000068d0 * POSFAC) *  REAL(NWX * NWY,r64) /  &
                     WindowSolidAngleDaylightPoint**0.8d0
            AVWLSUdisk(1,iHour) = XAVWL * TVISS * ObTransDisk  ! Bare window

            IF (ShType == WSC_ST_ExteriorBlind .OR. ShType == WSC_ST_InteriorBlind .OR.  &
                ShType == WSC_ST_BetweenGlassBlind) THEN
              DO JB = 1,MaxSlatAngs
                !IF (.NOT. SurfaceWindow(IWin)%MovableSlats .AND. JB > 1) EXIT
                AVWLSUdisk(JB+1,iHour) = XAVWL * TVISS * TransBmBmMult(JB) * ObTransDisk
                IF (.NOT. SurfaceWindow(IWin)%MovableSlats) EXIT
              END DO
            ELSE IF(ShType == WSC_ST_ExteriorScreen) THEN
              AVWLSUdisk(2,iHour) = XAVWL * TVISS * TransBmBmMult(1) * ObTransDisk
            END IF
          ENDIF ! Position Factor
        END IF ! Beam avoids all obstructions
      END IF ! Beam passes thru window
    END IF ! Sun on front side

    ! Beam solar reaching reference point after beam-beam (specular) reflection from
    ! an exterior surface

    IF(CalcSolRefl) THEN
      ! Receiving surface number corresponding this window
      RecSurfNum = Surface(IWin2)%ShadowSurfRecSurfNum
      IF (RecSurfNum > 0) THEN    ! interior windows do not apply
        IF(SolReflRecSurf(RecSurfNum)%NumPossibleObs > 0) THEN
          ! This window has associated obstructions that could reflect beam onto the window
          DO loop = 1,SolReflRecSurf(RecSurfNum)%NumPossibleObs
            ReflSurfNum = SolReflRecSurf(RecSurfNum)%PossibleObsSurfNums(loop)
            ReflSurfNumX = ReflSurfNum
            ! Each shadowing surface has a "mirror" duplicate surface facing in the opposite direction.
            ! The following gets the correct side of a shadowing surface for reflection.
            IF(Surface(ReflSurfNum)%ShadowingSurf) THEN
              IF(DOT_PRODUCT(RAYCOS,Surface(ReflSurfNum)%OutNormVec) < 0.d0)  ReflSurfNumX = &
              ReflSurfNum + 1
            END IF
            ! Require that the surface can have specular reflection
            IF(Surface(ReflSurfNum)%Class == SurfaceClass_Window .OR. &
               Surface(ReflSurfNum)%ShadowSurfGlazingFrac > 0.d0) THEN
              ReflNorm = Surface(ReflSurfNumX)%OutNormVec
              ! Vector to sun that is mirrored in obstruction
              SunVecMir = RAYCOS - 2.0d0*DOT_PRODUCT(RAYCOS,ReflNorm)*ReflNorm
              ! Skip if reflecting surface is not sunlit
              IF(SunlitFrac(ReflSurfNumX,iHour,1) < 0.01d0) CYCLE
              ! Skip if altitude angle of mirrored sun is negative since reflected sun cannot
              ! reach reference point in this case
              IF(SunVecMir(3) <= 0.d0) CYCLE
              ! Cosine of incidence angle of reflected beam on window
              CosIncAngRec = DOT_PRODUCT(Surface(IWin2)%OutNormVec,SunVecMir)
              IF(CosIncAngRec <= 0.d0) CYCLE
              ! Does ray from ref. pt. along SunVecMir pass through window?
              CALL DayltgPierceSurface(IWin2,RREF2,SunVecMir,IP,HP)
              IF(IP == 0) CYCLE  ! Ray did not pass through window
              ! Check if this ray hits interior obstructions
              CALL DayltgHitInteriorObstruction(IWin2,RREF2,HP,IHit)
              IF(IHit > 0) CYCLE ! Interior obstruction was hit
              ! Does ray hit this reflecting surface?
              CALL DayltgPierceSurface(ReflSurfNum,RREF2,SunVecMir,IHitRefl,HitPtRefl)
              IF(IHitRefl == 0) CYCLE  ! Ray did not hit this reflecting surface
              ReflDistance = SQRT(DOT_PRODUCT(HitPtRefl-RREF2,HitPtRefl-RREF2))
              ! Is ray from ref. pt. to reflection point (HitPtRefl) obstructed?
              IHitObsRefl = 0
              DO loop2 = 1,SolReflRecSurf(RecSurfNum)%NumPossibleObs
                ObsSurfNum = SolReflRecSurf(RecSurfNum)%PossibleObsSurfNums(loop2)
                IF(ObsSurfNum == ReflSurfNum .OR. ObsSurfNum == Surface(ReflSurfNum)%BaseSurf) CYCLE
                CALL DayltgPierceSurface(ObsSurfNum,RREF2,SunVecMir,IHitObs,HitPtObs)
                IF(IHitObs > 0) THEN
                  ObsDistance = SQRT(DOT_PRODUCT(HitPtObs-RREF2,HitPtObs-RREF2))
                  IF(ObsDistance < ReflDistance) THEN
                    IHitObsRefl = 1
                    EXIT
                  END IF
                END IF
              END DO
              IF(IHitObsRefl > 0) CYCLE  ! Obstruct'n closer than reflect'n pt. was hit; go to next obstruction
              ! There is no obstruction for this ray between ref pt and hit pt on reflecting surface.
              ! See if ray from hit pt on reflecting surface to original (unmirrored) sun position is obstructed
              IHitObs = 0
              IF(Surface(ReflSurfNum)%Class == SurfaceClass_Window) THEN
                ! Reflecting surface is a window.
                ! Receiving surface number for this reflecting window.
                ReflSurfRecNum = Surface(ReflSurfNum)%ShadowSurfRecSurfNum
                IF(ReflSurfRecNum > 0) THEN
                  ! Loop over possible obstructions for this reflecting window
                  DO loop2 = 1,SolReflRecSurf(ReflSurfRecNum)%NumPossibleObs
                    ObsSurfNum = SolReflRecSurf(ReflSurfRecNum)%PossibleObsSurfNums(loop2)
                    CALL DayltgPierceSurface(ObsSurfNum,HitPtRefl,RAYCOS,IHitObs,HitPtObs)
                    IF(IHitObs > 0) EXIT
                  END DO
                END IF
              ELSE
                ! Reflecting surface is a building shade
                DO ObsSurfNum = 1, TotSurfaces
                  IF(.NOT.Surface(ObsSurfNum)%ShadowSurfPossibleObstruction) CYCLE
                  IF(ObsSurfNum == ReflSurfNum) CYCLE
                  CALL DayltgPierceSurface(ObsSurfNum,HitPtRefl,RAYCOS,IHitObs,HitPtObs)
                  IF(IHitObs > 0) EXIT
                END DO
              END IF  ! End of check if reflector is a window or shadowing surface

              IF(IHitObs > 0) CYCLE ! Obstruct'n hit between reflect'n hit point and sun; go to next obstruction

              ! No obstructions. Calculate reflected beam illuminance at ref. pt. from this reflecting surface.
              SpecReflectance = 0.0d0
              CosIncAngRefl = ABS(DOT_PRODUCT(RAYCOS,ReflNorm))
              IF(Surface(ReflSurfNum)%Class == SurfaceClass_Window) THEN
                ConstrNumRefl = Surface(ReflSurfNum)%Construction
                IF(SurfaceWindow(ReflSurfNum)%StormWinFlag==1) &
                  ConstrNumRefl = Surface(ReflSurfNum)%StormWinConstruction
                SpecReflectance = POLYF(ABS(CosIncAngRefl),Construct(ConstrNumRefl)%ReflSolBeamFrontCoef(1:6))
              END IF
              IF(Surface(ReflSurfNum)%ShadowingSurf  &
                 .AND.Surface(ReflSurfNum)%ShadowSurfGlazingConstruct > 0)  &
                 SpecReflectance = Surface(ReflSurfNum)%ShadowSurfGlazingFrac * POLYF(ABS(CosIncAngRefl), &
                       Construct(Surface(ReflSurfNum)%ShadowSurfGlazingConstruct)%ReflSolBeamFrontCoef(1:6))
              TVisRefl = POLYF(CosIncAngRec,Construct(IConst)%TransVisBeamCoef(1)) *   &
                                  SurfaceWindow(IWin)%GlazedFrac * SurfaceWindow(IWin)%LightWellEff
              EDIRSUdisk(1,iHour) = EDIRSUdisk(1,iHour) + SunVecMir(3) * SpecReflectance * TVisRefl  ! Bare window

              TransBmBmMultRefl = 0.0d0
              IF (ShType == WSC_ST_ExteriorBlind .OR. ShType == WSC_ST_InteriorBlind .OR.  &
                  ShType == WSC_ST_BetweenGlassBlind) THEN
                CALL ProfileAngle(IWin,SunVecMir,Blind(BlNum)%SlatOrientation,ProfAng)
                ! Contribution of reflected beam passing through slats and reaching reference point
                DO JB = 1,MaxSlatAngs
                  !IF (.NOT.SurfaceWindow(IWin)%MovableSlats .AND. JB > 1) EXIT
                  IF (SurfaceWindow(IWin)%MovableSlats) THEN
                    SlatAng = REAL((JB - 1),r64) * PI / REAL((MaxSlatAngs - 1),r64)
                  ELSE
                    SlatAng = Blind(BlNum)%SlatAngle * DegToRadians
                  END IF
                  TransBmBmMultRefl(JB) = BlindBeamBeamTrans(ProfAng,SlatAng,Blind(BlNum)%SlatWidth, &
                                     Blind(BlNum)%SlatSeparation,Blind(BlNum)%SlatThickness)
                  EDIRSUdisk(JB+1,iHour) = EDIRSUdisk(JB+1,iHour) +  &
                     SunVecMir(3) * SpecReflectance * TVisRefl * TransBmBmMultRefl(JB)

                  IF (.NOT.SurfaceWindow(IWin)%MovableSlats) EXIT
                END DO
              ELSE IF(ShType == WSC_ST_ExteriorScreen) THEN
!                             pass angle from sun to window normal here using PHSUN and THSUN from above and surface angles
!                             SunAltitudeToWindowNormalAngle = PHSUN - SurfaceWindow(IWin)%Phi
!                             SunAzimuthToWindowNormalAngle = THSUN - SurfaceWindow(IWin)%Theta
                CALL CalcScreenTransmittance(IWin, Phi=(PHSUN - SurfaceWindow(IWin)%Phi), &
                                                 Theta=(THSUN - SurfaceWindow(IWin)%Theta))
                TransBmBmMultRefl(1) = SurfaceScreens(SurfaceWindow(IWin)%ScreenNumber)%BmBmTrans
                EDIRSUdisk(2,iHour) = EDIRSUdisk(2,iHour) +  &
                     SunVecMir(3) * SpecReflectance * TVisRefl * TransBmBmMultRefl(1)
              END IF  ! End of check if window has a blind or screen

              ! Glare from reflected solar disk

              PHSUNrefl = SunVecMir(3)
              THSUNrefl = ATAN2(SunVecMir(2),SunVecMir(1))
              XR = TAN(ABS(PiOvr2 - AZVIEW - THSUNrefl) + 0.001d0)
              YR = TAN(PHSUNrefl + 0.001d0)
              POSFAC=DayltgGlarePositionFactor(XR,YR)
              IF(POSFAC /= 0.0d0.AND. SurfaceWindow(IWin)%SolidAngAtRefPtWtd(iRefPoint) > 0.000001d0) THEN
                XAVWL = 14700.0d0 * SQRT(0.000068d0 * POSFAC) *  REAL(NWX * NWY,r64) /  &
                       SurfaceWindow(IWin)%SolidAngAtRefPtWtd(iRefPoint)**0.8d0
                AVWLSUdisk(1,iHour) = AVWLSUdisk(1,iHour) + XAVWL * TVisRefl * SpecReflectance  ! Bare window
                IF (ShType == WSC_ST_ExteriorBlind .OR. ShType == WSC_ST_InteriorBlind .OR.  &
                    ShType == WSC_ST_BetweenGlassBlind) THEN
                  DO JB = 1,MaxSlatAngs
                    !IF(.NOT. SurfaceWindow(IWin)%MovableSlats .AND. JB > 1) EXIT
                    AVWLSUdisk(JB+1,iHour) = AVWLSUdisk(JB+1,iHour) +  &
                          XAVWL * TVisRefl * SpecReflectance * TransBmBmMultRefl(JB)
                    IF(.NOT. SurfaceWindow(IWin)%MovableSlats) EXIT
                  END DO
                ELSE IF(ShType == WSC_ST_ExteriorScreen) THEN
                    AVWLSUdisk(2,iHour) = AVWLSUdisk(2,iHour) +  &
                          XAVWL * TVisRefl * SpecReflectance * TransBmBmMultRefl(1)
                END IF
              END IF
            END IF  ! End of check that obstruction can specularly reflect
          END DO  ! End of loop over obstructions associated with this window

        END IF  ! End of check if this window has associated obstructions
      END IF   ! End of check to see if this is exterior type window
    END IF  ! End of check if exterior reflection calculation is in effect

  END IF ! Last pass

  IF((ICtrl > 0 .AND.  &
     (ShType ==WSC_ST_InteriorShade.OR.ShType ==WSC_ST_ExteriorShade.OR.ShType ==WSC_ST_BetweenGlassShade.OR. &
      ShType ==WSC_ST_InteriorBlind.OR.ShType ==WSC_ST_ExteriorBlind.OR.ShType ==WSC_ST_BetweenGlassBlind .OR. &
      ShType ==WSC_ST_ExteriorScreen)).OR.SurfaceWindow(IWin)%SolarDiffusing) THEN

    ! ----- CASE II -- WINDOW WITH SCREEN, SHADE, BLIND, OR DIFFUSING WINDOW

    ! Interior window visible transmittance multiplier for exterior window in adjacent zone
    TVisIntWinMult = 1.d0
    TVisIntWinDiskMult = 1.d0
    IF(Surface(IWin)%Zone /= ZoneNum) THEN
      TVisIntWinMult = TVisIntWin
      TVisIntWinDiskMult = TVisIntWinDisk
    END IF

    DO ISKY = 1,4
      DO JB = 1,MaxSlatAngs
        !IF (.NOT.SurfaceWindow(IWin)%MovableSlats .AND. JB > 1) EXIT
        AVWLSK(ISKY,JB+1,iHour) = AVWLSK(ISKY,JB+1,iHour) + WLUMSK(ISKY,JB+1,iHour)*TVisIntWinMult
        IF (ISky == 1) THEN
          AVWLSU(JB+1,iHour) = AVWLSU(JB+1,iHour) + WLUMSU(JB+1,iHour)*TVisIntWinMult
          AVWLSUdisk(JB+1,iHour) = AVWLSUdisk(JB+1,iHour) + WLUMSUdisk(JB+1,iHour)*TVisIntWinDiskMult
        END IF
        IF (PHRAY > 0.d0) THEN
          EDIRSK(ISKY,JB+1,iHour) = EDIRSK(ISKY,JB+1,iHour) + WLUMSK(ISKY,JB+1,iHour)*DOMEGA*RAY(3)*TVisIntWinMult
          IF (ISky == 1) EDIRSU(JB+1,iHour) = EDIRSU(JB+1,iHour) + WLUMSU(JB+1,iHour)*DOMEGA*RAY(3)*TVisIntWinMult
        END IF
        IF (.NOT.SurfaceWindow(IWin)%MovableSlats) EXIT
      END DO
    END DO
  END IF

  RETURN

END SUBROUTINE FigureDayltgCoeffsAtPointsForSunPosition