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