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