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