Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2025, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <cassert>
50 : #include <cmath>
51 : #include <cstdint>
52 :
53 : // ObjexxFCL Headers
54 : #include <ObjexxFCL/Array.functions.hh>
55 : #include <ObjexxFCL/ArrayS.functions.hh>
56 :
57 : // EnergyPlus Headers
58 : #include <EnergyPlus/Construction.hh>
59 : #include <EnergyPlus/Data/EnergyPlusData.hh>
60 : #include <EnergyPlus/DataEnvironment.hh>
61 : #include <EnergyPlus/DataHeatBalSurface.hh>
62 : #include <EnergyPlus/DataHeatBalance.hh>
63 : #include <EnergyPlus/DataLoopNode.hh>
64 : #include <EnergyPlus/DataShadowingCombinations.hh>
65 : #include <EnergyPlus/DataSurfaces.hh>
66 : #include <EnergyPlus/DataSystemVariables.hh>
67 : #include <EnergyPlus/DataViewFactorInformation.hh>
68 : #include <EnergyPlus/DataZoneEquipment.hh>
69 : #include <EnergyPlus/HeatBalanceSurfaceManager.hh>
70 : #include <EnergyPlus/Material.hh>
71 : #include <EnergyPlus/PierceSurface.hh>
72 : #include <EnergyPlus/Psychrometrics.hh>
73 : #include <EnergyPlus/ScheduleManager.hh>
74 : #include <EnergyPlus/TARCOGGassesParams.hh>
75 : #include <EnergyPlus/TARCOGMain.hh>
76 : #include <EnergyPlus/TARCOGParams.hh>
77 : #include <EnergyPlus/UtilityRoutines.hh>
78 : #include <EnergyPlus/Vectors.hh>
79 : #include <EnergyPlus/WindowComplexManager.hh>
80 : #include <EnergyPlus/ZoneTempPredictorCorrector.hh>
81 :
82 : namespace EnergyPlus {
83 :
84 : namespace WindowComplexManager {
85 :
86 : // Module containing the routines dealing with complex fenestration
87 :
88 : // MODULE INFORMATION:
89 : // AUTHOR Joe Klems
90 : // DATE WRITTEN ???
91 : // MODIFIED November 2011, Simon Vidanovic
92 : // RE-ENGINEERED na
93 :
94 : // PURPOSE OF THIS MODULE:
95 : // Initialize data for solar and thermal calculations and also performs thermal calculations for BSDF window
96 :
97 : using namespace DataVectorTypes;
98 : using namespace DataBSDFWindow;
99 : using namespace DataSurfaces; // , ONLY: TotSurfaces,TotWindows,Surface,SurfaceWindow !update this later
100 : using namespace DataHeatBalance;
101 : using namespace DataShadowingCombinations;
102 : using namespace Vectors;
103 :
104 801 : void InitBSDFWindows(EnergyPlusData &state)
105 : {
106 :
107 : // SUBROUTINE INFORMATION:
108 : // AUTHOR Joe Klems
109 : // DATE WRITTEN August 2011
110 : // MODIFIED na
111 : // RE-ENGINEERED na
112 :
113 : // PURPOSE OF THIS SUBROUTINE:
114 : // Set up the overall optical geometry for a BSDF window
115 :
116 : using namespace Vectors;
117 :
118 : int BaseSurf; // base surface number (used in finding back surface)
119 : int NumStates; // Local variable for the number of states
120 801 : Array1D<Real64> Thetas; // temp array holding theta values
121 801 : Array1D_int NPhis; // temp array holding number of phis for a given theta
122 801 : Array1D<Real64> V(3); // vector array
123 : Real64 VLen; // Length of vector array
124 : int NHold; // No. values in the Temporary array
125 :
126 : struct TempBasisIdx
127 : {
128 : // Members
129 : int Basis; // Basis no in basis table
130 : int State; // State in which basis first occurs
131 :
132 : // Default Constructor
133 0 : TempBasisIdx() : Basis(0), State(0)
134 : {
135 0 : }
136 : };
137 :
138 : // Object Data
139 801 : Array1D<TempBasisIdx> IHold; // Temporary array
140 :
141 801 : if (state.dataBSDFWindow->TotComplexFenStates <= 0) {
142 791 : return; // Nothing to do if no complex fenestration states
143 : }
144 : // Construct Basis List
145 10 : state.dataWindowComplexManager->BasisList.allocate(state.dataBSDFWindow->TotComplexFenStates);
146 :
147 : // Note: Construction of the basis list contains the assumption of identical incoming and outgoing bases in
148 : // that the complex fenestration state definition contains only one basis description, hence
149 : // assumes square property matrices. If this assumption were relaxed through change of the
150 : // definition or additional definition of a state type with non-square matrices, then the loop
151 : // below should be modified to enter both of the bases into the basis list.
152 :
153 24 : for (int IConst = state.dataBSDFWindow->FirstBSDF; IConst <= state.dataBSDFWindow->FirstBSDF + state.dataBSDFWindow->TotComplexFenStates - 1;
154 : ++IConst) {
155 14 : state.dataWindowComplexManager->MatrixNo = state.dataConstruction->Construct(IConst).BSDFInput.BasisMatIndex;
156 14 : if (state.dataWindowComplexManager->NumBasis == 0) {
157 10 : state.dataWindowComplexManager->NumBasis = 1;
158 10 : ConstructBasis(state, IConst, state.dataWindowComplexManager->BasisList(1));
159 : } else {
160 11 : for (int IBasis = 1; IBasis <= state.dataWindowComplexManager->NumBasis; ++IBasis) {
161 7 : if (state.dataWindowComplexManager->MatrixNo == state.dataWindowComplexManager->BasisList(IBasis).BasisMatIndex) {
162 0 : goto BsLoop_loop;
163 : }
164 : }
165 4 : ++state.dataWindowComplexManager->NumBasis;
166 4 : ConstructBasis(state, IConst, state.dataWindowComplexManager->BasisList(state.dataWindowComplexManager->NumBasis));
167 : }
168 14 : BsLoop_loop:;
169 : }
170 10 : state.dataWindowComplexManager->BasisList.redimension(state.dataWindowComplexManager->NumBasis);
171 : // Proceed to set up geometry for complex fenestration states
172 10 : state.dataBSDFWindow->ComplexWind.allocate(state.dataSurface->TotSurfaces); // Set up companion array to SurfaceWindow to hold window
173 : // geometry for each state. This is an allocatable array of
174 : // geometries for the window states but only the complex
175 : // fenestration surfaces will have the arrays allocated
176 : // Search Thru Surfaces for Complex Fenestration State references
177 : // This will define the first complex fenestration state for that window, others will follow if there are
178 : // control specifications
179 10 : state.dataWindowComplexManager->WindowList.allocate(state.dataSurface->TotSurfaces); // Temporary allocation
180 20 : state.dataWindowComplexManager->WindowStateList.allocate(state.dataBSDFWindow->TotComplexFenStates,
181 10 : state.dataSurface->TotSurfaces); // Temporary allocation
182 125 : for (int ISurf = 1; ISurf <= state.dataSurface->TotSurfaces; ++ISurf) {
183 115 : int IConst = state.dataSurface->Surface(ISurf).Construction;
184 115 : if (IConst == 0) {
185 0 : continue; // This is true for overhangs (Shading:Zone:Detailed)
186 : }
187 115 : if (!(state.dataConstruction->Construct(IConst).TypeIsWindow && (state.dataConstruction->Construct(IConst).WindowTypeBSDF))) {
188 91 : continue; // Only BSDF windows
189 : }
190 : // Simon Check: Thermal construction removed
191 : // ThConst = Construct(IConst)%BSDFInput%ThermalConstruction
192 24 : state.dataSurface->SurfWinWindowModelType(ISurf) = WindowModel::BSDF;
193 24 : state.dataHeatBal->AnyBSDF = true;
194 24 : ++state.dataWindowComplexManager->NumComplexWind;
195 24 : NumStates = 1;
196 24 : state.dataWindowComplexManager->WindowList(state.dataWindowComplexManager->NumComplexWind).NumStates =
197 : 1; // Having found the construction reference in
198 : // the Surface array defines the first state for this window
199 24 : state.dataWindowComplexManager->WindowList(state.dataWindowComplexManager->NumComplexWind).SurfNo = ISurf;
200 : // WindowList( NumComplexWind ).Azimuth = DegToRadians * Surface( ISurf ).Azimuth;
201 : // WindowList( NumComplexWind ).Tilt = DegToRadians * Surface( ISurf ).Tilt;
202 24 : state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).InitInc =
203 24 : state.dataWindowComplexManager->Calculate_Geometry;
204 24 : state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).InitTrn =
205 24 : state.dataWindowComplexManager->Calculate_Geometry;
206 24 : state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).CopyIncState = 0;
207 24 : state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).CopyTrnState = 0;
208 24 : state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).Konst = IConst;
209 : // Simon Check: ThermalConstruction assigned to current construction
210 : // WindowStateList(NumComplexWind, NumStates)%ThermConst = ThConst
211 72 : for (int I = 1; I <= state.dataWindowComplexManager->NumBasis; ++I) { // Find basis in Basis List
212 48 : if (state.dataConstruction->Construct(IConst).BSDFInput.BasisMatIndex == state.dataWindowComplexManager->BasisList(I).BasisMatIndex) {
213 24 : state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).IncBasisIndx =
214 : I; // Note: square property matrices
215 24 : state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).TrnBasisIndx =
216 : I; // assumption
217 : }
218 : }
219 24 : if (state.dataWindowComplexManager->WindowStateList(NumStates, state.dataWindowComplexManager->NumComplexWind).IncBasisIndx <= 0) {
220 0 : ShowFatalError(state, "Complex Window Init: Window Basis not in BasisList.");
221 : }
222 : }
223 : // Should now have a WindowList with dataWindowComplexManager. NumComplexWind entries containing all the complex fenestrations
224 : // with a first state defined for each.
225 : // * * *
226 : // Here a search should be made for control specifications, which will give additional states for
227 : // controlled complex fenestrations. These should be added to the dataWindowComplexManager.WindowStateList, and
228 : // WindowList( )%NumStates incremented for each window for which states are added.
229 : // Added states should have WindowStateList ( , )%InitInc set to Calculate_Geometry
230 : // * * *
231 :
232 : // At this point, we have a complete WindowList and WindowStateList, with dataWindowComplexManager. NumComplexWind
233 : // defined, and NumStates for each complex window defined
234 : // Now sort through the window list to see that geometry will only be done once for each
235 : // window, basis combination
236 : // Note: code below assumes identical incoming and outgoing bases; following code will
237 : // need revision if this assumption relaxed
238 :
239 34 : for (int IWind = 1; IWind <= state.dataWindowComplexManager->NumComplexWind; ++IWind) { // Search window list for repeated bases
240 24 : if (state.dataWindowComplexManager->WindowList(IWind).NumStates > 1) {
241 0 : IHold.allocate(state.dataWindowComplexManager->WindowList(IWind).NumStates);
242 0 : NHold = 1;
243 0 : IHold(1).State = 1;
244 0 : IHold(1).Basis = state.dataWindowComplexManager->WindowStateList(1, IWind).IncBasisIndx;
245 : // If the Mth new basis found is basis B in the basis list, and it
246 : // first occurs in the WindowStateList in state N, then IHold(M)%Basis=B
247 : // and IHold(M)%State=N
248 0 : for (int K = 1; K <= state.dataWindowComplexManager->NumBasis; ++K) {
249 0 : if (K > NHold) {
250 0 : break;
251 : }
252 0 : int KBasis = IHold(K).Basis;
253 0 : int J = IHold(K).State;
254 0 : state.dataWindowComplexManager->InitBSDFWindowsOnce = true;
255 0 : for (int I = J + 1; I <= state.dataWindowComplexManager->WindowList(IWind).NumStates;
256 : ++I) { // See if subsequent states have the same basis
257 0 : if ((state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).InitInc ==
258 0 : state.dataWindowComplexManager->Calculate_Geometry) &&
259 0 : (state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).IncBasisIndx ==
260 : KBasis)) {
261 : // Note: square property matrices (same inc & trn bases) assumption
262 : // If same incident and outgoing basis assumption removed, following code will need to
263 : // be extended to treat the two bases separately
264 0 : state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).InitInc =
265 0 : state.dataWindowComplexManager->Copy_Geometry;
266 0 : state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).InitTrn =
267 0 : state.dataWindowComplexManager->Copy_Geometry;
268 0 : state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).CopyIncState = J;
269 0 : state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).CopyTrnState = J;
270 0 : } else if (state.dataWindowComplexManager->InitBSDFWindowsOnce) {
271 0 : state.dataWindowComplexManager->InitBSDFWindowsOnce = false; // First occurrence of a different basis
272 0 : ++NHold;
273 0 : IHold(NHold).State = I;
274 0 : IHold(NHold).Basis = state.dataWindowComplexManager->WindowStateList(I, IWind).IncBasisIndx;
275 0 : state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).InitTrn =
276 0 : state.dataWindowComplexManager->Calculate_Geometry;
277 0 : state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).CopyIncState = 0;
278 0 : state.dataWindowComplexManager->WindowStateList(I, state.dataWindowComplexManager->NumComplexWind).CopyTrnState = 0;
279 : }
280 : }
281 : }
282 0 : IHold.deallocate();
283 : }
284 : }
285 :
286 : // Now go through window list and window state list and calculate or copy the
287 : // geometry information for each window, state
288 34 : for (int IWind = 1; IWind <= state.dataWindowComplexManager->NumComplexWind; ++IWind) {
289 24 : int ISurf = state.dataWindowComplexManager->WindowList(IWind).SurfNo;
290 24 : NumStates = state.dataWindowComplexManager->WindowList(IWind).NumStates;
291 : // ALLOCATE(SurfaceWindow( ISurf )%ComplexFen) !activate the BSDF window description
292 : // for this surface
293 24 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.NumStates = NumStates;
294 24 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State.allocate(NumStates); // Allocate space for the states
295 24 : state.dataBSDFWindow->ComplexWind(ISurf).NumStates = NumStates;
296 24 : state.dataBSDFWindow->ComplexWind(ISurf).Geom.allocate(NumStates); // Allocate space for the geometries
297 : // Azimuth = WindowList( IWind ).Azimuth;
298 : // Tilt = WindowList( IWind ).Tilt;
299 : // Get the number of back surfaces for this window
300 24 : BaseSurf = state.dataSurface->Surface(ISurf).BaseSurf; // ShadowComb is organized by base surface
301 24 : int NBkSurf = state.dataShadowComb->ShadowComb(BaseSurf).NumBackSurf;
302 24 : state.dataBSDFWindow->ComplexWind(ISurf).NBkSurf = NBkSurf;
303 : // Define the back surface directions
304 24 : state.dataBSDFWindow->ComplexWind(ISurf).sWinSurf.allocate(NBkSurf);
305 24 : state.dataBSDFWindow->ComplexWind(ISurf).sdotN.allocate(NBkSurf);
306 : // Define the unit vectors pointing from the window center to the back surface centers
307 180 : for (int KBkSurf = 1; KBkSurf <= NBkSurf; ++KBkSurf) {
308 156 : BaseSurf = state.dataSurface->Surface(ISurf).BaseSurf; // ShadowComb is organized by base surface
309 156 : int JSurf = state.dataShadowComb->ShadowComb(BaseSurf).BackSurf(KBkSurf); // these are all proper back surfaces
310 156 : V = state.dataSurface->Surface(JSurf).Centroid - state.dataSurface->Surface(ISurf).Centroid;
311 156 : VLen = magnitude(V);
312 : // Define the unit vector from the window center to the back
313 156 : state.dataBSDFWindow->ComplexWind(ISurf).sWinSurf(KBkSurf) = V / VLen;
314 : // surface center
315 : // Define the back surface cosine(incident angle)
316 156 : state.dataBSDFWindow->ComplexWind(ISurf).sdotN(KBkSurf) = dot(V, state.dataSurface->Surface(JSurf).OutNormVec) / VLen;
317 : }
318 48 : for (int IState = 1; IState <= NumStates; ++IState) {
319 : // The following assumes identical incoming and outgoing bases. The logic will need to be
320 : // redesigned if this assumption is relaxed
321 24 : int IConst = state.dataWindowComplexManager->WindowStateList(IState, IWind).Konst;
322 : // ThConst = WindowStateList ( IWind , IState )%ThermConst
323 24 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState).Konst = IConst;
324 : // SurfaceWindow(ISurf)%ComplexFen%State(IState)%ThermConst = ThConst
325 24 : if (state.dataWindowComplexManager->WindowStateList(IState, IWind).InitInc == state.dataWindowComplexManager->Calculate_Geometry) {
326 48 : state.dataBSDFWindow->ComplexWind(ISurf).Geom(IState).Inc = state.dataWindowComplexManager->BasisList(
327 48 : state.dataWindowComplexManager->WindowStateList(IState, IWind).IncBasisIndx); // Put in the basis structure from the BasisList
328 24 : state.dataBSDFWindow->ComplexWind(ISurf).Geom(IState).Trn =
329 48 : state.dataWindowComplexManager->BasisList(state.dataWindowComplexManager->WindowStateList(IState, IWind).TrnBasisIndx);
330 :
331 24 : SetupComplexWindowStateGeometry(state,
332 : ISurf,
333 : IState,
334 : IConst,
335 24 : state.dataBSDFWindow->ComplexWind(ISurf),
336 24 : state.dataBSDFWindow->ComplexWind(ISurf).Geom(IState),
337 24 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState));
338 : // Note--setting up the state geometry will include constructing outgoing basis/surface
339 : // maps and those incoming maps that will not depend on shading.
340 : } else {
341 0 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState) = state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(
342 0 : state.dataWindowComplexManager->WindowStateList(IState, IWind).CopyIncState); // Note this overwrites Konst
343 0 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState).Konst = IConst; // so it has to be put back
344 : // SurfaceWindow (ISurf )%ComplexFen%State(IState)%ThermConst = ThConst !same for ThermConst
345 0 : state.dataBSDFWindow->ComplexWind(ISurf).Geom(IState) =
346 0 : state.dataBSDFWindow->ComplexWind(ISurf).Geom(state.dataWindowComplexManager->WindowStateList(IState, IWind).CopyIncState);
347 : }
348 :
349 : } // State loop
350 : } // Complex Window loop
351 : // Allocate all beam-dependent complex fenestration quantities
352 34 : for (int IWind = 1; IWind <= state.dataWindowComplexManager->NumComplexWind; ++IWind) {
353 24 : int ISurf = state.dataWindowComplexManager->WindowList(IWind).SurfNo;
354 24 : NumStates = state.dataWindowComplexManager->WindowList(IWind).NumStates;
355 48 : for (int IState = 1; IState <= NumStates; ++IState) {
356 24 : AllocateCFSStateHourlyData(state, ISurf, IState);
357 : } // State loop
358 : } // Complex Window loop
359 3174 : }
360 :
361 24 : void AllocateCFSStateHourlyData(EnergyPlusData &state,
362 : int const iSurf, // Surface number
363 : int const iState // Complex fenestration state number
364 : )
365 : {
366 :
367 : // SUBROUTINE INFORMATION:
368 : // AUTHOR Simon Vidanovic
369 : // DATE WRITTEN May 2013
370 : // MODIFIED na
371 : // RE-ENGINEERED na
372 :
373 : // PURPOSE OF THIS SUBROUTINE:
374 : // Allocate hourly data arrays for complex fenestration state
375 :
376 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
377 : int NLayers; // Number of complex fenestration layers
378 : int NBkSurf; // Number of back surfaces
379 : int KBkSurf; // Back surfaces counter
380 :
381 24 : NLayers = state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).NLayers;
382 24 : NBkSurf = state.dataBSDFWindow->ComplexWind(iSurf).NBkSurf;
383 :
384 48 : state.dataBSDFWindow->ComplexWind(iSurf).Geom(iState).SolBmGndWt.allocate(
385 24 : 24, state.dataGlobal->TimeStepsInHour, state.dataBSDFWindow->ComplexWind(iSurf).Geom(iState).NGnd);
386 24 : state.dataBSDFWindow->ComplexWind(iSurf).Geom(iState).SolBmIndex.allocate(24, state.dataGlobal->TimeStepsInHour);
387 24 : state.dataBSDFWindow->ComplexWind(iSurf).Geom(iState).ThetaBm.allocate(24, state.dataGlobal->TimeStepsInHour);
388 24 : state.dataBSDFWindow->ComplexWind(iSurf).Geom(iState).PhiBm.allocate(24, state.dataGlobal->TimeStepsInHour);
389 24 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).WinDirHemiTrans.allocate(24, state.dataGlobal->TimeStepsInHour);
390 24 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).WinDirSpecTrans.allocate(24, state.dataGlobal->TimeStepsInHour);
391 24 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).WinBmGndTrans.allocate(24, state.dataGlobal->TimeStepsInHour);
392 24 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).WinBmFtAbs.allocate(24, state.dataGlobal->TimeStepsInHour, NLayers);
393 24 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).WinBmGndAbs.allocate(24, state.dataGlobal->TimeStepsInHour, NLayers);
394 24 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).WinToSurfBmTrans.allocate(24, state.dataGlobal->TimeStepsInHour, NBkSurf);
395 24 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).BkSurf.allocate(NBkSurf);
396 180 : for (KBkSurf = 1; KBkSurf <= NBkSurf; ++KBkSurf) {
397 312 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).BkSurf(KBkSurf).WinDHBkRefl.allocate(24,
398 156 : state.dataGlobal->TimeStepsInHour);
399 312 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState).BkSurf(KBkSurf).WinDirBkAbs.allocate(
400 156 : 24, state.dataGlobal->TimeStepsInHour, NLayers);
401 : }
402 24 : }
403 :
404 0 : void ExpandComplexState(EnergyPlusData &state,
405 : int const iSurf, // Surface number
406 : int const iConst // Construction number
407 : )
408 : {
409 :
410 : // SUBROUTINE INFORMATION:
411 : // AUTHOR Simon Vidanovic
412 : // DATE WRITTEN May 2013
413 : // MODIFIED Simon Vidanovic (July 2013)
414 : // RE-ENGINEERED na
415 :
416 : // PURPOSE OF THIS SUBROUTINE:
417 : // When complex fenestration is controlled by EMS, program does not know in advance how many states are assigned to
418 : // ceratin surface. This information can be obtain only at runtime. Purpose of this routine is to extend number of states
419 : // used by complex fenestration in case that is necessary.
420 :
421 : // Expands states by one
422 0 : int NumOfStates = state.dataSurface->SurfaceWindow(iSurf).ComplexFen.NumStates;
423 :
424 0 : state.dataBSDFWindow->ComplexWind(iSurf).Geom.redimension(NumOfStates + 1);
425 0 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State.redimension(NumOfStates + 1);
426 :
427 : // Do daylighting geometry only in case it is initialized. If daylighting is not used then no need to expand state for that
428 0 : if (state.dataBSDFWindow->ComplexWind(iSurf).DaylightingInitialized) {
429 0 : state.dataBSDFWindow->ComplexWind(iSurf).DaylghtGeom.redimension(NumOfStates + 1);
430 0 : state.dataBSDFWindow->ComplexWind(iSurf).DaylightingInitialized = false;
431 : } else {
432 0 : state.dataBSDFWindow->ComplexWind(iSurf).DaylghtGeom.allocate(NumOfStates + 1);
433 : }
434 :
435 : // Increase number of states and insert new state
436 0 : ++NumOfStates;
437 0 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.NumStates = NumOfStates;
438 0 : state.dataBSDFWindow->ComplexWind(iSurf).NumStates = NumOfStates;
439 :
440 0 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(NumOfStates).Konst = iConst;
441 :
442 : // load basis and setup window state geometry
443 0 : ConstructBasis(state, iConst, state.dataBSDFWindow->ComplexWind(iSurf).Geom(NumOfStates).Inc);
444 0 : ConstructBasis(state, iConst, state.dataBSDFWindow->ComplexWind(iSurf).Geom(NumOfStates).Trn);
445 :
446 0 : SetupComplexWindowStateGeometry(state,
447 : iSurf,
448 : NumOfStates,
449 : iConst,
450 0 : state.dataBSDFWindow->ComplexWind(iSurf),
451 0 : state.dataBSDFWindow->ComplexWind(iSurf).Geom(NumOfStates),
452 0 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(NumOfStates));
453 :
454 : // allocation of memory for hourly data can be performed only after window state geometry has been setup
455 0 : AllocateCFSStateHourlyData(state, iSurf, NumOfStates);
456 :
457 : // calculate static properties for complex fenestration
458 0 : CalcWindowStaticProperties(state,
459 : iSurf,
460 : NumOfStates,
461 0 : state.dataBSDFWindow->ComplexWind(iSurf),
462 0 : state.dataBSDFWindow->ComplexWind(iSurf).Geom(NumOfStates),
463 0 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(NumOfStates));
464 :
465 : // calculate hourly data from complex fenestration
466 0 : CFSShadeAndBeamInitialization(state, iSurf, NumOfStates);
467 0 : }
468 :
469 39540 : void CheckCFSStates(EnergyPlusData &state, int const iSurf) // Surface number
470 : {
471 :
472 : // SUBROUTINE INFORMATION:
473 : // AUTHOR Simon Vidanovic
474 : // DATE WRITTEN May 2013
475 : // MODIFIED na
476 : // RE-ENGINEERED na
477 :
478 : // PURPOSE OF THIS SUBROUTINE:
479 : // Check if there are new states available for complex fenestration and performs proper initialization
480 :
481 : int NumOfStates; // number of states for current surface
482 : bool StateFound; // variable to indicate if state has been found
483 : int i; // Local counter
484 : int CurrentCFSState;
485 :
486 39540 : StateFound = false;
487 39540 : CurrentCFSState = state.dataSurface->SurfaceWindow(iSurf).ComplexFen.CurrentState;
488 :
489 : // Check if EMS changed construction number
490 39540 : if (state.dataSurface->Surface(iSurf).Construction != state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(CurrentCFSState).Konst) {
491 :
492 : // If construction number changed then take new state
493 : // First search for existing states. Maybe state is already added in previous timestep
494 0 : NumOfStates = state.dataSurface->SurfaceWindow(iSurf).ComplexFen.NumStates;
495 0 : for (i = 1; i <= NumOfStates; ++i) {
496 0 : if (state.dataSurface->Surface(iSurf).Construction == state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(i).Konst) {
497 0 : StateFound = true;
498 0 : CurrentCFSState = i;
499 0 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.CurrentState = i;
500 : }
501 : }
502 : } else {
503 39540 : StateFound = true;
504 : }
505 :
506 : // If new state is not found in the list of current states, then create new one, initialize and make it active
507 39540 : if (!StateFound) {
508 0 : ExpandComplexState(state, iSurf, state.dataSurface->Surface(iSurf).Construction);
509 0 : CurrentCFSState = state.dataSurface->SurfaceWindow(iSurf).ComplexFen.NumStates;
510 0 : state.dataSurface->SurfaceWindow(iSurf).ComplexFen.CurrentState = CurrentCFSState;
511 : }
512 39540 : }
513 :
514 801 : void InitComplexWindows(EnergyPlusData &state)
515 : {
516 :
517 : // SUBROUTINE INFORMATION:
518 : // AUTHOR Linda Lawrie
519 : // DATE WRITTEN November 2012
520 : // MODIFIED na
521 : // RE-ENGINEERED na
522 :
523 : // PURPOSE OF THIS SUBROUTINE:
524 : // Extract simple init for Complex Windows
525 :
526 : // One-time initialization
527 801 : if (state.dataWindowComplexManager->InitComplexWindowsOnce) {
528 801 : state.dataWindowComplexManager->InitComplexWindowsOnce = false;
529 801 : InitBSDFWindows(state);
530 801 : CalcStaticProperties(state);
531 : }
532 801 : }
533 :
534 10829 : void UpdateComplexWindows(EnergyPlusData &state)
535 : {
536 :
537 : // SUBROUTINE INFORMATION:
538 : // AUTHOR Joe Klems
539 : // DATE WRITTEN August 2011
540 : // MODIFIED B. Griffith, Nov. 2012 revised for detailed timestep integration mode
541 : // RE-ENGINEERED na
542 :
543 : // PURPOSE OF THIS SUBROUTINE:
544 : // Performs the shading-dependent initialization of the Complex Fenestration data;
545 : // On first call, calls the one-time initializition
546 :
547 : // LOGICAL,SAVE :: Once =.TRUE. !Flag for insuring things happen once
548 : int NumStates; // Number of states for a given complex fen
549 : int ISurf; // Index for sorting thru Surface array
550 : int IState; // Index identifying the window state for a particular window
551 : int IWind; // Index identifying a window in the WindowList
552 :
553 10829 : if (state.dataWindowComplexManager->NumComplexWind == 0) {
554 10805 : return;
555 : }
556 :
557 24 : if (state.dataGlobal->KickOffSizing || state.dataGlobal->KickOffSimulation) {
558 0 : return;
559 : }
560 :
561 : // Shading-dependent initialization; performed once for each shading period
562 :
563 : // Initialize the geometric quantities
564 :
565 76 : for (IWind = 1; IWind <= state.dataWindowComplexManager->NumComplexWind; ++IWind) {
566 52 : ISurf = state.dataWindowComplexManager->WindowList(IWind).SurfNo;
567 52 : NumStates = state.dataBSDFWindow->ComplexWind(ISurf).NumStates;
568 104 : for (IState = 1; IState <= NumStates; ++IState) {
569 52 : CFSShadeAndBeamInitialization(state, ISurf, IState);
570 : } // State loop
571 : } // window loop
572 : }
573 :
574 52 : void CFSShadeAndBeamInitialization(EnergyPlusData &state,
575 : int const iSurf, // Window surface number
576 : int const iState // Window state number
577 : )
578 : {
579 :
580 : // SUBROUTINE INFORMATION:
581 : // AUTHOR Simon Vidanovic
582 : // DATE WRITTEN May 2013
583 : // MODIFIED na
584 : // RE-ENGINEERED na
585 :
586 : // PURPOSE OF THIS SUBROUTINE:
587 : // Calculates shading properties of complex fenestration
588 : // Refactoring from Klems code
589 :
590 : using namespace Vectors;
591 :
592 52 : Vector SunDir(0.0, 0.0, 1.0); // unit vector pointing toward sun (world CS)
593 52 : Vector HitPt(0.0, 0.0, 1.0); // vector location of ray intersection with a surface
594 :
595 52 : if (state.dataGlobal->KickOffSizing || state.dataGlobal->KickOffSimulation) {
596 0 : return;
597 : }
598 :
599 : int IncRay; // Index of incident ray corresponding to beam direction
600 : Real64 Theta; // Theta angle of incident ray corresponding to beam direction
601 : Real64 Phi; // Phi angle of incident ray corresponding to beam direction
602 : bool hit; // hit flag
603 : int TotHits; // hit counter
604 52 : auto &complexWindow(state.dataBSDFWindow->ComplexWind(iSurf));
605 52 : auto &complexWindowGeom(complexWindow.Geom(iState));
606 52 : auto &surfaceWindowState(state.dataSurface->SurfaceWindow(iSurf).ComplexFen.State(iState));
607 :
608 52 : if (!state.dataSysVars->DetailedSolarTimestepIntegration) {
609 52 : std::size_t lHT(0); // Linear index for ( Hour, TS )
610 52 : std::size_t lHTI(0); // Linear index for ( Hour, TS, I )
611 1300 : for (int Hour = 1; Hour <= 24; ++Hour) {
612 6864 : for (int TS = 1; TS <= state.dataGlobal->TimeStepsInHour; ++TS, ++lHT) { // [ lHT ] == ( Hour, TS )
613 5616 : SunDir = state.dataBSDFWindow->SUNCOSTS(TS, Hour);
614 5616 : Theta = 0.0;
615 5616 : Phi = 0.0;
616 5616 : if (state.dataBSDFWindow->SUNCOSTS(TS, Hour)(3) > DataEnvironment::SunIsUpValue) {
617 2833 : IncRay = FindInBasis(state, SunDir, RayIdentificationType::Front_Incident, iSurf, iState, complexWindowGeom.Inc, Theta, Phi);
618 2833 : complexWindowGeom.ThetaBm[lHT] = Theta;
619 2833 : complexWindowGeom.PhiBm[lHT] = Phi;
620 : } else {
621 2783 : complexWindowGeom.ThetaBm[lHT] = 0.0;
622 2783 : complexWindowGeom.PhiBm[lHT] = 0.0;
623 2783 : IncRay = 0; // sundown can't have ray incident on window
624 : }
625 5616 : if (IncRay > 0) { // Sun may be incident on the window
626 1675 : complexWindowGeom.SolBmIndex[lHT] = IncRay;
627 : } else { // Window can't be sunlit, set front incidence ray index to zero
628 3941 : complexWindowGeom.SolBmIndex[lHT] = 0;
629 : }
630 212976 : for (int I = 1, nGnd = complexWindowGeom.NGnd; I <= nGnd; ++I, ++lHTI) { // Gnd pt loop
631 207360 : TotHits = 0;
632 207360 : Vector const gndPt(complexWindowGeom.GndPt(I));
633 3184002 : for (int JSurf = 1, eSurf = state.dataSurface->TotSurfaces; JSurf <= eSurf; ++JSurf) {
634 : // the following test will cycle on anything except exterior surfaces and shading surfaces
635 5977366 : if (state.dataSurface->Surface(JSurf).HeatTransSurf &&
636 2988683 : state.dataSurface->Surface(JSurf).ExtBoundCond != ExternalEnvironment) {
637 1036924 : continue;
638 : }
639 : // skip surfaces that face away from the ground point
640 1951759 : if (dot(SunDir, state.dataSurface->Surface(JSurf).NewellSurfaceNormalVector) >= 0.0) {
641 1430842 : continue;
642 : }
643 : // Looking for surfaces between GndPt and sun
644 520917 : hit = PierceSurface(state, JSurf, gndPt, SunDir, HitPt);
645 520917 : if (hit) {
646 : // Are not going into the details of whether a hit surface is transparent
647 : // Since this is ultimately simply weighting the transmittance, so great
648 : // detail is not warranted
649 12041 : ++TotHits;
650 12041 : break;
651 : }
652 : }
653 207360 : if (TotHits > 0) {
654 12041 : complexWindowGeom.SolBmGndWt[lHTI] = 0.0; // [ lHTI ] == ( Hour, TS, I )
655 : } else {
656 195319 : complexWindowGeom.SolBmGndWt[lHTI] = 1.0; // [ lHTI ] == ( Hour, TS, I )
657 : }
658 207360 : } // Gnd pt loop
659 :
660 : // update window beam properties
661 5616 : CalculateWindowBeamProperties(state, iSurf, iState, complexWindow, complexWindowGeom, surfaceWindowState, Hour, TS);
662 : } // Timestep loop
663 : } // Hour loop
664 : } else { // detailed timestep integration
665 : std::size_t const lHT(
666 0 : complexWindowGeom.ThetaBm.index(state.dataGlobal->HourOfDay, state.dataGlobal->TimeStep)); // [ lHT ] == ( HourOfDay, TimeStep )
667 0 : SunDir = state.dataBSDFWindow->SUNCOSTS(state.dataGlobal->TimeStep, state.dataGlobal->HourOfDay);
668 0 : Theta = 0.0;
669 0 : Phi = 0.0;
670 0 : if (state.dataBSDFWindow->SUNCOSTS(state.dataGlobal->TimeStep, state.dataGlobal->HourOfDay)(3) > DataEnvironment::SunIsUpValue) {
671 0 : IncRay = FindInBasis(state, SunDir, RayIdentificationType::Front_Incident, iSurf, iState, complexWindowGeom.Inc, Theta, Phi);
672 0 : complexWindowGeom.ThetaBm[lHT] = Theta;
673 0 : complexWindowGeom.PhiBm[lHT] = Phi;
674 : } else {
675 0 : complexWindowGeom.ThetaBm[lHT] = 0.0;
676 0 : complexWindowGeom.PhiBm[lHT] = 0.0;
677 0 : IncRay = 0; // sundown can't have ray incident on window
678 : }
679 :
680 0 : if (IncRay > 0) { // Sun may be incident on the window
681 0 : complexWindowGeom.SolBmIndex[lHT] = IncRay;
682 : } else { // Window can't be sunlit, set front incidence ray index to zero
683 0 : complexWindowGeom.SolBmIndex[lHT] = 0.0;
684 : }
685 0 : std::size_t lHTI(complexWindowGeom.SolBmGndWt.index(
686 0 : state.dataGlobal->HourOfDay, state.dataGlobal->TimeStep, 1)); // Linear index for ( HourOfDay, TimeStep, I )
687 0 : for (int I = 1, nGnd = complexWindowGeom.NGnd; I <= nGnd; ++I, ++lHTI) { // Gnd pt loop
688 0 : TotHits = 0;
689 0 : Vector const gndPt(complexWindowGeom.GndPt(I));
690 0 : for (int JSurf = 1; JSurf <= state.dataSurface->TotSurfaces; ++JSurf) {
691 : // the following test will cycle on anything except exterior surfaces and shading surfaces
692 0 : if (state.dataSurface->Surface(JSurf).HeatTransSurf && state.dataSurface->Surface(JSurf).ExtBoundCond != ExternalEnvironment) {
693 0 : continue;
694 : }
695 : // skip surfaces that face away from the ground point
696 0 : if (dot(SunDir, state.dataSurface->Surface(JSurf).NewellSurfaceNormalVector) >= 0.0) {
697 0 : continue;
698 : }
699 : // Looking for surfaces between GndPt and sun
700 0 : hit = PierceSurface(state, JSurf, gndPt, SunDir, HitPt);
701 0 : if (hit) {
702 : // Are not going into the details of whether a hit surface is transparent
703 : // Since this is ultimately simply weighting the transmittance, so great
704 : // detail is not warranted
705 0 : ++TotHits;
706 0 : break;
707 : }
708 : }
709 0 : if (TotHits > 0) {
710 0 : complexWindowGeom.SolBmGndWt[lHTI] = 0.0; // [ lHTI ] == ( HourOfDay, TimeStep, I )
711 : } else {
712 0 : complexWindowGeom.SolBmGndWt[lHTI] = 1.0; // [ lHTI ] == ( HourOfDay, TimeStep, I )
713 : }
714 0 : } // Gnd pt loop
715 :
716 : // Update window beam properties
717 0 : CalculateWindowBeamProperties(
718 0 : state, iSurf, iState, complexWindow, complexWindowGeom, surfaceWindowState, state.dataGlobal->HourOfDay, state.dataGlobal->TimeStep);
719 : } // solar calculation mode, average over days or detailed
720 52 : }
721 :
722 5616 : void CalculateWindowBeamProperties(EnergyPlusData &state,
723 : int const ISurf, // Window surface number
724 : int const IState, // Window state number
725 : BSDFWindowGeomDescr const &Window, // Window Geometry
726 : BSDFGeomDescr const &Geom, // State Geometry
727 : BSDFStateDescr &State, // State Description
728 : int const Hour, // Hour number
729 : int const TS // Timestep number
730 : )
731 : {
732 :
733 : // SUBROUTINE INFORMATION:
734 : // AUTHOR Joe Klems
735 : // DATE WRITTEN August 2011
736 : // MODIFIED na
737 : // RE-ENGINEERED na
738 :
739 : // PURPOSE OF THIS SUBROUTINE:
740 : // Calculates those optical properties of all the Complex Fenestrations that
741 : // depend on the beam direction (hence, on hour and time step)
742 :
743 : // METHODOLOGY EMPLOYED:
744 : // Locate the bidirectional property matrices in the BSDFInput structure
745 : // and use them to calculate the desired average properties.
746 :
747 : using namespace Vectors;
748 :
749 : int IConst; // State construction number
750 : int I; // general purpose index--Back surface
751 : int J; // general purpose index--ray
752 : int JRay; // ray index number
753 : Real64 Theta;
754 : Real64 Phi;
755 : int JSurf; // gen purpose surface no
756 : int BaseSurf; // base surface no
757 : int M; // general purpose index--ray
758 : int L; // general purpose index--layer
759 : int KBkSurf; // general purpose index--back surface
760 : Real64 Sum1; // general purpose sum
761 : Real64 Sum2; // general purpose sum
762 : int IBm; // index of beam ray in incoming basis
763 : int BkIncRay; // index of sun dir in back incidence basis
764 : bool RegWindFnd; // flag for regular exterior back surf window
765 5616 : Array1D_int RegWinIndex; // bk surf nos of reg windows
766 5616 : int NRegWin(0); // no reg windows found as back surfaces
767 5616 : int KRegWin(0); // index of reg window as back surface
768 : Real64 Refl; // temporary reflectance
769 5616 : Array1D<Real64> Absorb; // temporary layer absorptance
770 :
771 : // Object Data
772 5616 : Vector SunDir; // current sun direction
773 :
774 5616 : IConst = state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState).Konst;
775 :
776 : // Begin calculation
777 : // Calculate the Transmittance from a given beam direction to a given zone surface
778 :
779 5616 : IBm = Geom.SolBmIndex(Hour, TS);
780 5616 : if (IBm <= 0.0) { // Beam cannot be incident on window for this Hour, TS
781 3941 : State.WinToSurfBmTrans(Hour, TS, {1, Window.NBkSurf}) = 0.0;
782 3941 : State.WinDirHemiTrans(Hour, TS) = 0.0;
783 3941 : State.WinDirSpecTrans(Hour, TS) = 0.0;
784 3941 : State.WinBmFtAbs(Hour, TS, {1, State.NLayers}) = 0.0;
785 : } else {
786 12084 : for (I = 1; I <= Window.NBkSurf; ++I) { // Back surface loop
787 10409 : Sum1 = 0.0;
788 159996 : for (J = 1; J <= Geom.NSurfInt(I); ++J) { // Ray loop
789 149587 : Sum1 +=
790 149587 : Geom.Trn.Lamda(Geom.SurfInt(J, I)) * state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(IBm, Geom.SurfInt(J, I));
791 : } // Ray loop
792 10409 : State.WinToSurfBmTrans(Hour, TS, I) = Sum1;
793 : } // Back surface loop
794 : // Calculate the directional-hemispherical transmittance
795 1675 : Sum1 = 0.0;
796 151262 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
797 149587 : Sum1 += Geom.Trn.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(IBm, J);
798 : }
799 1675 : State.WinDirHemiTrans(Hour, TS) = Sum1;
800 : // Calculate the directional specular transmittance
801 : // Note: again using assumption that Inc and Trn basis have same structure
802 1675 : State.WinDirSpecTrans(Hour, TS) = Geom.Trn.Lamda(IBm) * state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(IBm, IBm);
803 : // Calculate the layer front absorptance for beam radiation
804 6230 : for (L = 1; L <= State.NLayers; ++L) {
805 4555 : State.WinBmFtAbs(Hour, TS, L) = state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).FrtAbs(IBm, 1);
806 : }
807 : }
808 : // Calculate, for a given beam direction, the transmittance into the zone
809 : // for ground-reflected radiation (transmitted radiation assumed uniformly diffuse)
810 :
811 5616 : Sum1 = 0.0;
812 5616 : Sum2 = 0.0;
813 212976 : for (J = 1; J <= Geom.NGnd; ++J) { // Incident ray loop
814 207360 : JRay = Geom.GndIndex(J);
815 207360 : if (Geom.SolBmGndWt(Hour, TS, J) > 0.0) {
816 195319 : Sum2 += Geom.SolBmGndWt(Hour, TS, J) * Geom.Inc.Lamda(JRay);
817 23724254 : for (M = 1; M <= Geom.Trn.NBasis; ++M) { // Outgoing ray loop
818 23528935 : Sum1 += Geom.SolBmGndWt(Hour, TS, J) * Geom.Inc.Lamda(JRay) * Geom.Trn.Lamda(M) *
819 23528935 : state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(JRay, M);
820 : } // Outgoing ray loop
821 : }
822 : } // Indcident ray loop
823 5616 : if (Sum2 > 0.0) {
824 5561 : State.WinBmGndTrans(Hour, TS) = Sum1 / Sum2;
825 : } else {
826 55 : State.WinBmGndTrans(Hour, TS) = 0.0; // No unshaded ground => no transmittance
827 : }
828 :
829 : // Calculate, for a given beam direction, the layer front absorptance
830 : // for ground-reflected radiation
831 :
832 21216 : for (L = 1; L <= State.NLayers; ++L) { // layer loop
833 15600 : Sum1 = 0.0;
834 15600 : Sum2 = 0.0;
835 599280 : for (J = 1; J <= Geom.NGnd; ++J) { // Incident ray loop
836 583680 : JRay = Geom.GndIndex(J);
837 583680 : if (Geom.SolBmGndWt(Hour, TS, J) > 0.0) {
838 548444 : Sum2 += Geom.SolBmGndWt(Hour, TS, J) * Geom.Inc.Lamda(JRay);
839 548444 : Sum1 += Geom.SolBmGndWt(Hour, TS, J) * Geom.Inc.Lamda(JRay) *
840 548444 : state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).FrtAbs(JRay, 1);
841 : }
842 : } // Incident ray loop
843 15600 : if (Sum2 > 0.0) {
844 15435 : State.WinBmGndAbs(Hour, TS, L) = Sum1 / Sum2;
845 : } else {
846 165 : State.WinBmGndAbs(Hour, TS, L) = 0.0; // No unshaded ground => no absorptance
847 : }
848 : } // layer loop
849 :
850 : // Check the back surfaces for exterior windows
851 5616 : RegWindFnd = false;
852 5616 : NRegWin = 0.0;
853 5616 : RegWinIndex.allocate(Window.NBkSurf);
854 41760 : for (KBkSurf = 1; KBkSurf <= Window.NBkSurf; ++KBkSurf) {
855 36144 : BaseSurf = state.dataSurface->Surface(ISurf).BaseSurf; // ShadowComb is organized by base surface
856 36144 : JSurf = state.dataShadowComb->ShadowComb(BaseSurf).BackSurf(KBkSurf);
857 36144 : if (state.dataSurface->SurfWinWindowModelType(JSurf) == WindowModel::BSDF) {
858 8064 : continue;
859 : }
860 56160 : if (!(state.dataSurface->Surface(JSurf).Class == SurfaceClass::Window ||
861 28080 : state.dataSurface->Surface(JSurf).Class == SurfaceClass::GlassDoor)) {
862 28080 : continue;
863 : }
864 0 : if (!(state.dataSurface->Surface(JSurf).HeatTransSurf && state.dataSurface->Surface(JSurf).ExtBoundCond == ExternalEnvironment &&
865 0 : state.dataSurface->Surface(JSurf).ExtSolar)) {
866 0 : continue;
867 : }
868 : // Back surface is an exterior window or door
869 0 : RegWindFnd = true;
870 0 : ++NRegWin;
871 0 : RegWinIndex(NRegWin) = KBkSurf;
872 : }
873 5616 : if (RegWindFnd) {
874 0 : Absorb.allocate(State.NLayers);
875 0 : SunDir = state.dataBSDFWindow->SUNCOSTS(TS, Hour);
876 0 : BkIncRay = FindInBasis(state,
877 : SunDir,
878 : RayIdentificationType::Back_Incident,
879 : ISurf,
880 : IState,
881 0 : state.dataBSDFWindow->ComplexWind(ISurf).Geom(IState).Trn,
882 : Theta,
883 : Phi);
884 0 : if (BkIncRay > 0) {
885 : // Here calculate the back incidence properties for the solar ray
886 : // this does not say whether or not the ray can pass through the
887 : // back surface window and hit this one!
888 0 : Sum1 = 0.0;
889 0 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
890 0 : Sum1 += Geom.Trn.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.SolBkRefl(BkIncRay, J);
891 : }
892 0 : Refl = Sum1;
893 0 : for (L = 1; L <= State.NLayers; ++L) {
894 0 : Absorb(L) = state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).BkAbs(BkIncRay, 1);
895 : }
896 : } else {
897 : // solar ray can't be incident on back, so set properties equal to zero
898 0 : Refl = 0.0;
899 0 : for (L = 1; L <= State.NLayers; ++L) {
900 0 : Absorb(L) = 0.0;
901 : }
902 : }
903 0 : for (KRegWin = 1; KRegWin <= NRegWin; ++KRegWin) {
904 0 : KBkSurf = RegWinIndex(KRegWin);
905 0 : State.BkSurf(KBkSurf).WinDHBkRefl(Hour, TS) = Refl;
906 0 : for (L = 1; L <= State.NLayers; ++L) {
907 0 : State.BkSurf(KBkSurf).WinDirBkAbs(Hour, TS, L) = Absorb(L);
908 : }
909 : }
910 : }
911 5616 : if (allocated(Absorb)) {
912 0 : Absorb.deallocate();
913 : }
914 5616 : RegWinIndex.deallocate();
915 5616 : }
916 :
917 801 : void CalcStaticProperties(EnergyPlusData &state)
918 : {
919 :
920 : // SUBROUTINE INFORMATION:
921 : // AUTHOR Joe Klems
922 : // DATE WRITTEN <date_written>
923 : // MODIFIED na
924 : // RE-ENGINEERED na
925 :
926 : // PURPOSE OF THIS SUBROUTINE:
927 : // Calculates those optical properties of all the Complex Fenestrations that
928 : // do not depend on the beam direction (hence, on hour and time step)
929 :
930 : using namespace Vectors;
931 :
932 : // ISurf(0); // Index for sorting thru Surface array
933 : // // static int IConst( 0 ); // Index for accessing Construct array
934 : // IState(0); // Index identifying the window state for a particular window
935 : // IWind(0); // Index identifying a window in the WindowList
936 : // NumStates(0); // local copy of no of states
937 :
938 825 : for (int IWind = 1; IWind <= state.dataWindowComplexManager->NumComplexWind; ++IWind) {
939 24 : int ISurf = state.dataWindowComplexManager->WindowList(IWind).SurfNo;
940 24 : int NumStates = state.dataWindowComplexManager->WindowList(IWind).NumStates;
941 48 : for (int IState = 1; IState <= NumStates; ++IState) {
942 : // IConst = WindowStateList ( IWind , IState )%Konst
943 24 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState).Konst =
944 24 : state.dataWindowComplexManager->WindowStateList(IState, IWind).Konst;
945 24 : CalcWindowStaticProperties(state,
946 : ISurf,
947 : IState,
948 24 : state.dataBSDFWindow->ComplexWind(ISurf),
949 24 : state.dataBSDFWindow->ComplexWind(ISurf).Geom(IState),
950 24 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState));
951 : }
952 : }
953 801 : }
954 :
955 14 : void CalculateBasisLength(EnergyPlusData &state,
956 : BSDFWindowInputStruct const &Input, // BSDF data input struct for this construction
957 : int const IConst, // Construction number of input
958 : int &NBasis // Calculated Basis length
959 : )
960 : {
961 :
962 : // SUBROUTINE INFORMATION:
963 : // AUTHOR Joe Klems
964 : // DATE WRITTEN August 2011
965 : // MODIFIED na
966 : // RE-ENGINEERED na
967 :
968 : // PURPOSE OF THIS SUBROUTINE:
969 : // Calculates the basis length for a Window6 Non-Symmetric or Axisymmetric basis
970 : // from the input basis matrix
971 :
972 14 : if (Input.BasisMatNcols == 1) {
973 : // Axisymmetric basis, No. rows is no. of thetas = basis length
974 0 : NBasis = Input.BasisMatNrows;
975 0 : return;
976 : }
977 14 : NBasis = 1;
978 106 : for (int I = 2; I <= Input.BasisMatNrows; ++I) {
979 92 : NBasis += std::floor(state.dataConstruction->Construct(IConst).BSDFInput.BasisMat(2, I) + 0.001);
980 : }
981 : }
982 :
983 14 : void ConstructBasis(EnergyPlusData &state,
984 : int const IConst, // Index for accessing Construct array
985 : BasisStruct &Basis)
986 : {
987 : // SUBROUTINE INFORMATION:
988 : // AUTHOR Joe Klems
989 : // DATE WRITTEN June 2011
990 : // MODIFIED na
991 : // RE-ENGINEERED na
992 :
993 : // PURPOSE OF THIS SUBROUTINE:
994 : // Set up a basis from the matrix information pointed to in Construction by ICons
995 :
996 14 : int I(0); // general purpose index
997 14 : int J(0); // general purpose index
998 14 : int NThetas(0); // Current number of theta values
999 14 : int NumElem(0); // Number of elements in current basis
1000 14 : int ElemNo(0); // Current basis element number
1001 : int MaxNPhis; // Max no of NPhis for any theta
1002 14 : Real64 Theta(0.0); // Current theta value
1003 14 : Real64 Phi(0.0); // Current phi value
1004 14 : Real64 DTheta(0.0); // Increment for theta value (Window6 type input)
1005 14 : Real64 DPhi(0.0); // Increment for phi value (Window6 type input)
1006 14 : Real64 HalfDTheta(0.0); // Half-width of all theta bins except first and last (W6 input)
1007 14 : Real64 Lamda(0.0); // Current 'Lamda' value (element weight)
1008 14 : Real64 SolAng(0.0); // Current element solid angle
1009 14 : Real64 NextTheta(0.0); // Next theta in the W6 basis after current
1010 14 : Real64 LastTheta(0.0); // Previous theta in the W6 basis before current
1011 14 : Real64 LowerTheta(0.0); // Lower theta boundary of the element
1012 14 : Real64 UpperTheta(0.0); // Upper theta boundary of the element
1013 14 : Array1D<Real64> Thetas; // temp array holding theta values
1014 14 : Array1D_int NPhis; // temp array holding number of phis for a given theta
1015 :
1016 14 : NThetas = state.dataConstruction->Construct(IConst).BSDFInput.BasisMatNrows; // Note here assuming row by row input
1017 14 : Basis.NThetas = NThetas;
1018 14 : Basis.BasisMatIndex = state.dataConstruction->Construct(IConst).BSDFInput.BasisMatIndex;
1019 14 : Basis.NBasis = state.dataConstruction->Construct(IConst).BSDFInput.NBasis;
1020 14 : Basis.Grid.allocate(Basis.NBasis);
1021 14 : Thetas.allocate(NThetas + 1); // Temp array
1022 : // By convention the Thetas array contains a final point at Pi/2 which is not a basis element
1023 14 : NPhis.allocate(NThetas + 1); // Temp array
1024 14 : Basis.Thetas.allocate(NThetas + 1);
1025 14 : Basis.NPhis.allocate(NThetas + 1);
1026 :
1027 14 : Basis.Lamda.allocate(state.dataConstruction->Construct(IConst).BSDFInput.NBasis);
1028 14 : Basis.SolAng.allocate(state.dataConstruction->Construct(IConst).BSDFInput.NBasis);
1029 14 : if (state.dataConstruction->Construct(IConst).BSDFInput.BasisType == DataBSDFWindow::Basis::WINDOW) {
1030 : // WINDOW6 Basis
1031 14 : Basis.BasisType = DataBSDFWindow::Basis::WINDOW;
1032 14 : if (state.dataConstruction->Construct(IConst).BSDFInput.BasisSymmetryType == DataBSDFWindow::BasisSymmetry::None) {
1033 : // No basis symmetry
1034 14 : Basis.BasisSymmetryType = DataBSDFWindow::BasisSymmetry::None;
1035 14 : Thetas(1) = 0.0; // By convention, the first basis point is at the center (theta=0,phi=0)
1036 14 : Thetas(NThetas + 1) = 0.5 * Constant::Pi; // and there is an N+1st point (not a basis element) at Pi/2
1037 14 : NPhis(1) = 1;
1038 14 : NumElem = 1;
1039 106 : for (I = 2; I <= NThetas; ++I) {
1040 92 : Thetas(I) = state.dataConstruction->Construct(IConst).BSDFInput.BasisMat(1, I) * Constant::DegToRad;
1041 92 : NPhis(I) = std::floor(state.dataConstruction->Construct(IConst).BSDFInput.BasisMat(2, I) + 0.001);
1042 92 : if (NPhis(I) <= 0) {
1043 0 : ShowFatalError(state, "WindowComplexManager: incorrect input, no. phis must be positive.");
1044 : }
1045 92 : NumElem += NPhis(I);
1046 : }
1047 14 : MaxNPhis = maxval(NPhis({1, NThetas}));
1048 14 : Basis.Phis.allocate(NThetas + 1, MaxNPhis + 1); // N+1st Phi point (not basis element) at 2Pi
1049 14 : Basis.BasisIndex.allocate(MaxNPhis, NThetas + 1);
1050 14 : Basis.Phis = 0.0; // Initialize so undefined elements will contain zero
1051 14 : Basis.BasisIndex = 0; // Initialize so undefined elements will contain zero
1052 14 : if (NumElem != state.dataConstruction->Construct(IConst).BSDFInput.NBasis) { // Constructed Basis must match property matrices
1053 0 : ShowFatalError(state, "WindowComplexManager: Constructed basis length does not match property matrices.");
1054 : }
1055 14 : Basis.Thetas = Thetas;
1056 14 : Basis.NPhis = NPhis;
1057 14 : ElemNo = 0;
1058 120 : for (I = 1; I <= NThetas; ++I) {
1059 106 : Theta = Thetas(I);
1060 106 : if (I == 1) { // First theta value must always be zero
1061 14 : HalfDTheta = 0.5 * Thetas(I + 1);
1062 14 : LastTheta = 0.0;
1063 14 : NextTheta = Thetas(I + 1);
1064 14 : LowerTheta = 0.0;
1065 14 : UpperTheta = HalfDTheta;
1066 92 : } else if (I > 1 && I < NThetas) {
1067 78 : LastTheta = Thetas(I - 1);
1068 78 : NextTheta = Thetas(I + 1);
1069 78 : LowerTheta = UpperTheta;
1070 78 : HalfDTheta = Theta - LowerTheta;
1071 78 : UpperTheta = Theta + HalfDTheta;
1072 14 : } else if (I == NThetas) {
1073 14 : LastTheta = Thetas(I - 1);
1074 14 : NextTheta = 0.5 * Constant::Pi;
1075 14 : LowerTheta = UpperTheta; // It is assumed that Thetas(N) is the mean between the previous
1076 : // UpperTheta and pi/2.
1077 14 : UpperTheta = 0.5 * Constant::Pi;
1078 : }
1079 106 : DPhi = 2.0 * Constant::Pi / NPhis(I);
1080 106 : if (I == 1) {
1081 14 : Lamda = Constant::Pi * pow_2(std::sin(UpperTheta));
1082 14 : SolAng = 2.0 * Constant::Pi * (1.0 - std::cos(UpperTheta));
1083 : } else {
1084 92 : Lamda = 0.5 * DPhi * (pow_2(std::sin(UpperTheta)) - pow_2(std::sin(LowerTheta))); // For W6 basis, lamda is funct of Theta and
1085 : // NPhis, not individual Phi
1086 92 : SolAng = DPhi * (std::cos(LowerTheta) - std::cos(UpperTheta));
1087 : }
1088 106 : DTheta = UpperTheta - LowerTheta;
1089 106 : Basis.Phis(I, NPhis(I) + 1) = 2.0 * Constant::Pi; // Non-basis-element Phi point for table searching in Phi
1090 1616 : for (J = 1; J <= NPhis(I); ++J) {
1091 1510 : ++ElemNo;
1092 1510 : Basis.BasisIndex(J, I) = ElemNo;
1093 1510 : Phi = (J - 1) * DPhi;
1094 1510 : Basis.Phis(I, J) = Phi; // Note: this ordering of I & J are necessary to allow Phis(Theta) to
1095 : // be searched as a one-dimensional table
1096 1510 : FillBasisElement(state,
1097 : Theta,
1098 : Phi,
1099 : ElemNo,
1100 : Basis.Grid(ElemNo),
1101 : LowerTheta,
1102 : UpperTheta,
1103 : DPhi,
1104 : DataBSDFWindow::Basis::WINDOW); // This gets all the simple grid characteristics
1105 1510 : Basis.Lamda(ElemNo) = Lamda;
1106 1510 : Basis.SolAng(ElemNo) = SolAng;
1107 : }
1108 : }
1109 : } else { // BST
1110 : // Axisymmetric basis symmetry (Note this only useful specular systems, where it allows shorter data input)
1111 0 : Basis.BasisSymmetryType = DataBSDFWindow::BasisSymmetry::Axisymmetric;
1112 0 : Thetas(1) = 0.0; // By convention, the first basis point is at the center (theta=0,phi=0)
1113 0 : Thetas(NThetas + 1) = 0.5 * Constant::Pi; // and there is an N+1st point (not a basis element) at Pi/2
1114 0 : NPhis = 1; // As insurance, define one phi for each theta
1115 0 : NumElem = 1;
1116 0 : for (I = 2; I <= NThetas; ++I) {
1117 0 : Thetas(I) = state.dataConstruction->Construct(IConst).BSDFInput.BasisMat(1, I) * Constant::DegToRad;
1118 0 : ++NumElem;
1119 : }
1120 0 : Basis.Phis.allocate(1, NThetas);
1121 0 : Basis.BasisIndex.allocate(1, NThetas);
1122 0 : Basis.Phis = 0.0; // Initialize so undefined elements will contain zero
1123 0 : Basis.BasisIndex = 0; // Initialize so undefined elements will contain zero
1124 0 : if (NumElem != state.dataConstruction->Construct(IConst).BSDFInput.NBasis) { // Constructed Basis must match property matrices
1125 0 : ShowFatalError(state, "WindowComplexManager: Constructed basis length does not match property matrices.");
1126 : }
1127 0 : Basis.Thetas = Thetas;
1128 0 : Basis.NPhis = NPhis;
1129 0 : ElemNo = 0;
1130 0 : DPhi = 2.0 * Constant::Pi;
1131 0 : for (I = 1; I <= NThetas; ++I) {
1132 0 : Theta = Thetas(I);
1133 0 : if (I == 1) { // First theta value must always be zero
1134 0 : HalfDTheta = 0.5 * Thetas(I + 1);
1135 0 : LastTheta = 0.0;
1136 0 : NextTheta = Thetas(I + 1);
1137 0 : LowerTheta = 0.0;
1138 0 : UpperTheta = HalfDTheta;
1139 0 : } else if (I > 1 && I < NThetas) {
1140 0 : LastTheta = Thetas(I - 1);
1141 0 : NextTheta = Thetas(I + 1);
1142 0 : LowerTheta = UpperTheta;
1143 0 : HalfDTheta = Theta - LowerTheta;
1144 0 : UpperTheta = Theta + HalfDTheta;
1145 0 : } else if (I == NThetas) {
1146 0 : LastTheta = Thetas(I - 1);
1147 0 : NextTheta = 0.5 * Constant::Pi;
1148 0 : LowerTheta = UpperTheta; // It is assumed that Thetas(N) is the mean between the previous
1149 : // UpperTheta and pi/2.
1150 0 : UpperTheta = 0.5 * Constant::Pi;
1151 : }
1152 0 : if (I == 1) {
1153 0 : Lamda = Constant::Pi * pow_2(std::sin(UpperTheta));
1154 0 : SolAng = 2.0 * Constant::Pi * (1.0 - std::cos(UpperTheta));
1155 : } else {
1156 0 : Lamda = 0.5 * DPhi * (pow_2(std::sin(UpperTheta)) - pow_2(std::sin(LowerTheta))); // For W6 basis, lamda is funct of Theta and
1157 : // NPhis, not individual Phi
1158 0 : SolAng = DPhi * (std::cos(LowerTheta) - std::cos(UpperTheta));
1159 : }
1160 0 : DTheta = UpperTheta - LowerTheta;
1161 0 : ++ElemNo;
1162 0 : Basis.BasisIndex(1, I) = ElemNo;
1163 0 : Phi = 0.0;
1164 0 : Basis.Phis(I, 1) = Phi; // Note: this ordering of I & J are necessary to allow Phis(Theta) to
1165 : // be searched as a one-dimensional table
1166 0 : FillBasisElement(state,
1167 : Theta,
1168 : Phi,
1169 : ElemNo,
1170 : Basis.Grid(ElemNo),
1171 : LowerTheta,
1172 : UpperTheta,
1173 : DPhi,
1174 : DataBSDFWindow::Basis::WINDOW); // This gets all the simple grid characteristics
1175 0 : Basis.Lamda(ElemNo) = Lamda;
1176 0 : Basis.SolAng(ElemNo) = SolAng;
1177 : }
1178 : } // BST
1179 : } else { // BTW
1180 0 : ShowFatalError(state, "WindowComplexManager: Non-Window6 basis type not yet implemented.");
1181 : } // BTW
1182 14 : Thetas.deallocate();
1183 14 : NPhis.deallocate();
1184 14 : }
1185 :
1186 1510 : void FillBasisElement(EnergyPlusData &state,
1187 : Real64 const Theta, // Central polar angle of element
1188 : Real64 const Phi, // Central azimuthal angle of element
1189 : int const Elem, // Index number of element in basis
1190 : BasisElemDescr &BasisElem,
1191 : Real64 const LowerTheta, // Lower edge of element (polar angle)
1192 : Real64 const UpperTheta, // Upper edge of element (polar angle)
1193 : Real64 const DPhi, // Width of element (azimuthal angle)
1194 : DataBSDFWindow::Basis const InputType // Basis type
1195 : )
1196 : {
1197 :
1198 : // SUBROUTINE INFORMATION:
1199 : // AUTHOR Joe Klems
1200 : // DATE WRITTEN August 2010
1201 : // MODIFIED na
1202 : // RE-ENGINEERED na
1203 :
1204 : // PURPOSE OF THIS SUBROUTINE:
1205 : // fill in values for all the components of a basis element
1206 :
1207 1510 : if (InputType == DataBSDFWindow::Basis::WINDOW) {
1208 : // WINDOW6 Type BASIS
1209 1510 : if (Elem == 1) {
1210 : // first element, theta=0, is special case
1211 14 : BasisElem.Theta = Theta;
1212 14 : BasisElem.Phi = 0.0;
1213 14 : BasisElem.dPhi = 2.0 * Constant::Pi;
1214 14 : BasisElem.UpprTheta = UpperTheta;
1215 14 : BasisElem.dTheta = BasisElem.UpprTheta - Theta;
1216 14 : BasisElem.LwrTheta = Theta;
1217 14 : BasisElem.LwrPhi = 0.0;
1218 14 : BasisElem.UpprPhi = 2.0 * Constant::Pi;
1219 : } else {
1220 1496 : BasisElem.Theta = Theta;
1221 1496 : BasisElem.Phi = Phi;
1222 1496 : BasisElem.dPhi = DPhi;
1223 1496 : BasisElem.LwrPhi = Phi - DPhi / 2.0;
1224 1496 : BasisElem.UpprPhi = Phi + DPhi / 2.0;
1225 1496 : BasisElem.LwrTheta = LowerTheta;
1226 1496 : BasisElem.UpprTheta = UpperTheta;
1227 1496 : BasisElem.dTheta = BasisElem.UpprTheta - BasisElem.LwrTheta;
1228 : }
1229 : } else {
1230 : // Non-WINDOW6 Type Basis
1231 : // Currently not implemented
1232 0 : ShowFatalError(state, "WindowComplexManager: Custom basis type not yet implemented.");
1233 : }
1234 1510 : }
1235 :
1236 24 : void SetupComplexWindowStateGeometry(EnergyPlusData &state,
1237 : int const ISurf, // Surface number of the complex fenestration
1238 : int const IState, // State number of the complex fenestration state
1239 : int const IConst, // Pointer to construction for this state
1240 : BSDFWindowGeomDescr &Window, // Window Geometry
1241 : BSDFGeomDescr &Geom, // State Geometry
1242 : [[maybe_unused]] BSDFStateDescr &State // State Description
1243 : )
1244 : {
1245 :
1246 : // SUBROUTINE INFORMATION:
1247 : // AUTHOR J. Klems
1248 : // DATE WRITTEN June 2011
1249 : // MODIFIED na
1250 : // RE-ENGINEERED na
1251 :
1252 : // PURPOSE OF THIS SUBROUTINE:
1253 : // Define all the geometric quantites for a complex fenestration state
1254 :
1255 : using namespace Vectors;
1256 :
1257 : // Locals
1258 : // SUBROUTINE ARGUMENT DEFINITIONS:
1259 : // INTEGER, INTENT(IN) :: IWind !Complex fenestration number (in window list)
1260 :
1261 : Real64 Azimuth; // Complex fenestration azimuth
1262 : Real64 Tilt; // Complex fenestration tilt
1263 : int ElemNo; // Grid index variable
1264 : bool hit; // Surface intersection flag
1265 : int I; // Temp Indices
1266 : int J;
1267 : int IRay; // Ray index variable
1268 : int IZone; // Zone containing the complex window
1269 : int JSurf; // Secondary Surface index
1270 : int BaseSurf; // base surface index
1271 : int KBkSurf; // Back surface index
1272 : int MaxHits; // Max no of hits found
1273 : int MaxInt; // Max no of intersections found
1274 : int NSky; // No of sky rays
1275 : int NGnd; // No of gnd rays
1276 : int NReflSurf; // No of rays striking ext surfaces
1277 : int NBkSurf; // No of back surfaces
1278 : int TotHits; // Current number of surface intersections
1279 : Real64 Theta; // Basis theta angle
1280 : Real64 Phi; // Basis phi angle
1281 : Real64 HitDsq; // Squared distance to current hit pt
1282 : Real64 LeastHitDsq; // Squared distance to closest hit pt
1283 24 : Array1D<Real64> V(3); // vector array
1284 24 : Array1D_int TmpRfSfInd; // Temporary RefSurfIndex
1285 24 : Array1D_int TmpRfRyNH; // Temporary RefRayNHits
1286 24 : Array2D_int TmpHSurfNo; // Temporary HitSurfNo
1287 24 : Array2D<Real64> TmpHSurfDSq; // Temporary HitSurfDSq
1288 24 : Array1D_int TmpSkyInd; // Temporary sky index list
1289 24 : Array1D_int TmpGndInd; // Temporary gnd index list
1290 24 : Array2D_int TmpSurfInt; // Temporary index of ray intersecing back surf
1291 24 : Array2D<Real64> TmpSjdotN; // Temporary dot prod of ray angle w bk surf norm
1292 24 : Array1D_int ITemp1D; // Temporary INT 1D array
1293 24 : Array2D<Real64> Temp2D; // Temporary real 2D array
1294 : Real64 TransRSurf; // Norminal transmittance of shading surface
1295 : Real64 WtSum; // Sum for normalizing various weights
1296 : Real64 DotProd; // Temporary variable for manipulating dot product .dot.
1297 :
1298 : struct BackHitList
1299 : {
1300 : // Members
1301 : int KBkSurf; // Back surface index of the hit surface
1302 : int HitSurf; // Surface number of the hit surface
1303 : Vector HitPt; // coords of hit pt (world syst)
1304 : Real64 HitDsq; // Squared distance to the current hit pt
1305 :
1306 : // Default Constructor
1307 24 : BackHitList() : KBkSurf(0), HitSurf(0), HitDsq(0.0)
1308 : {
1309 24 : }
1310 : };
1311 :
1312 : // Object Data
1313 24 : Vector HitPt; // coords of hit pt (world syst)
1314 24 : Vector X; // position vector
1315 24 : Vector VecNorm; // outer normal vector
1316 24 : Array1D<Vector> TmpGndPt; // Temporary ground intersection list
1317 24 : Array2D<Vector> TempV2D; // Temporary vector 2D array
1318 24 : Array2D<Vector> TmpHitPt; // Temporary HitPt
1319 24 : BackHitList BSHit; // Temp list of back surface hit quantities for a ray
1320 :
1321 : // This routine primarily fills in the BSDFGeomDescr type for a given window and state
1322 : // Note that on call the incoming and outgoing basis structure types have already been filled in
1323 : // Define the central ray directions (in world coordinate system)
1324 :
1325 24 : state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState).NLayers = state.dataConstruction->Construct(IConst).BSDFInput.NumLayers;
1326 24 : Azimuth = Constant::DegToRad * state.dataSurface->Surface(ISurf).Azimuth;
1327 24 : Tilt = Constant::DegToRad * state.dataSurface->Surface(ISurf).Tilt;
1328 :
1329 : // For incoming grid
1330 :
1331 24 : Geom.sInc.allocate(Geom.Inc.NBasis);
1332 24 : Geom.sInc = Vector(0.0, 0.0, 0.0);
1333 24 : Geom.pInc.allocate(Geom.Inc.NBasis);
1334 24 : Geom.CosInc.allocate(Geom.Inc.NBasis);
1335 24 : Geom.DAInc.allocate(Geom.Inc.NBasis);
1336 24 : Geom.pInc = BSDFDaylghtPosition(0.0, 0.0);
1337 1944 : for (ElemNo = 1; ElemNo <= Geom.Inc.NBasis; ++ElemNo) {
1338 1920 : Theta = Geom.Inc.Grid(ElemNo).Theta;
1339 1920 : Phi = Geom.Inc.Grid(ElemNo).Phi;
1340 : // The following puts in the vectors depending on
1341 : // window orientation
1342 1920 : Geom.sInc(ElemNo) = WorldVectFromW6(state, Theta, Phi, RayIdentificationType::Front_Incident, Tilt, Azimuth);
1343 1920 : Geom.pInc(ElemNo) = DaylghtAltAndAzimuth(Geom.sInc(ElemNo));
1344 :
1345 1920 : Geom.CosInc(ElemNo) = std::cos(Geom.Inc.Grid(ElemNo).Theta);
1346 : // Geom%DAInc(ElemNo) = COS(Geom%pInc(ElemNo)%Altitude) * Geom%Inc%Grid(ElemNo)%dTheta * Geom%Inc%Grid(ElemNo)%dPhi
1347 : // Geom%DAInc(ElemNo) = Geom%Inc%Grid(ElemNo)%dTheta * Geom%Inc%Grid(ElemNo)%dPhi
1348 1920 : Geom.DAInc(ElemNo) = std::cos(Geom.Inc.Grid(ElemNo).Theta) * Geom.Inc.Grid(ElemNo).dTheta * Geom.Inc.Grid(ElemNo).dPhi;
1349 : }
1350 : // For outgoing grid
1351 24 : Geom.sTrn.allocate(Geom.Trn.NBasis);
1352 24 : Geom.sTrn = Vector(0.0, 0.0, 0.0);
1353 24 : Geom.pTrn.allocate(Geom.Trn.NBasis);
1354 24 : Geom.pTrn = BSDFDaylghtPosition(0.0, 0.0);
1355 1944 : for (ElemNo = 1; ElemNo <= Geom.Trn.NBasis; ++ElemNo) {
1356 1920 : Theta = Geom.Trn.Grid(ElemNo).Theta;
1357 1920 : Phi = Geom.Trn.Grid(ElemNo).Phi;
1358 : // The following puts in the vectors depending on
1359 : // window orientation
1360 1920 : Geom.sTrn(ElemNo) = WorldVectFromW6(state, Theta, Phi, RayIdentificationType::Front_Transmitted, Tilt, Azimuth);
1361 1920 : Geom.pTrn(ElemNo) = DaylghtAltAndAzimuth(Geom.sTrn(ElemNo));
1362 : }
1363 : // Incident Basis:
1364 : // Construct sky and ground ray index maps, and list of rays intersecting exterior surfaces
1365 : // Sky, and ground ray index maps, and rays that are potentially beam radiation reflected from exterior surfaces
1366 24 : TmpRfSfInd.allocate(Geom.Inc.NBasis);
1367 24 : TmpRfRyNH.allocate(Geom.Inc.NBasis);
1368 24 : TmpHSurfNo.allocate(state.dataSurface->TotSurfaces, Geom.Inc.NBasis);
1369 24 : TmpHSurfDSq.allocate(state.dataSurface->TotSurfaces, Geom.Inc.NBasis);
1370 24 : TmpHitPt.allocate(state.dataSurface->TotSurfaces, Geom.Inc.NBasis);
1371 24 : TmpSkyInd.allocate(Geom.Inc.NBasis);
1372 24 : TmpGndInd.allocate(Geom.Inc.NBasis);
1373 24 : TmpGndPt.allocate(Geom.Inc.NBasis);
1374 24 : NSky = 0;
1375 24 : NGnd = 0;
1376 24 : NReflSurf = 0;
1377 24 : TmpRfRyNH = 0;
1378 24 : Geom.NSkyUnobs = 0;
1379 24 : Geom.NGndUnobs = 0;
1380 : // Note--this loop could be repeated for different positions in the window plane (as for detailed reflection
1381 : // calculations, varying the origin in the call to PierceSurface. Essentially, have set NsubV =1.
1382 1944 : for (IRay = 1; IRay <= Geom.Inc.NBasis; ++IRay) {
1383 1920 : if (Geom.sInc(IRay).z < 0.0) {
1384 : // A ground ray
1385 816 : ++Geom.NGndUnobs;
1386 : } else {
1387 : // A sky ray
1388 1104 : ++Geom.NSkyUnobs;
1389 : }
1390 : // Exterior reveal shadowing/reflection treatment should be inserted here
1391 1920 : TotHits = 0;
1392 37070 : for (JSurf = 1; JSurf <= state.dataSurface->TotSurfaces; ++JSurf) {
1393 : // the following test will cycle on anything except exterior surfaces and shading surfaces
1394 35150 : if (state.dataSurface->Surface(JSurf).HeatTransSurf && state.dataSurface->Surface(JSurf).ExtBoundCond != ExternalEnvironment) {
1395 12652 : continue;
1396 : }
1397 : // skip the base surface containing the window and any other subsurfaces of that surface
1398 43076 : if (JSurf == state.dataSurface->Surface(ISurf).BaseSurf ||
1399 20578 : state.dataSurface->Surface(JSurf).BaseSurf == state.dataSurface->Surface(ISurf).BaseSurf) {
1400 3840 : continue;
1401 : }
1402 : // skip surfaces that face away from the window
1403 18658 : DotProd = dot(Geom.sInc(IRay), state.dataSurface->Surface(JSurf).NewellSurfaceNormalVector);
1404 18658 : if (DotProd >= 0.0) {
1405 8862 : continue;
1406 : }
1407 9796 : hit = PierceSurface(state, JSurf, state.dataSurface->Surface(ISurf).Centroid, Geom.sInc(IRay), HitPt);
1408 9796 : if (!hit) {
1409 9796 : continue; // Miss: Try next surface
1410 : }
1411 0 : if (TotHits == 0) {
1412 : // First hit for this ray
1413 0 : TotHits = 1;
1414 0 : ++NReflSurf;
1415 0 : TmpRfSfInd(NReflSurf) = IRay;
1416 0 : TmpRfRyNH(NReflSurf) = 1;
1417 0 : TmpHSurfNo(1, NReflSurf) = JSurf;
1418 0 : TmpHitPt(1, NReflSurf) = HitPt;
1419 0 : V = HitPt - state.dataSurface->Surface(ISurf).Centroid; // vector array from window ctr to hit pt
1420 0 : LeastHitDsq = magnitude_squared(V); // dist^2 window ctr to hit pt
1421 0 : TmpHSurfDSq(1, NReflSurf) = LeastHitDsq;
1422 0 : if (!state.dataSurface->Surface(JSurf).HeatTransSurf && state.dataSurface->Surface(JSurf).shadowSurfSched != nullptr) {
1423 0 : TransRSurf = 1.0; // If a shadowing surface may have a scheduled transmittance,
1424 : // treat it here as completely transparent
1425 : } else {
1426 0 : TransRSurf = 0.0;
1427 : }
1428 : } else {
1429 0 : V = HitPt - state.dataSurface->Surface(ISurf).Centroid;
1430 0 : HitDsq = magnitude_squared(V);
1431 0 : if (HitDsq >= LeastHitDsq) {
1432 0 : if (TransRSurf > 0.0) { // forget the new hit if the closer hit is opaque
1433 0 : J = TotHits + 1;
1434 0 : if (TotHits > 1) {
1435 0 : for (I = 2; I <= TotHits; ++I) {
1436 0 : if (HitDsq < TmpHSurfDSq(I, NReflSurf)) {
1437 0 : J = I;
1438 0 : break;
1439 : }
1440 : }
1441 0 : if (!state.dataSurface->Surface(JSurf).HeatTransSurf &&
1442 0 : state.dataSurface->Surface(JSurf).shadowSurfSched == nullptr) {
1443 : // The new hit is opaque, so we can drop all the hits further away
1444 0 : TmpHSurfNo(J, NReflSurf) = JSurf;
1445 0 : TmpHitPt(J, NReflSurf) = HitPt;
1446 0 : TmpHSurfDSq(J, NReflSurf) = HitDsq;
1447 0 : TotHits = J;
1448 : } else {
1449 : // The new hit is scheduled (presumed transparent), so keep the more distant hits
1450 : // Note that all the hists in the list will be transparent except the last,
1451 : // which may be either transparent or opaque
1452 0 : if (TotHits >= J) {
1453 0 : for (I = TotHits; I >= J; --I) {
1454 0 : TmpHSurfNo(I + 1, NReflSurf) = TmpHSurfNo(I, NReflSurf);
1455 0 : TmpHitPt(I + 1, NReflSurf) = TmpHitPt(I, NReflSurf);
1456 0 : TmpHSurfDSq(I + 1, NReflSurf) = TmpHSurfDSq(I, NReflSurf);
1457 : }
1458 0 : TmpHSurfNo(J, NReflSurf) = JSurf;
1459 0 : TmpHitPt(J, NReflSurf) = HitPt;
1460 0 : TmpHSurfDSq(J, NReflSurf) = HitDsq;
1461 0 : ++TotHits;
1462 : }
1463 : }
1464 : }
1465 : }
1466 : } else {
1467 : // A new closest hit. If it is opaque, drop the current hit list,
1468 : // otherwise add it at the front
1469 0 : LeastHitDsq = HitDsq;
1470 0 : if (!state.dataSurface->Surface(JSurf).HeatTransSurf && state.dataSurface->Surface(JSurf).shadowSurfSched != nullptr) {
1471 0 : TransRSurf = 1.0; // New closest hit is transparent, keep the existing hit list
1472 0 : for (I = TotHits; I >= 1; --I) {
1473 0 : TmpHSurfNo(I + 1, NReflSurf) = TmpHSurfNo(I, NReflSurf);
1474 0 : TmpHitPt(I + 1, NReflSurf) = TmpHitPt(I, NReflSurf);
1475 0 : TmpHSurfDSq(I + 1, NReflSurf) = TmpHSurfDSq(I, NReflSurf);
1476 0 : ++TotHits;
1477 : }
1478 : } else {
1479 0 : TransRSurf = 0.0; // New closest hit is opaque, drop the existing hit list
1480 0 : TotHits = 1;
1481 : }
1482 0 : TmpHSurfNo(1, NReflSurf) = JSurf; // In either case the new hit is put in position 1
1483 0 : TmpHitPt(1, NReflSurf) = HitPt;
1484 0 : TmpHSurfDSq(1, NReflSurf) = LeastHitDsq;
1485 : }
1486 : }
1487 : } // End of loop over surfaces
1488 1920 : if (TotHits <= 0) {
1489 : // This ray reached the sky or ground unobstructed
1490 1920 : if (Geom.sInc(IRay).z < 0.0) {
1491 : // A ground ray
1492 816 : ++NGnd;
1493 816 : TmpGndInd(NGnd) = IRay;
1494 816 : TmpGndPt(NGnd).x = state.dataSurface->Surface(ISurf).Centroid.x -
1495 816 : (Geom.sInc(IRay).x / Geom.sInc(IRay).z) * state.dataSurface->Surface(ISurf).Centroid.z;
1496 816 : TmpGndPt(NGnd).y = state.dataSurface->Surface(ISurf).Centroid.y -
1497 816 : (Geom.sInc(IRay).y / Geom.sInc(IRay).z) * state.dataSurface->Surface(ISurf).Centroid.z;
1498 816 : TmpGndPt(NGnd).z = 0.0;
1499 : } else {
1500 : // A sky ray
1501 1104 : ++NSky;
1502 1104 : TmpSkyInd(NSky) = IRay;
1503 : }
1504 : } else {
1505 : // Save the number of hits for this ray
1506 0 : TmpRfRyNH(NReflSurf) = TotHits;
1507 : }
1508 : } // End of loop over basis rays
1509 : // Store results of indexing the incident basis for this window
1510 24 : Geom.NSky = NSky;
1511 24 : Geom.NGnd = NGnd;
1512 24 : Geom.NReflSurf = NReflSurf;
1513 24 : Geom.SkyIndex.allocate(NSky);
1514 24 : Geom.SkyIndex = TmpSkyInd({1, NSky});
1515 24 : TmpSkyInd.deallocate();
1516 24 : Geom.GndIndex.allocate(NGnd);
1517 24 : Geom.GndPt.allocate(NGnd);
1518 24 : Geom.GndIndex = TmpGndInd({1, NGnd});
1519 24 : Geom.GndPt = TmpGndPt({1, NGnd});
1520 24 : TmpGndInd.deallocate();
1521 24 : TmpGndPt.deallocate();
1522 24 : MaxHits = maxval(TmpRfRyNH);
1523 24 : Geom.RefSurfIndex.allocate(NReflSurf);
1524 24 : Geom.RefRayNHits.allocate(NReflSurf);
1525 24 : Geom.HitSurfNo.allocate(MaxHits, NReflSurf);
1526 24 : Geom.HitSurfDSq.allocate(MaxHits, NReflSurf);
1527 24 : Geom.HitPt.allocate(MaxHits, NReflSurf);
1528 24 : Geom.RefSurfIndex = TmpRfSfInd({1, NReflSurf});
1529 24 : Geom.RefRayNHits = TmpRfRyNH({1, NReflSurf});
1530 24 : Geom.HitSurfNo = 0;
1531 24 : Geom.HitSurfDSq = 0.0;
1532 24 : Geom.HitPt = Vector(0.0, 0.0, 0.0);
1533 24 : for (I = 1; I <= NReflSurf; ++I) {
1534 0 : TotHits = TmpRfRyNH(I);
1535 0 : Geom.HitSurfNo({1, TotHits}, I) = TmpHSurfNo({1, TotHits}, I);
1536 0 : Geom.HitSurfDSq({1, TotHits}, I) = TmpHSurfDSq({1, TotHits}, I);
1537 0 : Geom.HitPt({1, TotHits}, I) = TmpHitPt({1, TotHits}, I);
1538 : }
1539 24 : TmpRfRyNH.deallocate();
1540 24 : TmpRfSfInd.deallocate();
1541 24 : TmpHSurfNo.deallocate();
1542 24 : TmpHSurfDSq.deallocate();
1543 24 : TmpHitPt.deallocate();
1544 : // In above scheme sky and ground rays are those that intesect no exterior surfaces.
1545 : // The list of hit points is compiled for later (future?) calculation
1546 : // of reflections from these surfaces. The hit list for each ray includes all
1547 : // surfaces with schedulable transmittance intersected by the ray,
1548 : // in order of increasing distance, up to the first opaque surface.
1549 : // Rays that intesect one or more schedulable transmittance but no opaque
1550 : // surfaces (therefore may reach the sky or ground) are left out of the sky/ground
1551 : // calcuation. A correction for these rays could/should be made after the
1552 : // shading calculation.
1553 : // Now calculate weights for averaging the transmittance matrix
1554 : // Sky Weights
1555 24 : Geom.SolSkyWt.allocate(NSky);
1556 1128 : for (I = 1; I <= NSky; ++I) {
1557 1104 : J = Geom.SkyIndex(I);
1558 1104 : Geom.SolSkyWt(I) = SkyWeight(Geom.sInc(J));
1559 : }
1560 24 : WtSum = sum(Geom.SolSkyWt({1, NSky}));
1561 24 : if (WtSum > Constant::rTinyValue) {
1562 24 : Geom.SolSkyWt({1, NSky}) /= WtSum;
1563 : } else {
1564 0 : Geom.SolSkyWt({1, NSky}) = 0.0;
1565 : }
1566 : // SkyGround Weights
1567 24 : Geom.SolSkyGndWt.allocate(NGnd);
1568 840 : for (I = 1; I <= NGnd; ++I) {
1569 816 : Geom.SolSkyGndWt(I) = SkyGndWeight(Geom.GndPt(I));
1570 : }
1571 24 : WtSum = sum(Geom.SolSkyGndWt({1, NGnd}));
1572 24 : if (WtSum > Constant::rTinyValue) {
1573 24 : Geom.SolSkyGndWt({1, NGnd}) /= WtSum;
1574 : } else {
1575 0 : Geom.SolSkyGndWt({1, NGnd}) = 0.0;
1576 : }
1577 : // Weights for beam reflected from ground are calculated after shading
1578 : // interval is determined
1579 : // Transmitted Basis:
1580 : // Construct back surface intersection maps
1581 24 : IZone = state.dataSurface->Surface(ISurf).Zone;
1582 24 : NBkSurf = Window.NBkSurf;
1583 24 : Geom.NSurfInt.allocate(NBkSurf);
1584 24 : Geom.NSurfInt = 0; // Initialize the number of intersections to zero
1585 24 : TmpSurfInt.allocate(Geom.Trn.NBasis, NBkSurf);
1586 24 : TmpSjdotN.allocate(Geom.Trn.NBasis, NBkSurf);
1587 : // Find the intersections of the basis rays with the back surfaces
1588 1944 : for (IRay = 1; IRay <= Geom.Trn.NBasis; ++IRay) { // ray loop
1589 1920 : TotHits = 0;
1590 : // Insert treatment of intersection & reflection from interior reveals here
1591 14244 : for (KBkSurf = 1; KBkSurf <= NBkSurf; ++KBkSurf) { // back surf loop
1592 12324 : BaseSurf = state.dataSurface->Surface(ISurf).BaseSurf; // ShadowComb is organized by base surface
1593 12324 : JSurf = state.dataShadowComb->ShadowComb(BaseSurf).BackSurf(KBkSurf); // these are all proper back surfaces
1594 12324 : hit = PierceSurface(state, JSurf, state.dataSurface->Surface(ISurf).Centroid, Geom.sTrn(IRay), HitPt);
1595 12324 : if (!hit) {
1596 10132 : continue; // Miss: Try next surface
1597 : }
1598 2192 : if (TotHits == 0) {
1599 : // First hit for this ray
1600 1920 : TotHits = 1;
1601 1920 : BSHit.KBkSurf = KBkSurf;
1602 1920 : BSHit.HitSurf = JSurf;
1603 1920 : BSHit.HitPt = HitPt;
1604 1920 : V = HitPt - state.dataSurface->Surface(ISurf).Centroid;
1605 1920 : BSHit.HitDsq = magnitude_squared(V);
1606 272 : } else if (BSHit.HitSurf == state.dataSurface->Surface(JSurf).BaseSurf) {
1607 : // another hit, check whether this is a subsurface of a previously hit base surface
1608 : // (which would be listed first in the Surface array)
1609 : // if so, replace the previous hit with this one
1610 272 : ++TotHits;
1611 272 : BSHit.KBkSurf = KBkSurf;
1612 272 : BSHit.HitSurf = JSurf;
1613 272 : BSHit.HitPt = HitPt;
1614 272 : V = HitPt - state.dataSurface->Surface(ISurf).Centroid;
1615 272 : BSHit.HitDsq = magnitude_squared(V);
1616 : } else {
1617 0 : ++TotHits;
1618 : // is the new hit closer than the previous one (i.e., zone not strictly convex)?
1619 : // if so, take the closer hit
1620 0 : V = HitPt - state.dataSurface->Surface(ISurf).Centroid;
1621 0 : HitDsq = magnitude_squared(V);
1622 0 : if (HitDsq < BSHit.HitDsq) {
1623 0 : BSHit.KBkSurf = KBkSurf;
1624 0 : BSHit.HitSurf = JSurf;
1625 0 : BSHit.HitPt = HitPt;
1626 0 : BSHit.HitDsq = HitDsq;
1627 : }
1628 : }
1629 : } // back surf loop
1630 1920 : if (TotHits == 0) {
1631 : // this should not happen--means a ray has gotten lost
1632 : // ShowWarningError(state, "BSDF--Zone surfaces do not completely enclose zone--transmitted ray lost");
1633 : } else {
1634 1920 : KBkSurf = BSHit.KBkSurf;
1635 1920 : JSurf = BSHit.HitSurf;
1636 1920 : ++Geom.NSurfInt(KBkSurf);
1637 1920 : TmpSurfInt(Geom.NSurfInt(KBkSurf), KBkSurf) = IRay;
1638 1920 : VecNorm = state.dataSurface->Surface(JSurf).OutNormVec;
1639 1920 : TmpSjdotN(Geom.NSurfInt(KBkSurf), KBkSurf) = dot(Geom.sTrn(IRay), VecNorm);
1640 : }
1641 : } // ray loop
1642 : // All rays traced, now put away the results in the temporary arrays
1643 24 : MaxInt = maxval(Geom.NSurfInt);
1644 24 : Geom.SurfInt.allocate(MaxInt, Window.NBkSurf);
1645 24 : Geom.SjdotN.allocate(MaxInt, Window.NBkSurf);
1646 24 : Geom.SurfInt = 0;
1647 180 : for (I = 1; I <= Window.NBkSurf; ++I) {
1648 156 : Geom.SurfInt({1, Geom.NSurfInt(I)}, I) = TmpSurfInt({1, Geom.NSurfInt(I)}, I);
1649 156 : Geom.SjdotN({1, Geom.NSurfInt(I)}, I) = TmpSjdotN({1, Geom.NSurfInt(I)}, I);
1650 : }
1651 :
1652 24 : TmpSurfInt.deallocate();
1653 24 : TmpSjdotN.deallocate();
1654 24 : }
1655 :
1656 24 : void CalcWindowStaticProperties(EnergyPlusData &state,
1657 : int const ISurf, // Surface number of the complex fenestration
1658 : int const IState, // State number of the complex fenestration state
1659 : BSDFWindowGeomDescr &Window, // Window Geometry
1660 : BSDFGeomDescr &Geom, // State Geometry
1661 : BSDFStateDescr &State // State Description
1662 : )
1663 : {
1664 :
1665 : // SUBROUTINE INFORMATION:
1666 : // AUTHOR Joe Klems
1667 : // DATE WRITTEN <date_written>
1668 : // MODIFIED na
1669 : // RE-ENGINEERED na
1670 :
1671 : // PURPOSE OF THIS SUBROUTINE:
1672 : // Calculates those optical properties of all the Complex Fenestrations that
1673 : // do not depend on the beam direction (hence, on hour and time step)
1674 :
1675 : using namespace Vectors;
1676 :
1677 : int IConst; // Pointer to construction for this fenestration
1678 24 : int I(0); // general purpose index
1679 24 : int J(0); // general purpose index
1680 24 : int JJ(0); // general purpose index--ray
1681 24 : int L(0); // general purpose index--layer
1682 24 : int M(0); // general purpose index--ray
1683 : int KBkSurf; // back surface index
1684 : int JSurf; // surface number (used for back surface)
1685 : int BaseSurf; // base surface number (used for finding back surface)
1686 : Real64 Sum1; // general purpose temporary sum
1687 : Real64 Sum2; // general purpose temporary sum
1688 : Real64 Sum3; // general purpose temporary sum
1689 : Real64 Hold; // temp variable
1690 :
1691 24 : IConst = state.dataSurface->SurfaceWindow(ISurf).ComplexFen.State(IState).Konst;
1692 :
1693 : // Calculate the hemispherical-hemispherical transmittance
1694 :
1695 24 : Sum1 = 0.0;
1696 24 : Sum2 = 0.0;
1697 1944 : for (J = 1; J <= Geom.Inc.NBasis; ++J) { // Incident ray loop
1698 1920 : Sum2 += Geom.Inc.Lamda(J);
1699 216360 : for (M = 1; M <= Geom.Trn.NBasis; ++M) { // Outgoing ray loop
1700 214440 : Sum1 += Geom.Inc.Lamda(J) * Geom.Trn.Lamda(M) * state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(M, J);
1701 : } // Outgoing ray loop
1702 : } // Incident ray loop
1703 24 : if (Sum2 > 0) {
1704 24 : State.WinDiffTrans = Sum1 / Sum2;
1705 : } else {
1706 0 : State.WinDiffTrans = 0.0;
1707 0 : ShowWarningError(state, "BSDF--Inc basis has zero projected solid angle");
1708 : }
1709 :
1710 : // Calculate the hemispherical-hemispherical transmittance for visible spetrum
1711 :
1712 24 : Sum1 = 0.0;
1713 24 : Sum2 = 0.0;
1714 1944 : for (J = 1; J <= Geom.Inc.NBasis; ++J) { // Incident ray loop
1715 1920 : Sum2 += Geom.Inc.Lamda(J);
1716 216360 : for (M = 1; M <= Geom.Trn.NBasis; ++M) { // Outgoing ray loop
1717 214440 : Sum1 += Geom.Inc.Lamda(J) * Geom.Trn.Lamda(M) * state.dataConstruction->Construct(IConst).BSDFInput.VisFrtTrans(M, J);
1718 : } // Outgoing ray loop
1719 : } // Incident ray loop
1720 24 : if (Sum2 > 0.0) {
1721 24 : State.WinDiffVisTrans = Sum1 / Sum2;
1722 : } else {
1723 0 : State.WinDiffVisTrans = 0.0;
1724 0 : ShowWarningError(state, "BSDF--Inc basis has zero projected solid angle");
1725 : }
1726 :
1727 : // Set the nominal diffuse transmittance so the surface isn't mistaken as opaque
1728 : // Simon: Commented this out. We are not using TransDiff and it is already set to 0.1 in input routines.
1729 : // Construct( IConst ).TransDiff = SurfaceWindow( ISurf ).ComplexFen.State( IState ).WinDiffTrans;
1730 : // Calculate Window Sky Transmittance (transmitted radiation assumed diffuse)
1731 : // and Sky Absorptance (by layer)
1732 24 : Sum1 = 0.0;
1733 24 : Sum2 = 0.0;
1734 24 : Sum3 = 0.0;
1735 1128 : for (JJ = 1; JJ <= Geom.NSky; ++JJ) {
1736 122184 : for (M = 1; M <= Geom.Trn.NBasis; ++M) {
1737 121080 : J = Geom.SkyIndex(JJ);
1738 121080 : Sum1 +=
1739 121080 : Geom.SolSkyWt(JJ) * state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(M, J) * Geom.Inc.Lamda(J) * Geom.Trn.Lamda(M);
1740 : }
1741 : }
1742 1128 : for (JJ = 1; JJ <= Geom.NSky; ++JJ) {
1743 1104 : J = Geom.SkyIndex(JJ);
1744 1104 : Sum2 += Geom.SolSkyWt(JJ) * Geom.Inc.Lamda(J);
1745 : }
1746 :
1747 24 : if (Sum2 != 0.0) {
1748 24 : State.WinSkyTrans = Sum1 / Sum2;
1749 : } else {
1750 0 : State.WinSkyTrans = 0.0;
1751 : }
1752 :
1753 24 : State.WinSkyFtAbs.allocate(State.NLayers);
1754 : // Also allocate the beam quantities for this state
1755 92 : for (L = 1; L <= State.NLayers; ++L) {
1756 68 : Sum3 = 0.0;
1757 3224 : for (JJ = 1; JJ <= Geom.NSky; ++JJ) {
1758 3156 : J = Geom.SkyIndex(JJ);
1759 3156 : Sum3 += Geom.SolSkyWt(JJ) * Geom.Inc.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).FrtAbs(J, 1);
1760 : }
1761 :
1762 68 : if (Sum2 != 0.0) {
1763 68 : State.WinSkyFtAbs(L) = Sum3 / Sum2;
1764 : } else {
1765 0 : State.WinSkyFtAbs(L) = 0.0;
1766 : }
1767 : }
1768 :
1769 : // Calculate Window Sky/Ground Transmittance
1770 : //(applies to ground-reflected sky radiation, transmitted radiation assumed diffuse)
1771 : // This is the same calculation as the sky transmittance, except that the set of incident
1772 : // rays and the ray weights are different
1773 : // Also calculate Window Sky/Ground Absorptance (by layer)
1774 24 : Sum1 = 0.0;
1775 24 : Sum2 = 0.0;
1776 24 : Sum3 = 0.0;
1777 :
1778 840 : for (JJ = 1; JJ <= Geom.NGnd; ++JJ) {
1779 94176 : for (M = 1; M <= Geom.Trn.NBasis; ++M) {
1780 93360 : J = Geom.GndIndex(JJ);
1781 93360 : Sum1 += Geom.SolSkyGndWt(JJ) * state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(M, J) * Geom.Inc.Lamda(J) *
1782 93360 : Geom.Trn.Lamda(M);
1783 : }
1784 : }
1785 :
1786 840 : for (JJ = 1; JJ <= Geom.NGnd; ++JJ) {
1787 816 : J = Geom.GndIndex(JJ);
1788 816 : Sum2 += Geom.SolSkyGndWt(JJ) * Geom.Inc.Lamda(J);
1789 : }
1790 :
1791 24 : if (Sum2 != 0.0) {
1792 24 : State.WinSkyGndTrans = Sum1 / Sum2;
1793 : } else {
1794 0 : State.WinSkyGndTrans = 0.0;
1795 : }
1796 :
1797 24 : State.WinSkyGndAbs.allocate(State.NLayers);
1798 92 : for (L = 1; L <= State.NLayers; ++L) {
1799 68 : Sum3 = 0.0;
1800 2404 : for (JJ = 1; JJ <= Geom.NGnd; ++JJ) {
1801 2336 : J = Geom.GndIndex(JJ);
1802 2336 : Sum3 += Geom.SolSkyGndWt(JJ) * Geom.Inc.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).FrtAbs(J, 1);
1803 : }
1804 :
1805 68 : if (Sum2 != 0.0) {
1806 68 : State.WinSkyGndAbs(L) = Sum3 / Sum2;
1807 : } else {
1808 0 : State.WinSkyGndAbs(L) = 0.0;
1809 : }
1810 : }
1811 :
1812 : // Calculate Window Back Hemispherical Reflectance and Layer Back Hemispherical Absorptance
1813 24 : Sum1 = 0.0;
1814 24 : Sum2 = 0.0;
1815 24 : Sum3 = 0.0;
1816 : // Note this again assumes the equivalence Inc basis = transmission basis for back incidence and
1817 : // Trn basis = incident basis for back incidence
1818 1944 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
1819 216360 : for (M = 1; M <= Geom.Inc.NBasis; ++M) {
1820 214440 : Sum1 += state.dataConstruction->Construct(IConst).BSDFInput.SolBkRefl(M, J) * Geom.Trn.Lamda(J) * Geom.Inc.Lamda(M);
1821 : }
1822 : }
1823 1944 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
1824 1920 : Sum2 += Geom.Trn.Lamda(J);
1825 : }
1826 :
1827 24 : if (Sum2 != 0.0) {
1828 24 : State.WinBkHemRefl = Sum1 / Sum2;
1829 : } else {
1830 0 : State.WinBkHemRefl = 0.0;
1831 : }
1832 :
1833 24 : state.dataConstruction->Construct(IConst).ReflectSolDiffBack = State.WinBkHemRefl;
1834 :
1835 24 : State.WinBkHemAbs.allocate(State.NLayers);
1836 92 : for (L = 1; L <= State.NLayers; ++L) {
1837 5560 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
1838 5492 : Sum3 += Geom.Trn.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).BkAbs(J, 1);
1839 : }
1840 :
1841 68 : if (Sum2 != 0.0) {
1842 68 : State.WinBkHemAbs(L) = Sum3 / Sum2;
1843 : } else {
1844 0 : State.WinBkHemAbs(L) = 0.0;
1845 : }
1846 :
1847 : // Put this into the construction for use in non-detailed optical calculations
1848 68 : state.dataConstruction->Construct(IConst).AbsDiffBack(L) = State.WinBkHemAbs(L);
1849 : }
1850 :
1851 : // Calculate Window Layer Front Hemispherical Absorptance
1852 24 : Sum1 = 0.0;
1853 24 : Sum2 = 0.0;
1854 1944 : for (J = 1; J <= Geom.Inc.NBasis; ++J) {
1855 1920 : Sum2 += Geom.Inc.Lamda(J);
1856 : }
1857 24 : State.WinFtHemAbs.allocate(State.NLayers);
1858 92 : for (L = 1; L <= State.NLayers; ++L) {
1859 68 : Sum1 = 0.0;
1860 5560 : for (J = 1; J <= Geom.Inc.NBasis; ++J) {
1861 5492 : Sum1 += Geom.Inc.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).FrtAbs(J, 1);
1862 : }
1863 :
1864 68 : if (Sum2 != 0.0) {
1865 68 : State.WinFtHemAbs(L) = Sum1 / Sum2;
1866 : } else {
1867 0 : State.WinFtHemAbs(L) = 0.0;
1868 : }
1869 :
1870 : // Put this into the construction for use in non-detailed optical calculations
1871 68 : state.dataConstruction->Construct(IConst).AbsDiff(L) = State.WinFtHemAbs(L);
1872 : }
1873 :
1874 : // Calculate Window Back Hemispherical Visible Reflectance
1875 24 : Sum1 = 0.0;
1876 24 : Sum2 = 0.0;
1877 : // Note this again assumes the equivalence Inc basis = transmission basis for back incidence and
1878 : // Trn basis = incident basis for back incidence
1879 1944 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
1880 216360 : for (M = 1; M <= Geom.Inc.NBasis; ++M) {
1881 214440 : Sum1 += state.dataConstruction->Construct(IConst).BSDFInput.VisBkRefl(M, J) * Geom.Trn.Lamda(J) * Geom.Inc.Lamda(M);
1882 : }
1883 : }
1884 1944 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
1885 1920 : Sum2 += Geom.Trn.Lamda(J);
1886 : }
1887 :
1888 24 : if (Sum2 != 0.0) {
1889 24 : State.WinBkHemVisRefl = Sum1 / Sum2;
1890 : } else {
1891 0 : State.WinBkHemVisRefl = 0.0;
1892 : }
1893 :
1894 24 : state.dataConstruction->Construct(IConst).ReflectVisDiffBack = State.WinBkHemVisRefl;
1895 :
1896 : // * * * *
1897 : // Note potential problem if one relaxes the assumption that Inc and Trn basis have same structure:
1898 : // The following calculations are made for the set of ray numbers defined in the Trn basis that
1899 : // were determined to connect the center of the window to a particular back surface.
1900 : // Here it is assumed that one can reverse these rays and get an equivalent set in the Trn
1901 : // basis for back-incidence quantities: back transmittance and back layer absorptance
1902 : // This assumption may fail if the Inc and Trn bases are allowed to have different structure.
1903 : // Note also that in this case one would need to rethink the relationship of the basis
1904 : // definitions to back-incidence quantities: possibly this would
1905 : // also require that the basis for back incident quantities be
1906 : // different from the Trn basis, and similarly the basis for backward outgoing rays
1907 : // be different from the Inc basis.
1908 :
1909 : // * * * *
1910 : // Note that we are assuming that for back incidence the layer numberings are the same
1911 : // as for front incidence, i.e., from outside to inside when incidence is from inside
1912 : // * * * *
1913 : // For back surfaces that are complex fenestrations, calculate the directional-hemispherical back
1914 : // reflectance and the directional back absorptance by layer for this fenestration receiving
1915 : // radiation via the back surface
1916 : // Make this calculation only for cases where the back surface is a Complex Fenestration
1917 : // First allocate the back surface section of the state properties
1918 24 : if (!allocated(State.BkSurf)) {
1919 0 : State.BkSurf.allocate(Window.NBkSurf);
1920 : }
1921 180 : for (KBkSurf = 1; KBkSurf <= Window.NBkSurf; ++KBkSurf) { // back surface loop
1922 156 : BaseSurf = state.dataSurface->Surface(ISurf).BaseSurf; // ShadowComb is organized by base surface
1923 156 : JSurf = state.dataShadowComb->ShadowComb(BaseSurf).BackSurf(KBkSurf);
1924 156 : if (state.dataSurface->SurfWinWindowModelType(JSurf) != WindowModel::BSDF) {
1925 120 : continue;
1926 : }
1927 :
1928 : // Directional-hemispherical back reflectance
1929 36 : Sum1 = 0.0;
1930 36 : Sum2 = 0.0;
1931 308 : for (J = 1; J <= Geom.NSurfInt(KBkSurf); ++J) { // Inc Ray loop
1932 272 : Sum2 += Geom.Trn.Lamda(Geom.SurfInt(J, KBkSurf));
1933 33888 : for (M = 1; M <= Geom.Inc.NBasis; ++M) { // Outgoing Ray loop
1934 33616 : Sum1 += Geom.Trn.Lamda(Geom.SurfInt(J, KBkSurf)) * Geom.Inc.Lamda(M) *
1935 33616 : state.dataConstruction->Construct(IConst).BSDFInput.SolBkRefl(Geom.SurfInt(J, KBkSurf), M);
1936 : } // Outgoing Ray loop
1937 : } // Inc Ray loop
1938 36 : if (Sum2 > 0.0) {
1939 36 : Hold = Sum1 / Sum2;
1940 900 : for (I = 1; I <= 24; ++I) {
1941 4896 : for (J = 1; J <= state.dataGlobal->TimeStepsInHour; ++J) {
1942 4032 : State.BkSurf(KBkSurf).WinDHBkRefl(I, J) = Hold;
1943 : }
1944 : }
1945 : } else {
1946 0 : for (I = 1; I <= 24; ++I) {
1947 0 : for (J = 1; J <= state.dataGlobal->TimeStepsInHour; ++J) {
1948 0 : State.BkSurf(KBkSurf).WinDHBkRefl(I, J) = 0.0;
1949 : }
1950 : }
1951 : }
1952 :
1953 : // Directional layer back absorption
1954 147 : for (L = 1; L <= State.NLayers; ++L) { // layer loop
1955 111 : Sum1 = 0.0;
1956 111 : Sum2 = 0.0;
1957 962 : for (J = 1; J <= Geom.NSurfInt(KBkSurf); ++J) { // Inc Ray loop
1958 851 : Sum2 += Geom.Trn.Lamda(Geom.SurfInt(J, KBkSurf));
1959 851 : Sum1 += Geom.Trn.Lamda(Geom.SurfInt(J, KBkSurf)) *
1960 851 : state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).BkAbs(Geom.SurfInt(J, KBkSurf), 1);
1961 : } // Inc Ray loop
1962 111 : if (Sum2 > 0.0) {
1963 111 : Hold = Sum1 / Sum2;
1964 2775 : for (I = 1; I <= 24; ++I) {
1965 15192 : for (J = 1; J <= state.dataGlobal->TimeStepsInHour; ++J) {
1966 12528 : State.BkSurf(KBkSurf).WinDirBkAbs(I, J, L) = Hold;
1967 : }
1968 : }
1969 : } else {
1970 0 : for (I = 1; I <= 24; ++I) {
1971 0 : for (J = 1; J <= state.dataGlobal->TimeStepsInHour; ++J) {
1972 0 : State.BkSurf(KBkSurf).WinDirBkAbs(I, J, L) = 0.0;
1973 : }
1974 : }
1975 : }
1976 :
1977 : } // layer loop
1978 : } // back surface loop
1979 :
1980 : // ********************************************************************************
1981 : // Allocation and calculation of integrated values for front of window surface
1982 : // ********************************************************************************
1983 :
1984 : // Sum of front absorptances for each incident direction (integration of absorptances)
1985 24 : if (!allocated(State.IntegratedFtAbs)) {
1986 24 : State.IntegratedFtAbs.allocate(Geom.Inc.NBasis);
1987 : }
1988 1944 : for (J = 1; J <= Geom.Inc.NBasis; ++J) {
1989 1920 : Sum1 = 0.0;
1990 7412 : for (L = 1; L <= State.NLayers; ++L) { // layer loop
1991 5492 : Sum1 += state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).FrtAbs(J, 1);
1992 : }
1993 1920 : State.IntegratedFtAbs(J) = Sum1;
1994 : }
1995 :
1996 : // Integrating front transmittance
1997 24 : if (!allocated(State.IntegratedFtTrans)) {
1998 24 : State.IntegratedFtTrans.allocate(Geom.Inc.NBasis);
1999 : }
2000 1944 : for (J = 1; J <= Geom.Inc.NBasis; ++J) { // Incident ray loop
2001 1920 : Sum1 = 0.0;
2002 216360 : for (M = 1; M <= Geom.Trn.NBasis; ++M) { // Outgoing ray loop
2003 214440 : Sum1 += Geom.Trn.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.SolFrtTrans(J, M);
2004 : } // Outgoing ray loop
2005 1920 : State.IntegratedFtTrans(J) = Sum1;
2006 : } // Incident ray loop
2007 :
2008 24 : if (!allocated(State.IntegratedFtRefl)) {
2009 24 : State.IntegratedFtRefl.allocate(Geom.Inc.NBasis);
2010 : }
2011 : // Integrating front reflectance
2012 1944 : for (J = 1; J <= Geom.Inc.NBasis; ++J) { // Incoming ray loop
2013 1920 : State.IntegratedFtRefl(J) = 1 - State.IntegratedFtTrans(J) - State.IntegratedFtAbs(J);
2014 : } // Incoming ray loop
2015 :
2016 : // ********************************************************************************
2017 : // Allocation and calculation of integrated values for back of window surface
2018 : // ********************************************************************************
2019 :
2020 : // Sum of back absorptances for each incident direction (integration of absorptances)
2021 24 : if (!allocated(State.IntegratedBkAbs)) {
2022 24 : State.IntegratedBkAbs.allocate(Geom.Trn.NBasis);
2023 : }
2024 1944 : for (J = 1; J <= Geom.Trn.NBasis; ++J) {
2025 1920 : Sum1 = 0.0;
2026 7412 : for (L = 1; L <= State.NLayers; ++L) { // layer loop
2027 5492 : Sum1 += state.dataConstruction->Construct(IConst).BSDFInput.Layer(L).BkAbs(J, 1);
2028 : }
2029 1920 : State.IntegratedBkAbs(J) = Sum1;
2030 : }
2031 :
2032 : // Integrating back reflectance
2033 24 : if (!allocated(State.IntegratedBkRefl)) {
2034 24 : State.IntegratedBkRefl.allocate(Geom.Trn.NBasis);
2035 : }
2036 1944 : for (J = 1; J <= Geom.Trn.NBasis; ++J) { // Outgoing ray loop
2037 1920 : Sum1 = 0.0;
2038 216360 : for (M = 1; M <= Geom.Inc.NBasis; ++M) { // Incident ray loop
2039 214440 : Sum1 += Geom.Inc.Lamda(J) * state.dataConstruction->Construct(IConst).BSDFInput.SolBkRefl(J, M);
2040 : } // Incident ray loop
2041 1920 : State.IntegratedBkRefl(J) = Sum1;
2042 : } // Outgoing ray loop
2043 :
2044 24 : if (!allocated(State.IntegratedBkTrans)) {
2045 24 : State.IntegratedBkTrans.allocate(Geom.Trn.NBasis);
2046 : }
2047 : // Integrating back transmittance
2048 1944 : for (J = 1; J <= Geom.Trn.NBasis; ++J) { // Outgoing ray loop
2049 1920 : State.IntegratedBkTrans(J) = 1 - State.IntegratedBkRefl(J) - State.IntegratedBkAbs(J);
2050 : } // Outgoing ray loop
2051 24 : }
2052 :
2053 1104 : Real64 SkyWeight([[maybe_unused]] Vector const &DirVec) // Direction of the element to be weighted
2054 : {
2055 :
2056 : // FUNCTION INFORMATION:
2057 : // AUTHOR Joe Klems
2058 : // DATE WRITTEN June 2011
2059 : // MODIFIED na
2060 : // RE-ENGINEERED na
2061 :
2062 : // PURPOSE OF THIS FUNCTION:
2063 : // Search a one-dimensional array for a given value, returning the index of the element equal to the value, if
2064 : // found, or zero
2065 :
2066 : using namespace Vectors;
2067 :
2068 : // Return value
2069 : Real64 Wt; // Weight
2070 :
2071 1104 : Wt = 1.0;
2072 :
2073 : // To do: figure out how to weight sky elements to reproduce the current E+ assumptions
2074 : // Possibly one will need to calculated average DH transmittance for isotropic sky and
2075 : // horizon separately and then later average according to sky conditions. Then a completely
2076 : // different scheme for daylight. For now: rays that reach sky equally weighted in calculating
2077 : // transmittance, rays passing through surfaces with scheduled transmittance are neglected.
2078 :
2079 1104 : return Wt;
2080 : }
2081 :
2082 816 : Real64 SkyGndWeight([[maybe_unused]] Vector const &PosVec) // x,y,z(=0) of ground intersection pt
2083 : {
2084 :
2085 : // FUNCTION INFORMATION:
2086 : // AUTHOR Joe Klems
2087 : // DATE WRITTEN June 2011
2088 : // MODIFIED na
2089 : // RE-ENGINEERED na
2090 :
2091 : // PURPOSE OF THIS FUNCTION:
2092 : // Search a one-dimensional array for a given value, returning the index of the element equal to the value, if
2093 : // found, or zero
2094 :
2095 : using namespace Vectors;
2096 :
2097 : // Return value
2098 : Real64 Wt; // Weight
2099 :
2100 816 : Wt = 1.0;
2101 :
2102 : // At present, equally weights all ground rays for calculation of the complex window transmittance for
2103 : // sky radiation reflected from ground. This does not take into account shading of the ground.
2104 : // The correct procedure would be to generate a set of rays to the sky and see which do not intersect
2105 : // surfaces, as is done in the reflection manager. However, this would increase computational load.
2106 : // Given that equal weighting, by averaging the transmittance only over rays that come from the ground,
2107 : // already produces a more accurate ground transmittance than the existing method, it is at least questionable
2108 : // whether the more detailed procedure would produce enough improvement in accuracy to make up for
2109 : // the additional calculation time. Therefore a more detailed treatment is deferred until there is some
2110 : // experience with the new method to determine whether further detail is warranted.
2111 :
2112 816 : return Wt;
2113 : }
2114 :
2115 4074 : BSDFDaylghtPosition DaylghtAltAndAzimuth(Vector const &UnitVect) // vector which needs to be converted
2116 : {
2117 :
2118 : // SUBROUTINE INFORMATION:
2119 : // AUTHOR Simon Vidanovic
2120 : // DATE WRITTEN April 2013
2121 : // MODIFIED na
2122 : // RE-ENGINEERED na
2123 :
2124 : // PURPOSE OF THIS FUNCTION:
2125 : // Transform unit vector (given in world coordinates) into altitude and azimuth. Azimuth is measured from positive x-axe.
2126 : // Altitude range is from -pi/2 to pi/2. Vector laying in horizontal plane will have altitude equal to zero and vector
2127 : // pointing upward will have altitude equal to pi/2. Range for azimuth is calculated from -pi to +pi.
2128 :
2129 : using namespace DataBSDFWindow;
2130 : using namespace Vectors;
2131 : // Return value
2132 4074 : BSDFDaylghtPosition DayPos; // altitude and azimuth in world coordinates
2133 :
2134 4074 : if (UnitVect.x != 0.0) {
2135 3620 : if (UnitVect.x >= 0.0) {
2136 1828 : DayPos.Azimuth = std::atan(UnitVect.y / UnitVect.x);
2137 : } else {
2138 1792 : if (UnitVect.y >= 0.0) {
2139 922 : DayPos.Azimuth = Constant::Pi + std::atan(UnitVect.y / UnitVect.x);
2140 : } else {
2141 870 : DayPos.Azimuth = -Constant::Pi + std::atan(UnitVect.y / UnitVect.x);
2142 : }
2143 : }
2144 : } else {
2145 454 : if (UnitVect.y >= 0.0) {
2146 182 : DayPos.Azimuth = Constant::PiOvr2;
2147 : } else {
2148 272 : DayPos.Azimuth = -Constant::PiOvr2;
2149 : }
2150 : }
2151 :
2152 4074 : DayPos.Altitude = std::asin(UnitVect.z);
2153 :
2154 4074 : return DayPos;
2155 : }
2156 :
2157 3840 : Vector WorldVectFromW6([[maybe_unused]] EnergyPlusData &state,
2158 : Real64 const Theta, // Polar angle in W6 Coords
2159 : Real64 const Phi, // Azimuthal angle in W6 Coords
2160 : const RayIdentificationType RadType, // Type of radiation: Front_Incident, etc.
2161 : Real64 const Gamma, // Surface tilt angle, radians, world coordinate system
2162 : Real64 const Alpha // Surface azimuth, radians, world coordinate system
2163 : )
2164 : {
2165 :
2166 : // SUBROUTINE INFORMATION:
2167 : // AUTHOR Joe Klems
2168 : // DATE WRITTEN Aug 2010
2169 : // MODIFIED na
2170 : // RE-ENGINEERED na
2171 :
2172 : // PURPOSE OF THIS FUNCTION:
2173 : // Transform angular coordinates in the WINDOW6 coordinate system for
2174 : // a given surface into a unit vector in the world coordinate system,
2175 : // pointing to the radiation source (for incident radiation) or in
2176 : // the direction of propagation (for outgoing radiation)
2177 :
2178 : using namespace Vectors;
2179 :
2180 : // Return value
2181 3840 : Vector UnitVect; // unit vector direction in world CS
2182 :
2183 : // Error tolerance is used to make small numbers equal to zero. Due to precision of pi constant used in E+, performing
2184 : // trigonometric operations on those constant will not cause absolutely accurate results
2185 3840 : Real64 constexpr ErrorTolerance(1.e-10);
2186 :
2187 3840 : UnitVect = Vector(0.0, 0.0, 0.0);
2188 :
2189 3840 : Real64 const sin_Phi = std::sin(Phi);
2190 3840 : Real64 const cos_Phi = std::cos(Phi);
2191 :
2192 3840 : Real64 const sin_Gamma = std::sin(Gamma);
2193 3840 : Real64 const cos_Gamma = std::cos(Gamma);
2194 :
2195 3840 : Real64 const sin_Alpha = std::sin(Alpha);
2196 3840 : Real64 const cos_Alpha = std::cos(Alpha);
2197 :
2198 3840 : Real64 const sin_Theta = std::sin(Theta);
2199 3840 : Real64 const cos_Theta = std::cos(Theta);
2200 :
2201 3840 : switch (RadType) {
2202 1920 : case RayIdentificationType::Front_Incident: { // W6 vector will point in direction of propagation, must reverse to get world vector
2203 : // after the W6 vector has been rotated into the world CS
2204 1920 : UnitVect.x = sin_Theta * sin_Phi * cos_Gamma * sin_Alpha - sin_Theta * cos_Phi * cos_Alpha + cos_Theta * sin_Gamma * sin_Alpha;
2205 1920 : UnitVect.y = sin_Theta * cos_Phi * sin_Alpha + sin_Theta * sin_Phi * cos_Gamma * cos_Alpha + cos_Theta * sin_Gamma * cos_Alpha;
2206 1920 : UnitVect.z = -(sin_Theta * sin_Phi * sin_Gamma - cos_Theta * cos_Gamma);
2207 1920 : break;
2208 : }
2209 1920 : case RayIdentificationType::Front_Transmitted: {
2210 1920 : UnitVect.x = sin_Theta * cos_Phi * cos_Alpha - sin_Theta * sin_Phi * cos_Gamma * sin_Alpha - cos_Theta * sin_Gamma * sin_Alpha;
2211 1920 : UnitVect.y = -(sin_Theta * cos_Phi * sin_Alpha + sin_Theta * sin_Phi * cos_Gamma * cos_Alpha + cos_Theta * sin_Gamma * cos_Alpha);
2212 1920 : UnitVect.z = sin_Theta * sin_Phi * sin_Gamma - cos_Theta * cos_Gamma;
2213 1920 : break;
2214 : }
2215 0 : case RayIdentificationType::Front_Reflected: {
2216 0 : UnitVect.x = sin_Theta * cos_Phi * cos_Alpha - sin_Theta * sin_Phi * cos_Gamma * sin_Alpha + cos_Theta * sin_Gamma * sin_Alpha;
2217 0 : UnitVect.y = cos_Theta * sin_Gamma * cos_Alpha - sin_Theta * cos_Phi * sin_Alpha - sin_Theta * sin_Phi * cos_Gamma * cos_Alpha;
2218 0 : UnitVect.z = sin_Theta * sin_Phi * sin_Gamma + cos_Theta * cos_Gamma;
2219 0 : break;
2220 : }
2221 0 : case RayIdentificationType::Back_Incident: {
2222 0 : UnitVect.x = sin_Theta * sin_Phi * cos_Gamma * sin_Alpha - sin_Theta * cos_Phi * cos_Alpha - cos_Theta * sin_Gamma * sin_Alpha;
2223 0 : UnitVect.y = sin_Theta * cos_Phi * sin_Alpha + sin_Theta * sin_Phi * cos_Gamma * cos_Alpha - cos_Theta * sin_Gamma * cos_Alpha;
2224 0 : UnitVect.z = -cos_Theta * cos_Gamma - sin_Theta * sin_Phi * sin_Gamma;
2225 0 : break;
2226 : }
2227 0 : case RayIdentificationType::Back_Transmitted: { // This is same as front reflected
2228 0 : UnitVect.x = sin_Theta * cos_Phi * cos_Alpha - sin_Theta * sin_Phi * cos_Gamma * sin_Alpha + cos_Theta * sin_Gamma * sin_Alpha;
2229 0 : UnitVect.y = cos_Theta * sin_Gamma * cos_Alpha - sin_Theta * cos_Phi * sin_Alpha - sin_Theta * sin_Phi * cos_Gamma * cos_Alpha;
2230 0 : UnitVect.z = sin_Theta * sin_Phi * sin_Gamma + cos_Theta * cos_Gamma;
2231 0 : break;
2232 : }
2233 0 : case RayIdentificationType::Back_Reflected: { // This is same as front transmitted
2234 0 : UnitVect.x = sin_Theta * cos_Phi * cos_Alpha - sin_Theta * sin_Phi * cos_Gamma * cos_Alpha - cos_Theta * sin_Gamma * sin_Alpha;
2235 0 : UnitVect.y = -(sin_Theta * cos_Phi * sin_Alpha + sin_Theta * sin_Phi * cos_Gamma * cos_Alpha + cos_Theta * sin_Gamma * cos_Alpha);
2236 0 : UnitVect.z = sin_Theta * sin_Phi * sin_Gamma - cos_Theta * cos_Gamma;
2237 0 : break;
2238 : }
2239 0 : default:
2240 0 : break;
2241 : }
2242 :
2243 : // Remove small numbers from evaluation (due to limited decimal points for pi)
2244 3840 : if (std::abs(UnitVect.x) <= ErrorTolerance) {
2245 364 : UnitVect.x = 0.0;
2246 : }
2247 3840 : if (std::abs(UnitVect.y) <= ErrorTolerance) {
2248 212 : UnitVect.y = 0.0;
2249 : }
2250 3840 : if (std::abs(UnitVect.z) <= ErrorTolerance) {
2251 576 : UnitVect.z = 0.0;
2252 : }
2253 :
2254 3840 : return UnitVect;
2255 : }
2256 :
2257 2833 : int FindInBasis(EnergyPlusData &state,
2258 : Vector const &RayToFind, // Ray vector direction in world CS
2259 : const RayIdentificationType RadType, // Type of radiation: Front_Incident, etc.
2260 : int const ISurf, // Window Surface number
2261 : [[maybe_unused]] int const IState, // Complex Fenestration state number
2262 : BasisStruct const &Basis, // Complex Fenestration basis root
2263 : Real64 &Theta, // Theta value for ray
2264 : Real64 &Phi // Phi value for ray
2265 : )
2266 : {
2267 :
2268 : // SUBROUTINE INFORMATION:
2269 : // AUTHOR Joe Klems
2270 : // DATE WRITTEN August 2011
2271 : // MODIFIED na
2272 : // RE-ENGINEERED na
2273 :
2274 : using namespace Vectors;
2275 :
2276 : // Return value
2277 : int RayIndex; // Index of ray in basis, zero if ray not in hemisphere
2278 :
2279 : // INTEGER, INTENT(IN) :: IWind !window index in window list
2280 :
2281 : int ITheta; // Table index of Theta
2282 : int IPhi; // Table index of Phi, given ITheta
2283 : int IThDn; // Theta lower table index
2284 : int IThUp; // Theta upper table index
2285 : int IPhDn; // Phi lower table index
2286 : int IPhUp; // Phi upper table index
2287 : Real64 Gamma; // Gamma (tilt) angle of window
2288 : Real64 Alpha; // Alpha (azimuth) angle of window
2289 : Real64 DotProd;
2290 :
2291 2833 : Theta = 0.0;
2292 2833 : Phi = 0.0;
2293 :
2294 : // Check if surface and vector are pointing in different directions
2295 2833 : DotProd = dot(RayToFind, state.dataSurface->Surface(ISurf).NewellSurfaceNormalVector);
2296 2833 : if (DotProd <= 0.0) {
2297 1158 : RayIndex = 0;
2298 1158 : return RayIndex;
2299 : }
2300 :
2301 : // get window tilt and azimuth
2302 1675 : Gamma = Constant::DegToRad * state.dataSurface->Surface(ISurf).Tilt;
2303 1675 : Alpha = Constant::DegToRad * state.dataSurface->Surface(ISurf).Azimuth;
2304 : // get the corresponding local Theta, Phi for ray
2305 1675 : W6CoordsFromWorldVect(state, RayToFind, RadType, Gamma, Alpha, Theta, Phi);
2306 :
2307 1675 : if (Theta >= 0.5 * Constant::Pi) { // Ray was in not in correct hemisphere
2308 0 : RayIndex = 0;
2309 0 : return RayIndex;
2310 : }
2311 1675 : if (Basis.BasisSymmetryType == DataBSDFWindow::BasisSymmetry::None) {
2312 : // Search the basis thetas
2313 1675 : if (Theta <= 0.0) {
2314 : // Special case, Theta = 0.; this is always the first basis element
2315 0 : RayIndex = 1;
2316 0 : return RayIndex;
2317 : }
2318 : // So here Theta > 0
2319 : // Note the table searches always go to the limit point, which is not itself a basis element
2320 1675 : IThUp = SearchAscTable(Theta, Basis.NThetas + 1, Basis.Thetas);
2321 1675 : IThDn = IThUp - 1;
2322 : // Determine which of the theta basis points is closer to the Theta value
2323 1675 : if (Theta <= Basis.Grid(Basis.BasisIndex(1, IThDn)).UpprTheta) {
2324 : // Note this will take care of both the special cases IThUp=2 and IThUp=NThetas +1
2325 872 : ITheta = IThDn;
2326 : } else {
2327 803 : ITheta = IThUp;
2328 : }
2329 : // Now determine the Phi index
2330 1675 : if (Basis.NPhis(ITheta) == 1) {
2331 : // Note that for W6 basis this can only happen for the first basis element
2332 : // If later bases are introduced this logic may have to be redesigned
2333 0 : RayIndex = Basis.BasisIndex(1, ITheta);
2334 0 : return RayIndex;
2335 : }
2336 1675 : IPhUp = SearchAscTable(Phi, Basis.NPhis(ITheta) + 1, Basis.Phis(ITheta, _));
2337 1675 : IPhDn = IPhUp - 1;
2338 1675 : if (Phi <= Basis.Grid(Basis.BasisIndex(IPhDn, ITheta)).UpprPhi) {
2339 862 : IPhi = IPhDn;
2340 : } else {
2341 813 : if (IPhUp == Basis.NPhis(ITheta) + 1) {
2342 : // Phi is above upper limit for highest Phi basis element, meaning it is closer to 2Pi,
2343 : // i.e., the first element
2344 114 : IPhi = 1;
2345 : } else {
2346 699 : IPhi = IPhUp;
2347 : }
2348 : }
2349 1675 : RayIndex = Basis.BasisIndex(IPhi, ITheta);
2350 1675 : return RayIndex;
2351 0 : } else if (Basis.BasisSymmetryType == DataBSDFWindow::BasisSymmetry::Axisymmetric) {
2352 : // Search the basis thetas
2353 0 : if (Theta <= 0.0) {
2354 : // Special case, Theta = 0.; this is always the first basis element
2355 0 : RayIndex = 1;
2356 0 : return RayIndex;
2357 : }
2358 : // So here Theta > 0
2359 : // Note the table searches always go to the limit point, which is not itself a basis element
2360 0 : IThUp = SearchAscTable(Theta, Basis.NThetas + 1, Basis.Thetas);
2361 0 : IThDn = IThUp - 1;
2362 : // Determine which of the theta basis points is closer to the Theta value
2363 0 : if (Theta <= Basis.Grid(Basis.BasisIndex(1, IThDn)).UpprTheta) {
2364 : // Note this will take care of both the special cases IThUp=2 and IThUp=NThetas +1
2365 0 : ITheta = IThDn;
2366 : } else {
2367 0 : ITheta = IThUp;
2368 : }
2369 0 : RayIndex = Basis.BasisIndex(1, ITheta);
2370 0 : return RayIndex;
2371 : }
2372 : // No other type is implemented
2373 0 : RayIndex = 0;
2374 :
2375 0 : return RayIndex;
2376 : }
2377 :
2378 3327 : void W6CoordsFromWorldVect([[maybe_unused]] EnergyPlusData &state,
2379 : Vector const &RayVect, // Ray vector direction in world CS
2380 : const RayIdentificationType RadType, // Type of radiation: Front_Incident, etc.
2381 : Real64 const Gamma, // Surface tilt angle, world coordinate system
2382 : Real64 const Alpha, // Surface azimuth, world coordinate system
2383 : Real64 &Theta, // Polar angle in W6 Coords
2384 : Real64 &Phi // Azimuthal angle in W6 Coords
2385 : )
2386 : {
2387 :
2388 : // SUBROUTINE INFORMATION:
2389 : // AUTHOR Joe Klems
2390 : // DATE WRITTEN August 2011
2391 : // MODIFIED na
2392 : // RE-ENGINEERED na
2393 :
2394 : // PURPOSE OF THIS FUNCTION:
2395 : // Invert the transformation from W6 to world coordinates to
2396 : // calculate the theta, phi corresponding to a given ray direction
2397 : // in the world coordinate system, for a window with a
2398 : // given rotation and tilt (Gamma and Alpha)
2399 : // (needed for locating the sun direction in the local coordinate system)
2400 :
2401 : using namespace Vectors;
2402 :
2403 3327 : Real64 Cost(0.0); // Temp for cos theta
2404 : Real64 Sint; // Temp for sin theta
2405 : Real64 Psi; // Temp for phi before rotation adjustment
2406 : Real64 RdotX; // Temp variable for manipulating .dot. produt
2407 : Real64 RdotY; // Temp variable for manipulating .dot. produt
2408 : Real64 RdotZ; // Temp variable for manipulating .dot. produt
2409 :
2410 : // Object Data
2411 3327 : Vector W6x; // W6 x coordinate unit vector
2412 3327 : Vector W6y; // W6 y coordinate unit vector
2413 3327 : Vector W6z; // W6 z coordinate unit vector
2414 :
2415 : // define the local W6 coordinate vectors
2416 3327 : W6x.x = std::cos(Alpha);
2417 3327 : W6x.y = -std::sin(Alpha);
2418 3327 : W6x.z = 0.0;
2419 3327 : W6y.x = -std::cos(Gamma) * std::sin(Alpha);
2420 3327 : W6y.y = -std::cos(Gamma) * std::cos(Alpha);
2421 3327 : W6y.z = std::sin(Gamma);
2422 3327 : W6z.x = -std::sin(Gamma) * std::sin(Alpha);
2423 3327 : W6z.y = -std::sin(Gamma) * std::cos(Alpha);
2424 3327 : W6z.z = -std::cos(Gamma);
2425 :
2426 3327 : switch (RadType) {
2427 3327 : case RayIdentificationType::Front_Incident: {
2428 3327 : RdotZ = dot(W6z, RayVect);
2429 3327 : Cost = -RdotZ;
2430 3327 : Sint = std::sqrt(1.0 - pow_2(Cost));
2431 3327 : Theta = std::acos(Cost);
2432 3327 : RdotY = dot(W6y, RayVect);
2433 3327 : RdotX = dot(W6x, RayVect);
2434 3327 : Psi = std::atan2(-RdotY / Sint, -RdotX / Sint);
2435 3327 : if (Psi < 0.0) {
2436 3327 : Phi = 2.0 * Constant::Pi + Psi;
2437 : } else {
2438 0 : Phi = Psi;
2439 : }
2440 3327 : break;
2441 : }
2442 0 : case RayIdentificationType::Front_Transmitted: {
2443 0 : Cost = dot(W6z, RayVect);
2444 0 : Sint = std::sqrt(1.0 - pow_2(Cost));
2445 0 : Theta = std::acos(Cost);
2446 0 : RdotY = dot(W6y, RayVect);
2447 0 : RdotX = dot(W6x, RayVect);
2448 0 : Psi = std::atan2(RdotY / Sint, RdotX / Sint);
2449 0 : if (Psi < 0.0) {
2450 0 : Phi = 2.0 * Constant::Pi + Psi;
2451 : } else {
2452 0 : Phi = Psi;
2453 : }
2454 0 : break;
2455 : }
2456 0 : case RayIdentificationType::Front_Reflected: {
2457 0 : RdotZ = dot(W6z, RayVect);
2458 0 : Cost = -RdotZ;
2459 0 : Sint = std::sqrt(1.0 - pow_2(Cost));
2460 0 : Theta = std::acos(Cost);
2461 0 : RdotY = dot(W6y, RayVect);
2462 0 : RdotX = dot(W6x, RayVect);
2463 0 : Psi = std::atan2(RdotY / Sint, RdotX / Sint);
2464 0 : if (Psi < 0.0) {
2465 0 : Phi = 2.0 * Constant::Pi + Psi;
2466 : } else {
2467 0 : Phi = Psi;
2468 : }
2469 0 : break;
2470 : }
2471 0 : case RayIdentificationType::Back_Incident: {
2472 0 : Cost = dot(W6z, RayVect);
2473 0 : Sint = std::sqrt(1.0 - pow_2(Cost));
2474 0 : Theta = std::acos(Cost);
2475 0 : RdotY = dot(W6y, RayVect);
2476 0 : RdotX = dot(W6x, RayVect);
2477 0 : Psi = std::atan2(-RdotY / Sint, -RdotX / Sint);
2478 0 : if (Psi < 0.0) {
2479 0 : Phi = 2 * Constant::Pi + Psi;
2480 : } else {
2481 0 : Phi = Psi;
2482 : }
2483 0 : break;
2484 : }
2485 0 : case RayIdentificationType::Back_Transmitted: { // This is same as front reflected
2486 0 : RdotZ = dot(W6z, RayVect);
2487 0 : Cost = -RdotZ;
2488 0 : Sint = std::sqrt(1.0 - pow_2(Cost));
2489 0 : Theta = std::acos(Cost);
2490 0 : RdotY = dot(W6y, RayVect);
2491 0 : RdotX = dot(W6x, RayVect);
2492 0 : Psi = std::atan2(RdotY / Sint, RdotX / Sint);
2493 0 : if (Psi < 0.0) {
2494 0 : Phi = 2.0 * Constant::Pi + Psi;
2495 : } else {
2496 0 : Phi = Psi;
2497 : }
2498 0 : break;
2499 : }
2500 0 : case RayIdentificationType::Back_Reflected: { // This is same as front transmitted
2501 0 : Cost = dot(W6z, RayVect);
2502 0 : Sint = std::sqrt(1.0 - pow_2(Cost));
2503 0 : Theta = std::acos(Cost);
2504 0 : RdotY = dot(W6y, RayVect);
2505 0 : RdotX = dot(W6x, RayVect);
2506 0 : Psi = std::atan2(RdotY / Sint, RdotX / Sint);
2507 0 : if (Psi < 0.0) {
2508 0 : Phi = 2.0 * Constant::Pi + Psi;
2509 : } else {
2510 0 : Phi = Psi;
2511 : }
2512 0 : break;
2513 : }
2514 0 : default:
2515 0 : assert(false);
2516 : break;
2517 : }
2518 3327 : if (std::abs(Cost) < Constant::rTinyValue) {
2519 0 : Cost = 0.0;
2520 : }
2521 3327 : if (Cost < 0.0) {
2522 672 : Theta = Constant::Pi - Theta; // This signals ray out of hemisphere
2523 : }
2524 3327 : }
2525 :
2526 37990 : void CalcComplexWindowThermal(EnergyPlusData &state,
2527 : int const SurfNum, // Surface number
2528 : int &ConstrNum, // Construction number
2529 : Real64 const HextConvCoeff, // Outside air film conductance coefficient
2530 : Real64 &SurfInsideTemp, // Inside window surface temperature
2531 : Real64 &SurfOutsideTemp, // Outside surface temperature (C)
2532 : Real64 &SurfOutsideEmiss,
2533 : DataBSDFWindow::Condition const CalcCondition // Calucation condition (summer, winter or no condition)
2534 : )
2535 : {
2536 :
2537 : // SUBROUTINE INFORMATION:
2538 : // AUTHOR B. Griffith
2539 : // DATE WRITTEN October 2009
2540 : // MODIFIED Simon Vidanovic
2541 : // RE-ENGINEERED September 2011
2542 :
2543 : // PURPOSE OF THIS SUBROUTINE:
2544 : // wrapper between E+ and TARCOG
2545 :
2546 : // METHODOLOGY EMPLOYED:
2547 : // draft out an attempt for proof-of-concept, to reuse native TARCOG implementation
2548 : // based off of 1-26-2009 version of WinCOG/TARCOG solution from Carli, Inc.
2549 :
2550 : using namespace DataBSDFWindow;
2551 : using Psychrometrics::PsyCpAirFnW;
2552 : using Psychrometrics::PsyTdpFnWPb;
2553 : using TARCOGGassesParams::maxgas;
2554 : using TARCOGMain::TARCOG90;
2555 : using TARCOGParams::maxlay;
2556 : using TARCOGParams::maxlay1;
2557 :
2558 : // Locals
2559 : // SUBROUTINE ARGUMENT DEFINITIONS:
2560 : // (temperature of innermost face) [C]
2561 : // INTEGER, INTENT(IN) :: CurrentThermalModelNumber
2562 :
2563 : // SUBROUTINE PARAMETER DEFINITIONS:
2564 : // INTEGER, PARAMETER :: maxlay = 100 ! maximum number of layers (including laminates)
2565 : // INTEGER, PARAMETER :: maxgas = 10 ! maximum number of individual gasses
2566 : // INTEGER, PARAMETER :: maxlay1 = maxlay+1 ! maximum number of 'gaps', including in and out (maxlay+1)
2567 : // REAL(r64), PARAMETER :: StefanBoltzmannConst = 5.6697d-8 ! Stefan-Boltzmann constant in W/(m2*K4)
2568 :
2569 : // TARCOG Inputs:
2570 37990 : int nlayer(0); // Number of glazing layers
2571 37990 : int iwd(0); // Wind direction: 0 - windward, 1 - leeward
2572 37990 : Real64 tout(0.0); // Outdoor temperature [K]
2573 37990 : Real64 tind(0.0); // Indoor temperature [K]
2574 37990 : Real64 trmin(0.0); // Indoor mean radiant temperature [K]
2575 37990 : Real64 wso(0.0); // Outdoor wind speed [m/s]
2576 37990 : Real64 wsi(0.0); // Inside forced air speed [m/s]
2577 37990 : Real64 dir(0.0); // Direct solar radiation [W/m^2]
2578 37990 : int isky(0); // Flag for sky temperature (Tsky) and sky emittance (esky)
2579 : // 0 - both tsky and esky are specified
2580 : // 1 - tsky specified, esky = 1
2581 : // 2 - Swinbank model for effective sky emittance
2582 37990 : Real64 tsky(0.0); // Night sky temperature [K]
2583 37990 : Real64 esky(0.0); // Effective night sky emittance
2584 37990 : Real64 fclr(0.0); // Fraction of sky that is clear
2585 : Real64 VacuumPressure; // maximal pressure for gas to be considered as vacuum [Pa]
2586 : Real64 VacuumMaxGapThickness; // maximal gap thickness for which vacuum calculation will work without issuing
2587 : // warning message
2588 :
2589 37990 : auto &gap = state.dataWindowComplexManager->gap;
2590 37990 : auto &thick = state.dataWindowComplexManager->thick;
2591 37990 : auto &scon = state.dataWindowComplexManager->scon;
2592 37990 : auto &tir = state.dataWindowComplexManager->tir;
2593 37990 : auto &emis = state.dataWindowComplexManager->emis;
2594 37990 : auto &SupportPlr = state.dataWindowComplexManager->SupportPlr;
2595 37990 : auto &PillarSpacing = state.dataWindowComplexManager->PillarSpacing;
2596 37990 : auto &PillarRadius = state.dataWindowComplexManager->PillarRadius;
2597 37990 : auto &asol = state.dataWindowComplexManager->asol;
2598 37990 : auto &presure = state.dataWindowComplexManager->presure;
2599 37990 : auto &GapDefMax = state.dataWindowComplexManager->GapDefMax;
2600 37990 : auto &YoungsMod = state.dataWindowComplexManager->YoungsMod;
2601 37990 : auto &PoissonsRat = state.dataWindowComplexManager->PoissonsRat;
2602 37990 : auto &LayerDef = state.dataWindowComplexManager->LayerDef;
2603 37990 : auto &iprop = state.dataWindowComplexManager->iprop;
2604 37990 : auto &frct = state.dataWindowComplexManager->frct;
2605 37990 : auto &gcon = state.dataWindowComplexManager->gcon;
2606 37990 : auto &gvis = state.dataWindowComplexManager->gvis;
2607 37990 : auto &gcp = state.dataWindowComplexManager->gcp;
2608 37990 : auto &wght = state.dataWindowComplexManager->wght;
2609 37990 : auto &gama = state.dataWindowComplexManager->gama;
2610 37990 : auto &nmix = state.dataWindowComplexManager->nmix;
2611 37990 : auto &ibc = state.dataWindowComplexManager->ibc;
2612 37990 : auto &Atop = state.dataWindowComplexManager->Atop;
2613 37990 : auto &Abot = state.dataWindowComplexManager->Abot;
2614 37990 : auto &Al = state.dataWindowComplexManager->Al;
2615 37990 : auto &Ar = state.dataWindowComplexManager->Ar;
2616 37990 : auto &Ah = state.dataWindowComplexManager->Ah;
2617 37990 : auto &SlatThick = state.dataWindowComplexManager->SlatThick;
2618 37990 : auto &SlatWidth = state.dataWindowComplexManager->SlatWidth;
2619 37990 : auto &SlatAngle = state.dataWindowComplexManager->SlatAngle;
2620 37990 : auto &SlatCond = state.dataWindowComplexManager->SlatCond;
2621 37990 : auto &SlatSpacing = state.dataWindowComplexManager->SlatSpacing;
2622 37990 : auto &SlatCurve = state.dataWindowComplexManager->SlatCurve;
2623 37990 : auto &vvent = state.dataWindowComplexManager->vvent;
2624 37990 : auto &tvent = state.dataWindowComplexManager->tvent;
2625 37990 : auto &LayerType = state.dataWindowComplexManager->LayerType;
2626 37990 : auto &nslice = state.dataWindowComplexManager->nslice;
2627 37990 : auto &LaminateA = state.dataWindowComplexManager->LaminateA;
2628 37990 : auto &LaminateB = state.dataWindowComplexManager->LaminateB;
2629 37990 : auto &sumsol = state.dataWindowComplexManager->sumsol;
2630 37990 : auto &theta = state.dataWindowComplexManager->theta;
2631 37990 : auto &q = state.dataWindowComplexManager->q;
2632 37990 : auto &qv = state.dataWindowComplexManager->qv;
2633 37990 : auto &hcgap = state.dataWindowComplexManager->hcgap;
2634 37990 : auto &hrgap = state.dataWindowComplexManager->hrgap;
2635 37990 : auto &hg = state.dataWindowComplexManager->hg;
2636 37990 : auto &hr = state.dataWindowComplexManager->hr;
2637 37990 : auto &hs = state.dataWindowComplexManager->hs;
2638 37990 : auto &Ra = state.dataWindowComplexManager->Ra;
2639 37990 : auto &Nu = state.dataWindowComplexManager->Nu;
2640 37990 : auto &Keff = state.dataWindowComplexManager->Keff;
2641 37990 : auto &ShadeGapKeffConv = state.dataWindowComplexManager->ShadeGapKeffConv;
2642 :
2643 37990 : Real64 totsol(0.0); // Total solar transmittance of the IGU
2644 37990 : Real64 tilt(0.0); // Window tilt [degrees]
2645 37990 : Real64 height(0.0); // IGU cavity height [m]
2646 37990 : Real64 heightt(0.0); // Total window height [m]
2647 37990 : Real64 width(0.0); // Window width [m]
2648 :
2649 : // Deflection
2650 : // Tarcog requires deflection as input parameters. Deflection is NOT used in EnergyPlus simulations
2651 : TARCOGParams::DeflectionCalculation CalcDeflection; // Deflection calculation flag:
2652 : // 0 - no deflection calculations
2653 : // 1 - perform deflection calculation (input is Pressure/Temp)
2654 : // 2 - perform deflection calculation (input is measured deflection)
2655 : Real64 Pa; // Atmospheric (outside/inside) pressure (used onlu if CalcDeflection = 1)
2656 : Real64 Pini; // Initial presssure at time of fabrication (used only if CalcDeflection = 1)
2657 : Real64 Tini; // Initial temperature at time of fabrication (used only if CalcDeflection = 1)
2658 37990 : Real64 hin(0.0); // Indoor combined film coefficient (if non-zero) [W/m^2.K]
2659 37990 : Real64 hout(0.0); // Outdoor combined film coefficient (if non-zero) [W/m^2.K]
2660 37990 : TARCOGGassesParams::Stdrd standard(TARCOGGassesParams::Stdrd::ISO15099); // Calculation standard switch:
2661 : // 1 - ISO 15099,
2662 : // 2 - EN673 / ISO 10292 Declared,
2663 : // 3 - EN673 / ISO 10292 Design.
2664 37990 : TARCOGParams::TARCOGThermalModel ThermalMod = TARCOGParams::TARCOGThermalModel::ISO15099; // Thermal model:
2665 : // 0 - ISO15099
2666 : // 1 - Scaled Cavity Width (SCW)
2667 : // 2 - Convective Scalar Model (CSM)
2668 37990 : std::string Debug_dir; // Target directory for debug files (pointer to a character array)
2669 37990 : std::string Debug_file("Test"); // Template file name used to create debug output files
2670 37990 : std::int32_t Window_ID(-1); // ID of the window (long integer value, passed by W6)
2671 37990 : std::int32_t IGU_ID(-1); // ID of the IGU (long integer value, passed by W6)
2672 37990 : Real64 SDScalar(0.0); // SD convection factor (value between 0 and 1)
2673 : // 0.0 - No SD layer
2674 : // 1.0 - Closed SD
2675 : // Notes: * vvent, tvent, Atop, Abot, Al, Ar and Ah are considered for SD layers only.
2676 : // ** SlatThick, SlatWidth, SlatAngle, SlatCond, SlatSpacing, SlatCurve
2677 : // are used for Venetian blind layers only.
2678 : // *** For vvent & tvent: vvent(1) - exterior, vvent(nlayer+1) - interior.
2679 : // **** Forced ventilation calculation is not active at this time.
2680 : // TARCOG Output:
2681 :
2682 37990 : Real64 ufactor(0.0); // Center of glass U-value [W/m^2.K]
2683 37990 : Real64 sc(0.0); // Shading Coefficient
2684 37990 : Real64 hflux(0.0); // Net heat flux between room and window [W/m^2]
2685 37990 : Real64 hcin(0.0); // Indoor convective surface heat transfer coefficient [W/m^2.K]
2686 37990 : Real64 hcout(0.0); // Outdoor convective surface heat transfer coefficient [W/m^2.K]
2687 37990 : Real64 hrin(0.0); // Indoor radiative surface heat transfer coefficient [W/m^2.K]
2688 37990 : Real64 hrout(0.0); // Outdoor radiative surface heat transfer coefficient [W/m^2.K]
2689 37990 : Real64 shgc(0.0); // Solar heat gain coefficient - per ISO 15099
2690 37990 : Real64 shgct(0.0); // Solar heat gain coefficient - per old procedure
2691 37990 : Real64 tamb(0.0); // Outdoor environmental temperature [K]
2692 37990 : Real64 troom(0.0); // Indoor environmental temperature [K]
2693 37990 : Real64 he(0.0); // External heat transfer coefficient [W/m^2.K] - EN673 and ISO 10292 procedure
2694 37990 : Real64 hi(0.0); // Internal heat transfer coefficient [W/m^2.K] - EN673 and ISO 10292 procedure
2695 37990 : int nperr(0); // Error code
2696 37990 : Real64 ShadeEmisRatioOut(0.0); // Ratio of modified to glass emissivity at the outermost glazing surface
2697 37990 : Real64 ShadeEmisRatioIn(0.0); // Ratio of modified to glass emissivity at the innermost glazing surface
2698 37990 : Real64 ShadeHcRatioOut(0.0); // Ratio of modified to unshaded Hc at the outermost glazing surface
2699 37990 : Real64 ShadeHcRatioIn(0.0); // Ratio of modified to unshaded Hc at the innermost glazing surface
2700 37990 : Real64 HcUnshadedOut(0.0); // Hc value at outdoor surface of an unshaded subsystem [W/m^2.K]
2701 37990 : Real64 HcUnshadedIn(0.0); // Hc value at indoor surface of an unshaded subsystem [W/m^2.K]
2702 :
2703 : int ZoneNum; // Zone number corresponding to SurfNum
2704 :
2705 : int TotLay; // Total number of layers in a construction
2706 : // (sum of solid layers and gap layers)
2707 : int Lay; // Layer number
2708 : int LayPtr; // Material number for a layer
2709 : int IGlass; // glass layer number (1,2,3,...)
2710 : int IGap; // Gap layer number (1,2,...)
2711 : int TotGlassLay; // Total number of glass layers in a construction
2712 : int k; // Layer counter
2713 : int SurfNumAdj; // An interzone surface's number in the adjacent zone
2714 : WinShadingType ShadeFlag; // Flag indicating whether shade or blind is on, and shade/blind position
2715 : int IMix;
2716 :
2717 : Real64 IncidentSolar; // Solar incident on outside of window (W)
2718 : Real64 ConvHeatFlowNatural; // Convective heat flow from gap between glass and interior shade or blind (W)
2719 : Real64 ShadeArea; // shade/blind area (m2)
2720 : Real64 sconsh; // shade/blind conductance (W/m2-K)
2721 : Real64 CondHeatGainShade; // Conduction through shade/blind, outside to inside (W)
2722 :
2723 : Real64 ShGlReflFacIR; // Factor for long-wave inter-reflection between shade/blind and adjacent glass
2724 : // Real64 RhoGlIR1; // Long-wave reflectance of glass surface facing shade/blind; 1=exterior shade/blind,
2725 : Real64 RhoGlIR2;
2726 : // 2=interior shade/blind
2727 : Real64 RhoShIR1; // Long-wave reflectance of shade/blind surface facing glass; 1=interior shade/blind,
2728 : Real64 RhoShIR2;
2729 : // 2=exterior shade/blind
2730 : Real64 EpsShIR1; // Long-wave emissivity of shade/blind surface facing glass; 1=interior shade/blind,
2731 : Real64 EpsShIR2;
2732 : // 2=exterior shade/blind
2733 : Real64 TauShIR; // Long-wave transmittance of isolated shade/blind
2734 : Real64 NetIRHeatGainShade; // Net IR heat gain to zone from interior shade/blind (W)
2735 : Real64 NetIRHeatGainGlass; // Net IR heat gain to zone from shade/blind side of glass when interior
2736 : // shade/blind is present. Zero if shade/blind has zero IR transmittance (W)
2737 : Real64 ConvHeatGainFrZoneSideOfShade; // Convective heat gain to zone from side of interior shade facing zone (W)
2738 : Real64 ConvHeatGainFrZoneSideOfGlass; // Convective heat gain to zone from side of glass facing zone when
2739 : // no interior shade/blind is present (W)
2740 : Real64 CondHeatGainGlass; // Conduction through inner glass layer, outside to inside (W)
2741 : Real64 TotAirflowGap; // Total volumetric airflow through window gap (m3/s)
2742 : Real64 TAirflowGapOutlet; // Temperature of air leaving airflow gap between glass panes (K)
2743 : Real64 TAirflowGapOutletC; // Temperature of air leaving airflow gap between glass panes (C)
2744 : Real64 ConvHeatFlowForced; // Convective heat flow from forced airflow gap (W)
2745 : Real64 InletAirHumRat; // Humidity ratio of air from window gap entering fan
2746 : Real64 ZoneTemp; // Zone air temperature (C)
2747 : Real64 CpAirOutlet; // Heat capacity of air from window gap (J/kg-K)
2748 : Real64 CpAirZone; // Heat capacity of zone air (J/kg-K)
2749 : Real64 ConvHeatGainToZoneAir; // Convective heat gain to zone air from window gap airflow (W)
2750 : // int ConstrNumSh; // Construction number with shading device
2751 37990 : int CalcSHGC(0); // SHGC calculations are not necessary for E+ run
2752 37990 : int NumOfIterations(0);
2753 :
2754 37990 : std::string tarcogErrorMessage; // store error text from tarcog
2755 :
2756 : // Simon: locally used variables
2757 : int ngllayer;
2758 : int nglface;
2759 : int nglfacep;
2760 : int ThermalModelNum;
2761 : Real64 rmir; // IR radiance of window's interior surround (W/m2)
2762 : Real64 outir;
2763 : Real64 Ebout;
2764 : Real64 dominantGapWidth; // store value for dominant gap width. Used for airflow calculations
2765 : Real64 edgeGlCorrFac;
2766 :
2767 : Real64 SrdSurfTempAbs; // Absolute temperature of a surrounding surface
2768 : Real64 OutSrdIR;
2769 :
2770 37990 : auto &s_mat = state.dataMaterial;
2771 : // fill local vars
2772 :
2773 37990 : CalcDeflection = TARCOGParams::DeflectionCalculation::NONE;
2774 37990 : CalcSHGC = 0;
2775 :
2776 37990 : if (CalcCondition == DataBSDFWindow::Condition::Invalid) {
2777 37968 : ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
2778 37968 : SurfNumAdj = state.dataSurface->Surface(SurfNum).ExtBoundCond;
2779 37968 : ShadeFlag = state.dataSurface->SurfWinShadingFlag(SurfNum);
2780 : }
2781 :
2782 37990 : TotGlassLay = state.dataConstruction->Construct(ConstrNum).TotGlassLayers;
2783 37990 : ngllayer = state.dataConstruction->Construct(ConstrNum).TotGlassLayers;
2784 37990 : nglface = 2 * ngllayer;
2785 37990 : nglfacep = nglface;
2786 37990 : hrin = 0.0;
2787 37990 : hcin = 0.0;
2788 37990 : hrout = 0.0;
2789 37990 : hcout = 0.0;
2790 :
2791 37990 : Pa = state.dataEnvrn->OutBaroPress;
2792 :
2793 37990 : ThermalModelNum = state.dataConstruction->Construct(ConstrNum).BSDFInput.ThermalModel;
2794 37990 : standard = state.dataMaterial->WindowThermalModel(ThermalModelNum).CalculationStandard;
2795 37990 : ThermalMod = state.dataMaterial->WindowThermalModel(ThermalModelNum).ThermalModel;
2796 37990 : CalcDeflection = state.dataMaterial->WindowThermalModel(ThermalModelNum).DeflectionModel;
2797 37990 : SDScalar = state.dataMaterial->WindowThermalModel(ThermalModelNum).SDScalar;
2798 37990 : VacuumPressure = state.dataMaterial->WindowThermalModel(ThermalModelNum).VacuumPressureLimit;
2799 37990 : Tini = state.dataMaterial->WindowThermalModel(ThermalModelNum).InitialTemperature - Constant::Kelvin;
2800 37990 : Pini = state.dataMaterial->WindowThermalModel(ThermalModelNum).InitialPressure;
2801 :
2802 37990 : nlayer = state.dataConstruction->Construct(ConstrNum).TotSolidLayers;
2803 37990 : isky = 3; // IR radiation is provided from external source
2804 37990 : iwd = 0; // assume windward for now. TODO compare surface normal with wind direction
2805 :
2806 37990 : if (CalcCondition == DataBSDFWindow::Condition::Invalid) {
2807 37968 : ZoneNum = state.dataSurface->Surface(SurfNum).Zone;
2808 37968 : Real64 RefAirTemp = state.dataSurface->Surface(SurfNum).getInsideAirTemperature(state, SurfNum);
2809 37968 : tind = RefAirTemp + Constant::Kelvin; // Inside air temperature
2810 :
2811 : // now get "outside" air temperature
2812 37968 : if (SurfNumAdj > 0) { // Interzone window
2813 :
2814 0 : int enclNumAdj = state.dataSurface->Surface(SurfNumAdj).RadEnclIndex;
2815 0 : RefAirTemp = state.dataSurface->Surface(SurfNumAdj).getInsideAirTemperature(state, SurfNumAdj);
2816 0 : tout = RefAirTemp + Constant::Kelvin; // outside air temperature
2817 :
2818 0 : tsky = state.dataViewFactor->EnclRadInfo(enclNumAdj).MRT +
2819 : Constant::Kelvin; // TODO this misses IR from sources such as high temp radiant and baseboards
2820 :
2821 : // ! Add long-wave radiation from adjacent zone absorbed by glass layer closest to the adjacent zone.
2822 : // AbsRadGlassFace(1) = AbsRadGlassFace(1) + SurfQRadThermInAbs(SurfNumAdj)
2823 : // ! The IR radiance of this window's "exterior" surround is the IR radiance
2824 : // ! from surfaces and high-temp radiant sources in the adjacent zone
2825 0 : outir = state.dataSurface->SurfWinIRfromParentZone(SurfNumAdj) + state.dataHeatBalSurf->SurfQdotRadHVACInPerArea(SurfNumAdj);
2826 :
2827 : } else { // Exterior window (ExtBoundCond = 0)
2828 : // Calculate LWR from surrounding surfaces if defined for an exterior window
2829 37968 : OutSrdIR = 0;
2830 37968 : if (state.dataGlobal->AnyLocalEnvironmentsInModel) {
2831 0 : if (state.dataSurface->Surface(SurfNum).SurfHasSurroundingSurfProperty) {
2832 0 : SrdSurfTempAbs = state.dataSurface->Surface(SurfNum).SrdSurfTemp + Constant::Kelvin;
2833 0 : OutSrdIR = Constant::StefanBoltzmann * state.dataSurface->Surface(SurfNum).ViewFactorSrdSurfs * pow_4(SrdSurfTempAbs);
2834 : }
2835 : }
2836 37968 : if (state.dataSurface->Surface(SurfNum).ExtWind) { // Window is exposed to wind (and possibly rain)
2837 37968 : if (state.dataEnvrn->IsRain) { // Raining: since wind exposed, outside window surface gets wet
2838 0 : tout = state.dataSurface->SurfOutWetBulbTemp(SurfNum) + Constant::Kelvin;
2839 : } else { // Dry
2840 37968 : tout = state.dataSurface->SurfOutDryBulbTemp(SurfNum) + Constant::Kelvin;
2841 : }
2842 : } else { // Window not exposed to wind
2843 0 : tout = state.dataSurface->SurfOutDryBulbTemp(SurfNum) + Constant::Kelvin;
2844 : }
2845 : // tsky = SkyTemp + TKelvin
2846 37968 : tsky = state.dataEnvrn->SkyTempKelvin;
2847 37968 : Ebout = state.dataWindowComplexManager->sigma * pow_4(tout);
2848 37968 : outir = state.dataSurface->Surface(SurfNum).ViewFactorSkyIR *
2849 37968 : (state.dataSurface->SurfAirSkyRadSplit(SurfNum) * state.dataWindowComplexManager->sigma * pow_4(tsky) +
2850 37968 : (1.0 - state.dataSurface->SurfAirSkyRadSplit(SurfNum)) * Ebout) +
2851 37968 : state.dataSurface->Surface(SurfNum).ViewFactorGroundIR * Ebout + OutSrdIR;
2852 : }
2853 :
2854 37968 : hin = state.dataHeatBalSurf->SurfHConvInt(SurfNum); // Room-side surface convective film conductance
2855 37968 : ibc(2) = 0; // convective coefficient on indoor side will be recalculated (like in Winkelmann routines)
2856 :
2857 : // hcout=HextConvCoeff ! Exterior convection coefficient is passed in from outer routine
2858 37968 : hout = HextConvCoeff; // Exterior convection coefficient is passed in from outer routine
2859 37968 : ibc(1) = 2; // prescribed convective film coeff on outdoor side
2860 37968 : tilt = state.dataSurface->Surface(SurfNum).Tilt;
2861 37968 : height = state.dataSurface->Surface(SurfNum).Height;
2862 37968 : heightt = height; // for now put same window and glazing pocket hights
2863 37968 : width = state.dataSurface->Surface(SurfNum).Width;
2864 :
2865 : // indoor mean radiant temperature.
2866 : // IR incident on window from zone surfaces and high-temp radiant sources
2867 37968 : rmir = state.dataSurface->SurfWinIRfromParentZone(SurfNum) + state.dataHeatBalSurf->SurfQdotRadHVACInPerArea(SurfNum);
2868 37968 : trmin = root_4(rmir / Constant::StefanBoltzmann); // TODO check model equation.
2869 :
2870 : // outdoor wind speed
2871 37968 : if (!state.dataSurface->Surface(SurfNum).ExtWind) {
2872 0 : wso = 0.0; // No wind exposure
2873 : // ELSE IF (Surface(SurfNum)%Class == SurfaceClass::Window .AND. SurfaceWindow(SurfNum)%ShadingFlag ==
2874 : // WinShadingType::ExtShade) THEN
2875 : // wso = 0.0 ! Assume zero wind speed at outside glass surface of window with exterior shade
2876 : } else {
2877 37968 : wso = state.dataSurface->SurfOutWindSpeed(SurfNum);
2878 : }
2879 :
2880 : // indoor wind speed
2881 37968 : wsi = 0.0; // assumuption (TODO, what to use for inside air velocity?)
2882 :
2883 37968 : fclr = 1.0 - state.dataEnvrn->CloudFraction;
2884 : }
2885 :
2886 37990 : TotLay = state.dataConstruction->Construct(ConstrNum).TotLayers;
2887 37990 : IGap = 0;
2888 :
2889 : //****************************************************************************************************
2890 : // Inside and outside gas coefficients
2891 : //****************************************************************************************************
2892 37990 : iprop(1, 1) = 1; // air on outdoor side
2893 37990 : frct(1, 1) = 1.0; // pure air on outdoor side
2894 37990 : nmix(1) = 1; // pure air on outdoor side
2895 :
2896 37990 : iprop(1, nlayer + 1) = 1; // air on indoor side
2897 37990 : frct(1, nlayer + 1) = 1.0; // pure air on indoor side
2898 37990 : nmix(nlayer + 1) = 1; // pure air on indoor side
2899 :
2900 : // Simon: feed gas coefficients with air. This is necessary for tarcog because it is used on indoor and outdoor sides
2901 37990 : auto const &gasAir = Material::gases[(int)Material::GasType::Air];
2902 37990 : wght(iprop(1, 1)) = gasAir.wght;
2903 37990 : gama(iprop(1, 1)) = gasAir.specHeatRatio;
2904 :
2905 37990 : gcon(1, iprop(1, 1)) = gasAir.con.c0;
2906 37990 : gcon(2, iprop(1, 1)) = gasAir.con.c1;
2907 37990 : gcon(3, iprop(1, 1)) = gasAir.con.c2;
2908 37990 : gvis(1, iprop(1, 1)) = gasAir.vis.c0;
2909 37990 : gvis(2, iprop(1, 1)) = gasAir.vis.c1;
2910 37990 : gvis(3, iprop(1, 1)) = gasAir.vis.c2;
2911 37990 : gcp(1, iprop(1, 1)) = gasAir.cp.c0;
2912 37990 : gcp(2, iprop(1, 1)) = gasAir.cp.c1;
2913 37990 : gcp(3, iprop(1, 1)) = gasAir.cp.c2;
2914 :
2915 : // Fill window layer properties needed for window layer heat balance calculation
2916 37990 : IGlass = 0;
2917 37990 : IGap = 0;
2918 213144 : for (Lay = 1; Lay <= TotLay; ++Lay) {
2919 175154 : LayPtr = state.dataConstruction->Construct(ConstrNum).LayerPoint(Lay);
2920 175154 : auto const *mat = s_mat->materials(LayPtr);
2921 :
2922 175154 : if ((mat->group == Material::Group::Glass) || (mat->group == Material::Group::GlassSimple)) {
2923 76654 : auto const *matGlass = dynamic_cast<Material::MaterialGlass const *>(mat);
2924 76654 : assert(matGlass != nullptr);
2925 :
2926 76654 : ++IGlass;
2927 76654 : LayerType(IGlass) = TARCOGParams::TARCOGLayerType::SPECULAR; // this marks specular layer type
2928 76654 : thick(IGlass) = matGlass->Thickness;
2929 76654 : scon(IGlass) = matGlass->Conductivity;
2930 76654 : emis(2 * IGlass - 1) = matGlass->AbsorpThermalFront;
2931 76654 : emis(2 * IGlass) = matGlass->AbsorpThermalBack;
2932 76654 : tir(2 * IGlass - 1) = matGlass->TransThermal;
2933 76654 : tir(2 * IGlass) = matGlass->TransThermal;
2934 76654 : YoungsMod(IGlass) = matGlass->YoungModulus;
2935 76654 : PoissonsRat(IGlass) = matGlass->PoissonsRatio;
2936 :
2937 175154 : } else if (mat->group == Material::Group::ComplexShade) {
2938 29918 : auto const *matShade = dynamic_cast<Material::MaterialComplexShade const *>(mat);
2939 29918 : ++IGlass;
2940 29918 : LayerType(IGlass) = matShade->LayerType;
2941 :
2942 29918 : thick(IGlass) = matShade->Thickness;
2943 29918 : scon(IGlass) = matShade->Conductivity;
2944 29918 : emis(2 * IGlass - 1) = matShade->FrontEmissivity;
2945 29918 : emis(2 * IGlass) = matShade->BackEmissivity;
2946 29918 : tir(2 * IGlass - 1) = matShade->TransThermal;
2947 29918 : tir(2 * IGlass) = matShade->TransThermal;
2948 :
2949 : // This needs to be converted into correct areas. That can be done only after loading complete window data
2950 29918 : Atop(IGlass) = matShade->topOpeningMult;
2951 29918 : Abot(IGlass) = matShade->bottomOpeningMult;
2952 29918 : Al(IGlass) = matShade->leftOpeningMult;
2953 29918 : Ar(IGlass) = matShade->rightOpeningMult;
2954 29918 : Ah(IGlass) = matShade->frontOpeningMult;
2955 :
2956 29918 : SlatThick(IGlass) = matShade->SlatThickness;
2957 29918 : SlatWidth(IGlass) = matShade->SlatWidth;
2958 29918 : SlatAngle(IGlass) = matShade->SlatAngle;
2959 29918 : SlatCond(IGlass) = matShade->SlatConductivity;
2960 29918 : SlatSpacing(IGlass) = matShade->SlatSpacing;
2961 29918 : SlatCurve(IGlass) = matShade->SlatCurve;
2962 :
2963 68582 : } else if (mat->group == Material::Group::ComplexWindowGap) {
2964 68582 : auto const *matGap = dynamic_cast<Material::MaterialComplexWindowGap const *>(mat);
2965 68582 : ++IGap;
2966 68582 : gap(IGap) = matGap->Thickness;
2967 68582 : presure(IGap) = matGap->Pressure;
2968 :
2969 68582 : GapDefMax(IGap) = matGap->deflectedThickness;
2970 :
2971 68582 : PillarSpacing(IGap) = matGap->pillarSpacing;
2972 68582 : PillarRadius(IGap) = matGap->pillarRadius;
2973 68582 : SupportPlr(IGap) = matGap->pillarSpacing != 0.0 && matGap->pillarRadius != 0.0;
2974 :
2975 68582 : nmix(IGap + 1) = matGap->numGases;
2976 149272 : for (IMix = 1; IMix <= nmix(IGap + 1); ++IMix) {
2977 80690 : auto const &gas = matGap->gases[IMix - 1];
2978 80690 : frct(IMix, IGap + 1) = matGap->gasFracts[IMix - 1];
2979 :
2980 : // Now has to build-up gas coefficients arrays. All used gasses should be stored into these arrays and
2981 : // to be correctly referenced by gap arrays
2982 :
2983 : // First check if gas coefficients are already part of array. Duplicates are not necessary
2984 80690 : bool feedData(false);
2985 80690 : CheckGasCoefs(gas.wght, iprop(IMix, IGap + 1), wght, feedData);
2986 80690 : if (feedData) {
2987 2 : wght(iprop(IMix, IGap + 1)) = gas.wght;
2988 2 : gama(iprop(IMix, IGap + 1)) = gas.specHeatRatio;
2989 :
2990 2 : gcon(1, iprop(IMix, IGap + 1)) = gas.con.c0;
2991 2 : gcon(2, iprop(IMix, IGap + 1)) = gas.con.c1;
2992 2 : gcon(3, iprop(IMix, IGap + 1)) = gas.con.c2;
2993 :
2994 2 : gvis(1, iprop(IMix, IGap + 1)) = gas.vis.c0;
2995 2 : gvis(2, iprop(IMix, IGap + 1)) = gas.vis.c1;
2996 2 : gvis(3, iprop(IMix, IGap + 1)) = gas.vis.c2;
2997 :
2998 2 : gcp(1, iprop(IMix, IGap + 1)) = gas.cp.c0;
2999 2 : gcp(2, iprop(IMix, IGap + 1)) = gas.cp.c1;
3000 2 : gcp(3, iprop(IMix, IGap + 1)) = gas.cp.c2;
3001 : } // IF feedData THEN
3002 : }
3003 : } else {
3004 0 : ShowContinueError(state, "Illegal layer type in Construction:ComplexFenestrationState.");
3005 0 : ShowContinueError(state, "Allowed object are:");
3006 0 : ShowContinueError(state, " - WindowMaterial:Glazing");
3007 0 : ShowContinueError(state, " - WindowMaterial:ComplexShade");
3008 0 : ShowContinueError(state, " - WindowMaterial:Gap");
3009 0 : ShowFatalError(state, "halting because of error in layer definition for Construction:ComplexFenestrationState");
3010 : }
3011 :
3012 : } // End of loop over glass, gap and blind/shade layers in a window construction
3013 :
3014 37990 : if (CalcCondition == DataBSDFWindow::Condition::Invalid) {
3015 : // now calculate correct areas for multipliers
3016 144480 : for (Lay = 1; Lay <= nlayer; ++Lay) {
3017 106512 : if (LayerType(Lay) != TARCOGParams::TARCOGLayerType::SPECULAR) { // Layer is shading
3018 : // before changing multipliers, need to determine which one is dominant gap width
3019 29904 : if (Lay == 1) { // Exterior shading device
3020 10080 : dominantGapWidth = gap(Lay);
3021 19824 : } else if (Lay == nlayer) { // Interior shading device
3022 15792 : dominantGapWidth = gap(Lay - 1);
3023 : } else { // In-between shading device
3024 4032 : dominantGapWidth = min(gap(Lay - 1), gap(Lay));
3025 : }
3026 29904 : Atop(Lay) *= dominantGapWidth * width;
3027 29904 : Abot(Lay) *= dominantGapWidth * width;
3028 29904 : Al(Lay) *= dominantGapWidth * height;
3029 29904 : Ar(Lay) *= dominantGapWidth * height;
3030 29904 : Ah(Lay) *= width * height;
3031 : }
3032 : }
3033 : }
3034 :
3035 : // ThermalMod = 0
3036 37990 : CalcSHGC = 0;
3037 :
3038 37990 : Window_ID = ConstrNum;
3039 :
3040 : // vector of absorbed solar energy fractions for each layer.
3041 37990 : asol = 0.0;
3042 : // direct solar radiation
3043 37990 : if (CalcCondition == DataBSDFWindow::Condition::Invalid) {
3044 37968 : ShadeFlag = state.dataSurface->SurfWinShadingFlag(SurfNum);
3045 37968 : dir = state.dataHeatBal->SurfQRadSWOutIncident(SurfNum) +
3046 37968 : state.dataHeatBal->EnclSolQSWRad(state.dataSurface->Surface(SurfNum).SolarEnclIndex); // TODO, check , !
3047 : // currently using Exterior beam plus diffuse solar incident on surface
3048 : // plus zone short wave. CHECK
3049 : // if (dir.ne.0.0d0) then
3050 144480 : for (IGlass = 1; IGlass <= nlayer; ++IGlass) {
3051 : // IF (dir > 0.0D0 ) THEN
3052 106512 : asol(IGlass) = state.dataHeatBal->SurfWinQRadSWwinAbs(SurfNum, IGlass);
3053 : // ELSE
3054 : // asol(IGLASS) = 0.0D0
3055 : // ENDIF
3056 : }
3057 : // end if
3058 :
3059 : // Add contribution of IR from zone internal gains (lights, equipment and people). This is absorbed in zone-side layer and it
3060 : // is assumed that nothing is transmitted through
3061 37968 : asol(nlayer) += state.dataHeatBal->SurfQdotRadIntGainsInPerArea(SurfNum);
3062 :
3063 37968 : presure = state.dataEnvrn->OutBaroPress;
3064 :
3065 : // Instead of doing temperature guess get solution from previous iteration. That should be much better than guess
3066 250992 : for (k = 1; k <= 2 * nlayer; ++k) {
3067 213024 : theta(k) = state.dataSurface->SurfaceWindow(SurfNum).thetaFace[k];
3068 : }
3069 : }
3070 :
3071 : // Standard conditions run (winter and summer)
3072 37990 : if (CalcCondition == DataBSDFWindow::Condition::Winter) {
3073 11 : tind = 294.15;
3074 11 : tout = 255.15;
3075 11 : hcout = 26.0;
3076 11 : wso = 5.5;
3077 11 : dir = 0.0;
3078 37979 : } else if (CalcCondition == DataBSDFWindow::Condition::Summer) {
3079 11 : tind = 297.15;
3080 11 : tout = 305.15;
3081 11 : hcout = 15.0;
3082 11 : wso = 2.75;
3083 11 : dir = 783.0;
3084 11 : CalcSHGC = 1;
3085 : }
3086 :
3087 : // Common condition data
3088 37990 : if (CalcCondition != DataBSDFWindow::Condition::Invalid) {
3089 22 : trmin = tind;
3090 22 : outir = 0.0;
3091 22 : tsky = tout;
3092 22 : wsi = 0.0;
3093 22 : fclr = 1.0;
3094 22 : ibc(1) = 0;
3095 22 : ibc(2) = 0;
3096 22 : presure = 101325.0;
3097 22 : iwd = 0; // Windward wind direction
3098 22 : isky = 0;
3099 22 : esky = 1.0;
3100 22 : height = 1.0;
3101 22 : heightt = 1.0;
3102 22 : width = 1.0;
3103 22 : tilt = 90.0;
3104 : // Just to make initial quess different from absolute zero
3105 22 : theta = 273.15;
3106 : }
3107 :
3108 37990 : edgeGlCorrFac = (SurfNum != 0) ? state.dataSurface->SurfaceWindow(SurfNum).edgeGlassCorrFac : 1;
3109 :
3110 : // call TARCOG
3111 37990 : int constexpr Debug_mode = 0;
3112 37990 : TARCOG90(state,
3113 : nlayer,
3114 : iwd,
3115 : tout,
3116 : tind,
3117 : trmin,
3118 : wso,
3119 : wsi,
3120 : dir,
3121 : outir,
3122 : isky,
3123 : tsky,
3124 : esky,
3125 : fclr,
3126 : VacuumPressure,
3127 : VacuumMaxGapThickness,
3128 : CalcDeflection,
3129 : Pa,
3130 : Pini,
3131 : Tini,
3132 : gap,
3133 : GapDefMax,
3134 : thick,
3135 : scon,
3136 : YoungsMod,
3137 : PoissonsRat,
3138 : tir,
3139 : emis,
3140 : totsol,
3141 : tilt,
3142 : asol,
3143 : height,
3144 : heightt,
3145 : width,
3146 : presure,
3147 : iprop,
3148 : frct,
3149 : gcon,
3150 : gvis,
3151 : gcp,
3152 : wght,
3153 : gama,
3154 : nmix,
3155 : SupportPlr,
3156 : PillarSpacing,
3157 : PillarRadius,
3158 : theta,
3159 : LayerDef,
3160 : q,
3161 : qv,
3162 : ufactor,
3163 : sc,
3164 : hflux,
3165 : hcin,
3166 : hcout,
3167 : hrin,
3168 : hrout,
3169 : hin,
3170 : hout,
3171 : hcgap,
3172 : hrgap,
3173 : shgc,
3174 : nperr,
3175 : tarcogErrorMessage,
3176 : shgct,
3177 : tamb,
3178 : troom,
3179 : ibc,
3180 : Atop,
3181 : Abot,
3182 : Al,
3183 : Ar,
3184 : Ah,
3185 : SlatThick,
3186 : SlatWidth,
3187 : SlatAngle,
3188 : SlatCond,
3189 : SlatSpacing,
3190 : SlatCurve,
3191 : vvent,
3192 : tvent,
3193 : LayerType,
3194 : nslice,
3195 : LaminateA,
3196 : LaminateB,
3197 : sumsol,
3198 : hg,
3199 : hr,
3200 : hs,
3201 : he,
3202 : hi,
3203 : Ra,
3204 : Nu,
3205 : standard,
3206 : ThermalMod,
3207 : Debug_mode,
3208 : Debug_dir,
3209 : Debug_file,
3210 : Window_ID,
3211 : IGU_ID,
3212 : ShadeEmisRatioOut,
3213 : ShadeEmisRatioIn,
3214 : ShadeHcRatioOut,
3215 : ShadeHcRatioIn,
3216 : HcUnshadedOut,
3217 : HcUnshadedIn,
3218 : Keff,
3219 : ShadeGapKeffConv,
3220 : SDScalar,
3221 : CalcSHGC,
3222 : NumOfIterations,
3223 : edgeGlCorrFac);
3224 :
3225 : // process results from TARCOG
3226 37990 : if ((nperr > 0) && (nperr < 1000)) { // process error signal from tarcog
3227 :
3228 0 : ShowSevereError(state, "Window tarcog returned an error");
3229 0 : tarcogErrorMessage = "message = \"" + tarcogErrorMessage + "\"";
3230 0 : ShowContinueErrorTimeStamp(state, tarcogErrorMessage);
3231 0 : if (CalcCondition == DataBSDFWindow::Condition::Invalid) {
3232 0 : ShowContinueError(state, format("surface name = {}", state.dataSurface->Surface(SurfNum).Name));
3233 : }
3234 0 : ShowContinueError(state, format("construction name = {}", state.dataConstruction->Construct(ConstrNum).Name));
3235 0 : ShowFatalError(state, "halting because of error in tarcog");
3236 37990 : } else if (CalcCondition == DataBSDFWindow::Condition::Winter) {
3237 11 : state.dataHeatBal->NominalU(ConstrNum) = ufactor;
3238 37979 : } else if (CalcCondition == DataBSDFWindow::Condition::Summer) {
3239 : // tempInt = SurfaceWindow(SurfNum)%ComplexFen%CurrentState
3240 : // tempReal = SurfaceWindow(SurfNum)%ComplexFen%State(tempInt)%WinDiffTrans
3241 :
3242 : // Sum1 = 0.0d0
3243 : // Sum2 = 0.0d0
3244 : // do j = 1 , ComplexWind(SurfNum)%Geom%Inc%NBasis !Incident ray loop
3245 : // Sum2 = Sum2 + ComplexWind(SurfNum)%Geom%Inc%Lamda (j)
3246 : // do m = 1 , ComplexWind(SurfNum)%Geom%Trn%NBasis !Outgoing ray loop
3247 : // Sum1 =Sum1 + ComplexWind(SurfNum)%Geom%Inc%Lamda(j) * ComplexWind(SurfNum)%Geom%Trn%Lamda(m) * &
3248 : // & Construct(ConstrNum)%BSDFInput%SolFrtTrans ( j , m )
3249 : // end do !Outgoing ray loop
3250 : // end do !Incident ray loop
3251 : // if (Sum2 > 0.0d0) THEN
3252 : // tempReal = Sum1/Sum2
3253 : // else
3254 : // tempReal = 0.0d0
3255 : // CALL ShowWarningError(state, 'BSDF--Inc basis has zero projected solid angle')
3256 : // endif
3257 :
3258 11 : state.dataConstruction->Construct(ConstrNum).SummerSHGC = shgc;
3259 :
3260 : // Construct(SurfNum)%VisTransNorm = SurfaceWindow(SurfNum)%ComplexFen%State(tempInt)%WinDiffVisTrans
3261 37968 : } else if (CalcCondition == DataBSDFWindow::Condition::Invalid) { // expect converged results...
3262 : // Window heat balance solution has converged.
3263 :
3264 37968 : state.dataSurface->SurfWinWindowCalcIterationsRep(SurfNum) = NumOfIterations;
3265 37968 : state.dataHeatBalSurf->SurfHConvInt(SurfNum) = hcin;
3266 :
3267 : // For interior shade, add convective gain from glass/shade gap air flow to zone convective gain;
3268 : // For all cases, get total window heat gain for reporting. See CalcWinFrameAndDividerTemps for
3269 : // contribution of frame and divider.
3270 :
3271 37968 : SurfInsideTemp = theta(2 * nlayer) - Constant::Kelvin;
3272 37968 : state.dataSurface->SurfWinEffInsSurfTemp(SurfNum) = SurfInsideTemp;
3273 37968 : SurfOutsideTemp = theta(1) - Constant::Kelvin;
3274 37968 : SurfOutsideEmiss = emis(1);
3275 :
3276 37968 : IncidentSolar = state.dataSurface->Surface(SurfNum).Area * state.dataHeatBal->SurfQRadSWOutIncident(SurfNum);
3277 37968 : if (ANY_INTERIOR_SHADE_BLIND(ShadeFlag)) {
3278 : // Interior shade or blind
3279 15792 : ConvHeatFlowNatural = -qv(nlayer) * height * width;
3280 :
3281 15792 : state.dataSurface->SurfWinConvHeatFlowNatural(SurfNum) = ConvHeatFlowNatural;
3282 15792 : state.dataSurface->SurfWinGapConvHtFlowRep(SurfNum) = ConvHeatFlowNatural;
3283 15792 : state.dataSurface->SurfWinGapConvHtFlowRepEnergy(SurfNum) =
3284 15792 : state.dataSurface->SurfWinGapConvHtFlowRep(SurfNum) * state.dataGlobal->TimeStepZoneSec;
3285 : // Window heat gain from glazing and shade/blind to zone. Consists of transmitted solar, convection
3286 : // from air exiting gap, convection from zone-side of shade/blind, net IR to zone from shade and net IR to
3287 : // zone from the glass adjacent to the shade/blind (zero if shade/blind IR transmittance is zero).
3288 : // Following assumes glazed area = window area (i.e., dividers ignored) in calculating
3289 : // IR to zone from glass when interior shade/blind is present.
3290 15792 : ShadeArea = state.dataSurface->Surface(SurfNum).Area + state.dataSurface->SurfWinDividerArea(SurfNum);
3291 15792 : sconsh = scon(ngllayer + 1) / thick(ngllayer + 1);
3292 15792 : nglfacep = nglface + 2;
3293 15792 : CondHeatGainShade = ShadeArea * sconsh * (theta(nglfacep - 1) - theta(nglfacep));
3294 15792 : EpsShIR1 = emis(nglface + 1);
3295 15792 : EpsShIR2 = emis(nglface + 2);
3296 15792 : TauShIR = tir(nglface + 1);
3297 15792 : RhoShIR1 = max(0.0, 1.0 - TauShIR - EpsShIR1);
3298 15792 : RhoShIR2 = max(0.0, 1.0 - TauShIR - EpsShIR2);
3299 15792 : RhoGlIR2 = 1.0 - emis(2 * ngllayer);
3300 15792 : ShGlReflFacIR = 1.0 - RhoGlIR2 * RhoShIR1;
3301 15792 : NetIRHeatGainShade =
3302 15792 : ShadeArea * EpsShIR2 * (state.dataWindowComplexManager->sigma * pow_4(theta(nglfacep)) - rmir) +
3303 15792 : EpsShIR1 * (state.dataWindowComplexManager->sigma * pow_4(theta(nglfacep - 1)) - rmir) * RhoGlIR2 * TauShIR / ShGlReflFacIR;
3304 15792 : NetIRHeatGainGlass = ShadeArea * (emis(2 * ngllayer) * TauShIR / ShGlReflFacIR) *
3305 15792 : (state.dataWindowComplexManager->sigma * pow_4(theta(2 * ngllayer)) - rmir);
3306 15792 : ConvHeatGainFrZoneSideOfShade = ShadeArea * hcin * (theta(nglfacep) - tind);
3307 31584 : state.dataSurface->SurfWinHeatGain(SurfNum) = state.dataSurface->SurfWinTransSolar(SurfNum) + ConvHeatFlowNatural +
3308 15792 : ConvHeatGainFrZoneSideOfShade + NetIRHeatGainGlass + NetIRHeatGainShade;
3309 : // store components for reporting
3310 15792 : state.dataSurface->SurfWinGainConvShadeToZoneRep(SurfNum) = ConvHeatGainFrZoneSideOfShade;
3311 15792 : state.dataSurface->SurfWinGainIRGlazToZoneRep(SurfNum) = NetIRHeatGainGlass;
3312 15792 : state.dataSurface->SurfWinGainIRShadeToZoneRep(SurfNum) = NetIRHeatGainShade;
3313 : } else {
3314 : // Interior shade or blind not present; innermost layer is glass
3315 22176 : CondHeatGainGlass =
3316 22176 : state.dataSurface->Surface(SurfNum).Area * scon(nlayer) / thick(nlayer) * (theta(2 * nlayer - 1) - theta(2 * nlayer));
3317 22176 : NetIRHeatGainGlass = state.dataSurface->Surface(SurfNum).Area * emis(2 * nlayer) *
3318 22176 : (state.dataWindowComplexManager->sigma * pow_4(theta(2 * nlayer)) - rmir);
3319 22176 : ConvHeatGainFrZoneSideOfGlass = state.dataSurface->Surface(SurfNum).Area * hcin * (theta(2 * nlayer) - tind);
3320 22176 : state.dataSurface->SurfWinHeatGain(SurfNum) =
3321 22176 : state.dataSurface->SurfWinTransSolar(SurfNum) + ConvHeatGainFrZoneSideOfGlass + NetIRHeatGainGlass;
3322 : // store components for reporting
3323 22176 : state.dataSurface->SurfWinGainConvGlazToZoneRep(SurfNum) = ConvHeatGainFrZoneSideOfGlass;
3324 22176 : state.dataSurface->SurfWinGainIRGlazToZoneRep(SurfNum) = NetIRHeatGainGlass;
3325 : // need to report convective heat flow from the gap in case of exterior shade
3326 22176 : if (ShadeFlag == WinShadingType::ExtShade) {
3327 10080 : ConvHeatFlowNatural = -qv(2) * height * width; // qv(1) is exterior environment
3328 :
3329 10080 : state.dataSurface->SurfWinGapConvHtFlowRep(SurfNum) = ConvHeatFlowNatural;
3330 10080 : state.dataSurface->SurfWinGapConvHtFlowRepEnergy(SurfNum) =
3331 10080 : state.dataSurface->SurfWinGapConvHtFlowRep(SurfNum) * state.dataGlobal->TimeStepZoneSec;
3332 : }
3333 : }
3334 :
3335 : // Add convective heat gain from airflow window
3336 : // Note: effect of fan heat on gap outlet temperature is neglected since fan power (based
3337 : // on pressure drop through the gap) is extremely small
3338 :
3339 : // WinGapConvHtFlowRep(SurfNum) = 0.0d0
3340 : // WinGapConvHtFlowRepEnergy(SurfNum) = 0.0d0
3341 37968 : TotAirflowGap = state.dataSurface->SurfWinAirflowThisTS(SurfNum) * state.dataSurface->Surface(SurfNum).Width;
3342 37968 : TAirflowGapOutlet = Constant::Kelvin; // TODO Need to calculate this
3343 37968 : TAirflowGapOutletC = TAirflowGapOutlet - Constant::Kelvin;
3344 37968 : state.dataSurface->SurfWinTAirflowGapOutlet(SurfNum) = TAirflowGapOutletC;
3345 37968 : if (state.dataSurface->SurfWinAirflowThisTS(SurfNum) > 0.0) {
3346 0 : ConvHeatFlowForced = sum(qv); // TODO. figure forced ventilation heat flow in Watts
3347 :
3348 0 : state.dataSurface->SurfWinGapConvHtFlowRep(SurfNum) = ConvHeatFlowForced;
3349 0 : state.dataSurface->SurfWinGapConvHtFlowRepEnergy(SurfNum) =
3350 0 : state.dataSurface->SurfWinGapConvHtFlowRep(SurfNum) * state.dataGlobal->TimeStepZoneSec;
3351 : // Add heat from gap airflow to zone air if destination is inside air; save the heat gain to return
3352 : // air in case it needs to be sent to the zone (due to no return air determined in HVAC simulation)
3353 0 : if (state.dataSurface->SurfWinAirflowDestination(SurfNum) == WindowAirFlowDestination::Indoor ||
3354 0 : state.dataSurface->SurfWinAirflowDestination(SurfNum) == WindowAirFlowDestination::Return) {
3355 0 : auto &thisZoneHB = state.dataZoneTempPredictorCorrector->zoneHeatBalance(ZoneNum);
3356 0 : if (state.dataSurface->SurfWinAirflowSource(SurfNum) == WindowAirFlowSource::Indoor) {
3357 0 : InletAirHumRat = thisZoneHB.airHumRat;
3358 : } else { // AirflowSource = outside air
3359 0 : InletAirHumRat = state.dataEnvrn->OutHumRat;
3360 : }
3361 0 : ZoneTemp = thisZoneHB.MAT; // this should be Tin (account for different reference temps)
3362 0 : CpAirOutlet = PsyCpAirFnW(InletAirHumRat);
3363 0 : CpAirZone = PsyCpAirFnW(thisZoneHB.airHumRat);
3364 0 : ConvHeatGainToZoneAir = TotAirflowGap * (CpAirOutlet * (TAirflowGapOutletC)-CpAirZone * ZoneTemp);
3365 0 : if (state.dataSurface->SurfWinAirflowDestination(SurfNum) == WindowAirFlowDestination::Indoor) {
3366 0 : state.dataSurface->SurfWinConvHeatGainToZoneAir(SurfNum) = ConvHeatGainToZoneAir;
3367 0 : state.dataSurface->SurfWinHeatGain(SurfNum) += ConvHeatGainToZoneAir;
3368 : } else {
3369 0 : state.dataSurface->SurfWinRetHeatGainToZoneAir(SurfNum) = ConvHeatGainToZoneAir;
3370 : }
3371 : }
3372 : // For AirflowDestination = ReturnAir in a controlled (i.e., conditioned) zone with return air, see CalcZoneLeavingConditions
3373 : // for calculation of modification of return-air temperature due to airflow from window gaps into return air.
3374 : }
3375 :
3376 : // Correct WinHeatGain for interior diffuse shortwave (solar and shortwave from lights) transmitted
3377 : // back out window
3378 37968 : ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
3379 : // ConstrNumSh = Surface(SurfNum)%ShadedConstruction
3380 : // IF(SurfaceWindow(SurfNum)%StormWinFlag==1) THEN
3381 : // ConstrNum = Surface(SurfNum)%StormWinConstruction
3382 : // ConstrNumSh = Surface(SurfNum)%StormWinShadedConstruction
3383 : // END IF
3384 : // IF(ShadeFlag <= 0) THEN
3385 : // TransDiff = Construct(ConstrNum).TransDiff;
3386 37968 : int IState = state.dataSurface->SurfaceWindow(SurfNum).ComplexFen.NumStates;
3387 37968 : Real64 ReflDiff = state.dataSurface->SurfaceWindow(SurfNum).ComplexFen.State(IState).WinBkHemRefl;
3388 : // ELSE IF(ShadeFlag==WinShadingType::IntShade .OR. ShadeFlag==WinShadingType::ExtShade) THEN
3389 : // TransDiff = Construct(ConstrNum)%TransDiff
3390 : // ELSE IF(ShadeFlag==WinShadingType::IntBlind .OR. ShadeFlag==WinShadingType::ExtBlind .OR.ShadeFlag==WinShadingType::BGBlind) THEN
3391 : // TransDiff = InterpSlatAng(SurfaceWindow(SurfNum)%SlatAngThisTS,SurfaceWindow(SurfNum)%MovableSlats, &
3392 : // Construct(ConstrNumSh)%BlTransDiff)
3393 : // ELSE IF(ShadeFlag == SwitchableGlazing) THEN
3394 : // TransDiff = InterpSW(SurfaceWindow(SurfNum)%SwitchingFactor,Construct(ConstrNum)%TransDiff, &
3395 : // Construct(ConstrNumSh)%TransDiff)
3396 : // END IF
3397 37968 : state.dataSurface->SurfWinLossSWZoneToOutWinRep(SurfNum) =
3398 37968 : state.dataHeatBal->EnclSolQSWRad(state.dataSurface->Surface(SurfNum).SolarEnclIndex) * state.dataSurface->Surface(SurfNum).Area *
3399 75936 : (1 - ReflDiff) +
3400 37968 : state.dataHeatBalSurf->SurfWinInitialBeamSolInTrans(SurfNum);
3401 37968 : state.dataSurface->SurfWinHeatGain(SurfNum) -=
3402 37968 : (state.dataSurface->SurfWinLossSWZoneToOutWinRep(SurfNum) +
3403 37968 : state.dataHeatBalSurf->SurfWinInitialDifSolInTrans(SurfNum) * state.dataSurface->Surface(SurfNum).Area);
3404 :
3405 37968 : if (ShadeFlag == WinShadingType::IntShade || ShadeFlag == WinShadingType::ExtShade) {
3406 25872 : state.dataSurface->SurfWinShadingAbsorbedSolar(SurfNum) =
3407 25872 : (state.dataSurface->SurfWinExtBeamAbsByShade(SurfNum) + state.dataSurface->SurfWinExtDiffAbsByShade(SurfNum)) *
3408 25872 : (state.dataSurface->Surface(SurfNum).Area + state.dataSurface->SurfWinDividerArea(SurfNum));
3409 25872 : state.dataSurface->SurfWinShadingAbsorbedSolarEnergy(SurfNum) =
3410 25872 : state.dataSurface->SurfWinShadingAbsorbedSolar(SurfNum) * state.dataGlobal->TimeStepZoneSec;
3411 : }
3412 37968 : if (state.dataEnvrn->SunIsUp) {
3413 19152 : state.dataSurface->SurfWinSysSolTransmittance(SurfNum) =
3414 19152 : state.dataSurface->SurfWinTransSolar(SurfNum) /
3415 19152 : (state.dataHeatBal->SurfQRadSWOutIncident(SurfNum) *
3416 19152 : (state.dataSurface->Surface(SurfNum).Area + state.dataSurface->SurfWinDividerArea(SurfNum)) +
3417 : 0.0001);
3418 19152 : state.dataSurface->SurfWinSysSolAbsorptance(SurfNum) =
3419 19152 : (state.dataHeatBal->SurfWinQRadSWwinAbsTot(SurfNum) + state.dataSurface->SurfWinShadingAbsorbedSolar(SurfNum)) /
3420 19152 : (state.dataHeatBal->SurfQRadSWOutIncident(SurfNum) *
3421 19152 : (state.dataSurface->Surface(SurfNum).Area + state.dataSurface->SurfWinDividerArea(SurfNum)) +
3422 : 0.0001);
3423 19152 : state.dataSurface->SurfWinSysSolReflectance(SurfNum) =
3424 19152 : 1.0 - state.dataSurface->SurfWinSysSolTransmittance(SurfNum) - state.dataSurface->SurfWinSysSolAbsorptance(SurfNum);
3425 : } else {
3426 18816 : state.dataSurface->SurfWinSysSolTransmittance(SurfNum) = 0.0;
3427 18816 : state.dataSurface->SurfWinSysSolAbsorptance(SurfNum) = 0.0;
3428 18816 : state.dataSurface->SurfWinSysSolReflectance(SurfNum) = 0.0;
3429 : }
3430 :
3431 : // Save hcv for use in divider calc with interior or exterior shade (see CalcWinFrameAndDividerTemps)
3432 37968 : if (ShadeFlag == WinShadingType::IntShade) {
3433 15792 : state.dataSurface->SurfWinConvCoeffWithShade(SurfNum) = 0.0;
3434 : }
3435 :
3436 37968 : if (ShadeFlag == WinShadingType::IntShade) {
3437 15792 : auto const &surfShade = state.dataSurface->surfShades(SurfNum);
3438 15792 : SurfInsideTemp = theta(2 * ngllayer + 2) - Constant::Kelvin;
3439 :
3440 : // // Get properties of inside shading layer
3441 :
3442 15792 : Real64 EffShEmiss = surfShade.effShadeEmi;
3443 15792 : Real64 EffGlEmiss = surfShade.effGlassEmi;
3444 15792 : state.dataSurface->SurfWinEffInsSurfTemp(SurfNum) =
3445 15792 : (EffShEmiss * SurfInsideTemp + EffGlEmiss * (theta(2 * ngllayer) - Constant::Kelvin)) / (EffShEmiss + EffGlEmiss);
3446 :
3447 : } else {
3448 22176 : SurfOutsideTemp = theta(1) - Constant::Kelvin;
3449 : }
3450 :
3451 144480 : for (k = 1; k <= nlayer; ++k) {
3452 106512 : state.dataSurface->SurfaceWindow(SurfNum).thetaFace[2 * k - 1] = theta(2 * k - 1);
3453 106512 : state.dataSurface->SurfaceWindow(SurfNum).thetaFace[2 * k] = theta(2 * k);
3454 :
3455 : // temperatures for reporting
3456 106512 : state.dataHeatBal->SurfWinFenLaySurfTempFront(SurfNum, k) = theta(2 * k - 1) - Constant::Kelvin;
3457 106512 : state.dataHeatBal->SurfWinFenLaySurfTempBack(SurfNum, k) = theta(2 * k) - Constant::Kelvin;
3458 : // thetas(k) = theta(k)
3459 : }
3460 : }
3461 37990 : }
3462 :
3463 : // This function check if gas with molecular weight has already been feed into coefficients and
3464 : // feed arrays
3465 :
3466 80690 : void CheckGasCoefs(Real64 const currentWeight, int &indexNumber, Array1D<Real64> &wght, bool &feedData)
3467 : {
3468 :
3469 : // Using/Aliasing
3470 : using TARCOGGassesParams::maxgas;
3471 :
3472 : // Argument array dimensioning
3473 80690 : EP_SIZE_CHECK(wght, maxgas);
3474 :
3475 80690 : feedData = false;
3476 80690 : bool coeffFound = false;
3477 80690 : int counter = 1;
3478 173486 : while ((counter <= maxgas) && (wght(counter) != 0) && (!coeffFound)) {
3479 92796 : if (std::abs(currentWeight - wght(counter)) < 1.0e-5) {
3480 80688 : coeffFound = true;
3481 : } else {
3482 12108 : ++counter;
3483 : }
3484 : } // DO WHILE((counter.LE.maxgas).AND.(wght(couner).NE.0).AND.(.NOT.coeffFound))
3485 :
3486 : // In case coefficient is not found data needs to be stored in gas coefficients arrays
3487 80690 : if ((!coeffFound) && (counter < maxgas)) {
3488 2 : feedData = true;
3489 : }
3490 :
3491 80690 : indexNumber = counter;
3492 80690 : }
3493 :
3494 3350 : int SearchAscTable(Real64 const y, // Value to be found in the table
3495 : int const n, // Number of values in the table
3496 : Array1S<Real64> const ytab // Table of values, monotonic, ascending order
3497 : )
3498 : {
3499 :
3500 : // FUNCTION INFORMATION:
3501 : // AUTHOR Joe Klems
3502 : // DATE WRITTEN Feb 2011
3503 : // MODIFIED na
3504 : // RE-ENGINEERED na
3505 :
3506 : // PURPOSE OF THIS FUNCTION:
3507 : // Given an ascending monotonic table with n entries, find an index i
3508 : // such that ytab(i-1) < y <= ytab(i)
3509 :
3510 : // METHODOLOGY EMPLOYED:
3511 : // binary search
3512 :
3513 : // Return value
3514 : int SearchAscTable;
3515 :
3516 : int Ih; // Intex for upper end of interval
3517 : int Il; // Index for lower end of interval
3518 : int Im; // Index for midpoint of interval
3519 : Real64 Yh; // Table value for upper end of interval
3520 : Real64 Yl; // Table value for lower end of interval
3521 : Real64 Ym; // Table value for midpoint of interval
3522 :
3523 3350 : Yh = ytab(n);
3524 3350 : Yl = ytab(1);
3525 3350 : Ih = n;
3526 3350 : Il = 1;
3527 3350 : if (y < Yl) {
3528 0 : SearchAscTable = 1;
3529 0 : return SearchAscTable;
3530 3350 : } else if (y > Yh) {
3531 0 : SearchAscTable = n;
3532 0 : return SearchAscTable;
3533 : }
3534 : while (true) {
3535 14487 : if (Ih - Il <= 1) {
3536 3350 : break;
3537 : }
3538 11137 : Im = (Ih + Il) / 2;
3539 11137 : Ym = ytab(Im);
3540 11137 : if (y <= Ym) {
3541 3709 : Yh = Ym;
3542 3709 : Ih = Im;
3543 : } else {
3544 7428 : Yl = Ym;
3545 7428 : Il = Im;
3546 : }
3547 : }
3548 :
3549 3350 : SearchAscTable = Ih;
3550 :
3551 3350 : return SearchAscTable;
3552 : }
3553 :
3554 : } // namespace WindowComplexManager
3555 :
3556 : } // namespace EnergyPlus
|
