LCOV - code coverage report
Current view: top level - EnergyPlus - HeatBalanceIntRadExchange.cc (source / functions) Hit Total Coverage
Test: lcov.output.filtered Lines: 762 967 78.8 %
Date: 2023-01-17 19:17:23 Functions: 16 17 94.1 %

          Line data    Source code
       1             : // EnergyPlus, Copyright (c) 1996-2023, The Board of Trustees of the University of Illinois,
       2             : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
       3             : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
       4             : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
       5             : // contributors. All rights reserved.
       6             : //
       7             : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
       8             : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
       9             : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
      10             : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
      11             : // derivative works, and perform publicly and display publicly, and to permit others to do so.
      12             : //
      13             : // Redistribution and use in source and binary forms, with or without modification, are permitted
      14             : // provided that the following conditions are met:
      15             : //
      16             : // (1) Redistributions of source code must retain the above copyright notice, this list of
      17             : //     conditions and the following disclaimer.
      18             : //
      19             : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
      20             : //     conditions and the following disclaimer in the documentation and/or other materials
      21             : //     provided with the distribution.
      22             : //
      23             : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
      24             : //     the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
      25             : //     used to endorse or promote products derived from this software without specific prior
      26             : //     written permission.
      27             : //
      28             : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
      29             : //     without changes from the version obtained under this License, or (ii) Licensee makes a
      30             : //     reference solely to the software portion of its product, Licensee must refer to the
      31             : //     software as "EnergyPlus version X" software, where "X" is the version number Licensee
      32             : //     obtained under this License and may not use a different name for the software. Except as
      33             : //     specifically required in this Section (4), Licensee shall not use in a company name, a
      34             : //     product name, in advertising, publicity, or other promotional activities any name, trade
      35             : //     name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
      36             : //     similar designation, without the U.S. Department of Energy's prior written consent.
      37             : //
      38             : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
      39             : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
      40             : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
      41             : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
      42             : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
      43             : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
      44             : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
      45             : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
      46             : // POSSIBILITY OF SUCH DAMAGE.
      47             : 
      48             : // C++ Headers
      49             : #include <cassert>
      50             : #include <cmath>
      51             : 
      52             : // ObjexxFCL Headers
      53             : #include <ObjexxFCL/Array.functions.hh>
      54             : #include <ObjexxFCL/ArrayS.functions.hh>
      55             : #include <ObjexxFCL/Fmath.hh>
      56             : 
      57             : // EnergyPlus Headers
      58             : #include <EnergyPlus/Construction.hh>
      59             : #include <EnergyPlus/Data/EnergyPlusData.hh>
      60             : #include <EnergyPlus/DataHeatBalSurface.hh>
      61             : #include <EnergyPlus/DataHeatBalance.hh>
      62             : #include <EnergyPlus/DataIPShortCuts.hh>
      63             : #include <EnergyPlus/DataSurfaces.hh>
      64             : #include <EnergyPlus/DataSystemVariables.hh>
      65             : #include <EnergyPlus/DataViewFactorInformation.hh>
      66             : #include <EnergyPlus/DisplayRoutines.hh>
      67             : #include <EnergyPlus/General.hh>
      68             : #include <EnergyPlus/HeatBalanceIntRadExchange.hh>
      69             : #include <EnergyPlus/InputProcessing/InputProcessor.hh>
      70             : #include <EnergyPlus/Material.hh>
      71             : #include <EnergyPlus/UtilityRoutines.hh>
      72             : #include <EnergyPlus/WindowEquivalentLayer.hh>
      73             : 
      74             : namespace EnergyPlus {
      75             : 
      76             : namespace HeatBalanceIntRadExchange {
      77             :     // Module containing the routines dealing with the interior radiant exchange
      78             :     // between surfaces.
      79             : 
      80             :     // MODULE INFORMATION:
      81             :     //       AUTHOR         Rick Strand
      82             :     //       DATE WRITTEN   September 2000
      83             :     //       MODIFIED       Aug 2001, FW: recalculate ScriptF for a zone if window interior
      84             :     //                       shade/blind status is different from previous time step. This is
      85             :     //                       because ScriptF, which is used to calculate interior LW
      86             :     //                       exchange between surfaces, depends on inside surface emissivities,
      87             :     //                       which, for a window, depends on whether or not an interior
      88             :     //                       shade or blind is in place.
      89             :     //       RE-ENGINEERED  na
      90             : 
      91             :     // PURPOSE OF THIS MODULE:
      92             :     // Part of the heat balance modularization/re-engineering.  Purpose of this
      93             :     // module is to replace the MRT with RBAL method of modeling radiant exchange
      94             :     // between interior surfaces.
      95             : 
      96             :     // METHODOLOGY EMPLOYED:
      97             :     // Standard EnergyPlus methodology
      98             : 
      99             :     // REFERENCES:
     100             :     // ASHRAE Loads Toolkit "Script F" routines by Curt Pedersen
     101             :     // Hottel, H.C., and A.F. Sarofim. "Radiative Transfer" (mainly chapter 3),
     102             :     //  McGraw-Hill, Inc., New York, 1967.
     103             : 
     104             :     // OTHER NOTES: none
     105             : 
     106             :     // Using/Aliasing
     107             :     using namespace DataHeatBalance;
     108             :     using namespace DataSurfaces;
     109             :     using namespace DataSystemVariables;
     110             :     using namespace DataViewFactorInformation;
     111             : 
     112    19227093 :     void CalcInteriorRadExchange(EnergyPlusData &state,
     113             :                                  Array1S<Real64> const SurfaceTemp,   // Current surface temperatures
     114             :                                  int const SurfIterations,            // Number of iterations in calling subroutine
     115             :                                  Array1D<Real64> &NetLWRadToSurf,     // Net long wavelength radiant exchange from other surfaces
     116             :                                  Optional_int_const ZoneToResimulate, // if passed in, then only calculate for this zone
     117             : #ifdef EP_Count_Calls
     118             :                                  std::string_view const CalledFrom)
     119             : #else
     120             :                                  [[maybe_unused]] std::string_view const CalledFrom)
     121             : #endif
     122             :     {
     123             : 
     124             :         // SUBROUTINE INFORMATION:
     125             :         //       AUTHOR         Rick Strand
     126             :         //       DATE WRITTEN   September 2000
     127             :         //       MODIFIED       6/18/01, FCW: calculate IR on windows
     128             :         //                      Jan 2002, FCW: add blinds with movable slats
     129             :         //                      Sep 2011 LKL/BG - resimulate only zones needing it for Radiant systems
     130             : 
     131             :         // PURPOSE OF THIS SUBROUTINE:
     132             :         // Determines the interior radiant exchange between surfaces using
     133             :         // Hottel's ScriptF method for the grey interchange between surfaces
     134             :         // in an enclosure.
     135             : 
     136             :         // REFERENCES:
     137             :         // Hottel, H. C. and A. F. Sarofim, Radiative Transfer, Ch 3, McGraw Hill, 1967.
     138             : 
     139             :         // Types
     140             :         typedef Array1D<Real64>::size_type size_type;
     141             : 
     142             :         // Using/Aliasing
     143             :         using WindowEquivalentLayer::EQLWindowInsideEffectiveEmiss;
     144             : 
     145    19227093 :         Real64 constexpr StefanBoltzmannConst(5.6697e-8); // Stefan-Boltzmann constant in W/(m2*K4)
     146             : 
     147             :         bool IntShadeOrBlindStatusChanged; // True if status of interior shade or blind on at least
     148             :         // one window in a zone has changed from previous time step
     149             : 
     150    19227093 :         auto &SurfaceTempRad(state.dataHeatBalIntRadExchg->SurfaceTempRad);
     151    19227093 :         auto &SurfaceTempInKto4th(state.dataHeatBalIntRadExchg->SurfaceTempInKto4th);
     152    19227093 :         auto &SurfaceEmiss(state.dataHeatBalIntRadExchg->SurfaceEmiss);
     153             : 
     154    19227093 :         if (state.dataHeatBalIntRadExchg->CalcInteriorRadExchangefirstTime) {
     155         771 :             SurfaceTempRad.allocate(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs);
     156         771 :             SurfaceTempInKto4th.allocate(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs);
     157         771 :             SurfaceEmiss.allocate(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs);
     158         771 :             state.dataHeatBalIntRadExchg->CalcInteriorRadExchangefirstTime = false;
     159         771 :             if (state.dataSysVars->DeveloperFlag) {
     160           0 :                 DisplayString(state, " OMP turned off, HBIRE loop executed in serial");
     161             :             }
     162             :         }
     163             : 
     164    19227093 :         if (state.dataGlobal->KickOffSimulation || state.dataGlobal->KickOffSizing) return;
     165             : 
     166    19060283 :         bool const PartialResimulate(present(ZoneToResimulate));
     167             : 
     168             : #ifdef EP_Count_Calls
     169             :         if (!PartialResimulate) {
     170             :             ++state.dataTimingsData->NumIntRadExchange_Calls;
     171             :         } else {
     172             :             ++state.dataTimingsData->NumIntRadExchangeZ_Calls;
     173             :         }
     174             :         if (CalledFrom.empty()) {
     175             :             // do nothing
     176             :         } else if (CalledFrom == "Main") {
     177             :             ++state.dataTimingsData->NumIntRadExchangeMain_Calls;
     178             :         } else if (CalledFrom == "Outside") {
     179             :             ++state.dataTimingsData->NumIntRadExchangeOSurf_Calls;
     180             :         } else if (CalledFrom == "Inside") {
     181             :             ++state.dataTimingsData->NumIntRadExchangeISurf_Calls;
     182             :         }
     183             : #endif
     184             : 
     185    19060283 :         int startEnclosure = 1;
     186    19060283 :         int endEnclosure = state.dataViewFactor->NumOfRadiantEnclosures;
     187    19060283 :         if (PartialResimulate) {
     188             :             // ToDo: For now, use min and max enclosure numbers associated with this zone, this could include unrelated enclosures
     189     3226987 :             startEnclosure = state.dataHeatBal->Zone(ZoneToResimulate).zoneRadEnclosureFirst;
     190     3226987 :             endEnclosure = state.dataHeatBal->Zone(ZoneToResimulate).zoneRadEnclosureLast;
     191     6453974 :             for (int enclosureNum = startEnclosure; enclosureNum <= endEnclosure; ++enclosureNum) {
     192     3226987 :                 auto const &enclosure(state.dataViewFactor->EnclRadInfo(enclosureNum));
     193    26869398 :                 for (int i : enclosure.SurfacePtr) {
     194    23642411 :                     NetLWRadToSurf(i) = 0.0;
     195    23642411 :                     state.dataSurface->SurfWinIRfromParentZone(i) = 0.0;
     196             :                 }
     197             :             }
     198             :         } else {
     199    15833296 :             NetLWRadToSurf = 0.0;
     200  1117925858 :             for (int SurfNum = 1; SurfNum <= state.dataSurface->TotSurfaces; SurfNum++)
     201  1102092562 :                 state.dataSurface->SurfWinIRfromParentZone(SurfNum) = 0.0;
     202             :         }
     203             : 
     204   145871301 :         for (int enclosureNum = startEnclosure; enclosureNum <= endEnclosure; ++enclosureNum) {
     205             : 
     206   126811018 :             auto &zone_info(state.dataViewFactor->EnclRadInfo(enclosureNum));
     207   126811018 :             auto &zone_ScriptF(zone_info.ScriptF); // Tuned Transposed
     208   126811018 :             auto &zone_SurfacePtr(zone_info.SurfacePtr);
     209   126811018 :             int const n_zone_Surfaces(zone_info.NumOfSurfaces);
     210   126811018 :             size_type const s_zone_Surfaces(n_zone_Surfaces);
     211             : 
     212             :             // Calculate ScriptF if first time step in environment and surface heat-balance iterations not yet started;
     213             :             // recalculate ScriptF if status of window interior shades or blinds has changed from
     214             :             // previous time step. This recalculation is required since ScriptF depends on the inside
     215             :             // emissivity of the inside surfaces, which, for windows, is (1) the emissivity of the
     216             :             // inside face of the inside glass layer if there is no interior shade/blind, or (2) the effective
     217             :             // emissivity of the shade/blind if the shade/blind is in place. (The "effective emissivity"
     218             :             // in this case is (1) the shade/blind emissivity if the shade/blind IR transmittance is zero,
     219             :             // or (2) a weighted average of the shade/blind emissivity and inside glass emissivity if the
     220             :             // shade/blind IR transmittance is not zero (which is sometimes the case for a "shade" and
     221             :             // usually the case for a blind). It assumed for switchable glazing that the inside surface
     222             :             // emissivity does not change if the glazing is switched on or off.
     223             : 
     224             :             // Determine if status of interior shade/blind on one or more windows in the zone has changed
     225             :             // from previous time step.  Also make a check for any changes in interior movable insulation.
     226             : 
     227   126811018 :             if (SurfIterations == 0) {
     228             : 
     229             :                 bool IntMovInsulChanged; // True if the status of interior movable insulation has changed
     230             : 
     231    57532998 :                 IntShadeOrBlindStatusChanged = false;
     232    57532998 :                 IntMovInsulChanged = false;
     233             : 
     234    57532998 :                 if (!state.dataGlobal->BeginEnvrnFlag) { // Check for change in shade/blind status
     235   551923422 :                     for (int const SurfNum : zone_SurfacePtr) {
     236   494452195 :                         if (IntShadeOrBlindStatusChanged || IntMovInsulChanged)
     237             :                             break; // Need only check if one window's status or one movable insulation status has changed
     238   494451046 :                         if (state.dataConstruction->Construct(state.dataSurface->Surface(SurfNum).Construction).TypeIsWindow) {
     239    67058296 :                             WinShadingType ShadeFlag = state.dataSurface->SurfWinShadingFlag(SurfNum);
     240    67058296 :                             WinShadingType ShadeFlagPrev = state.dataSurface->SurfWinExtIntShadePrevTS(SurfNum);
     241    67058296 :                             if (ShadeFlagPrev != ShadeFlag && (ANY_INTERIOR_SHADE_BLIND(ShadeFlagPrev) || ANY_INTERIOR_SHADE_BLIND(ShadeFlag)))
     242       25581 :                                 IntShadeOrBlindStatusChanged = true;
     243    67082452 :                             if (state.dataSurface->SurfWinWindowModelType(SurfNum) == WindowModel::EQL &&
     244             :                                 state.dataWindowEquivLayer
     245       24156 :                                     ->CFS(state.dataConstruction->Construct(state.dataSurface->Surface(SurfNum).Construction).EQLConsPtr)
     246       24156 :                                     .ISControlled) {
     247           0 :                                 IntShadeOrBlindStatusChanged = true;
     248             :                             }
     249             :                         } else {
     250   427392750 :                             if (state.dataSurface->AnyMovableInsulation) UpdateMovableInsulationFlag(state, IntMovInsulChanged, SurfNum);
     251             :                         }
     252             :                     }
     253             :                 }
     254             : 
     255   115040415 :                 if (IntShadeOrBlindStatusChanged || IntMovInsulChanged ||
     256    57507417 :                     state.dataGlobal->BeginEnvrnFlag) { // Calc inside surface emissivities for this time step
     257      795676 :                     for (int ZoneSurfNum = 1; ZoneSurfNum <= n_zone_Surfaces; ++ZoneSurfNum) {
     258      709473 :                         int const SurfNum = zone_SurfacePtr(ZoneSurfNum);
     259      709473 :                         int const ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
     260      709473 :                         zone_info.Emissivity(ZoneSurfNum) = state.dataHeatBalSurf->SurfAbsThermalInt(SurfNum);
     261      811866 :                         if (state.dataConstruction->Construct(ConstrNum).TypeIsWindow &&
     262      102393 :                             ANY_INTERIOR_SHADE_BLIND(state.dataSurface->SurfWinShadingFlag(SurfNum))) {
     263       14367 :                             zone_info.Emissivity(ZoneSurfNum) = state.dataHeatBalSurf->SurfAbsThermalInt(SurfNum);
     264             :                         }
     265      709509 :                         if (state.dataSurface->SurfWinWindowModelType(SurfNum) == WindowModel::EQL &&
     266          36 :                             state.dataWindowEquivLayer->CFS(state.dataConstruction->Construct(ConstrNum).EQLConsPtr).ISControlled) {
     267           0 :                             zone_info.Emissivity(ZoneSurfNum) = EQLWindowInsideEffectiveEmiss(state, ConstrNum);
     268             :                         }
     269             :                     }
     270             : 
     271       86203 :                     if (state.dataHeatBalIntRadExchg->CarrollMethod) {
     272         456 :                         CalcFp(n_zone_Surfaces, zone_info.Emissivity, zone_info.FMRT, zone_info.Fp);
     273             :                     } else {
     274       85747 :                         CalcScriptF(state, n_zone_Surfaces, zone_info.Area, zone_info.F, zone_info.Emissivity, zone_ScriptF);
     275             :                         // precalc - multiply by StefanBoltzmannConstant
     276       85747 :                         zone_ScriptF *= StefanBoltzmannConst;
     277             :                     }
     278             :                 }
     279             : 
     280             :             } // End of check if SurfIterations = 0
     281             : 
     282             :             // Set surface emissivities and temperatures
     283             :             // Also, for Carroll method, calculate numerators and denominators of radiant temperature
     284   126811018 :             Real64 CarrollMRTNumerator(0.0);
     285   126811018 :             Real64 CarrollMRTDenominator(0.0);
     286             :             Real64 CarrollMRTInKTo4th; // Carroll MRT
     287  1222703011 :             for (size_type ZoneSurfNum = 0; ZoneSurfNum < s_zone_Surfaces; ++ZoneSurfNum) {
     288  1095891993 :                 int const SurfNum = zone_SurfacePtr[ZoneSurfNum];
     289  1095891993 :                 auto const &surface_window(state.dataSurface->SurfaceWindow(SurfNum));
     290  1095891993 :                 int const ConstrNum = state.dataSurface->Surface(SurfNum).Construction;
     291  1095891993 :                 auto const &construct(state.dataConstruction->Construct(ConstrNum));
     292  1095891993 :                 if (construct.WindowTypeEQL) {
     293       37770 :                     SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
     294       37770 :                     SurfaceEmiss[ZoneSurfNum] = EQLWindowInsideEffectiveEmiss(state, ConstrNum);
     295  1095854223 :                 } else if (construct.WindowTypeBSDF && state.dataSurface->SurfWinShadingFlag(SurfNum) == WinShadingType::IntShade) {
     296      100650 :                     SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
     297      100650 :                     SurfaceEmiss[ZoneSurfNum] = surface_window.EffShBlindEmiss[0] + surface_window.EffGlassEmiss[0];
     298  1095753573 :                 } else if (construct.WindowTypeBSDF) {
     299      146661 :                     SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
     300      146661 :                     SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
     301  1095606912 :                 } else if (construct.TypeIsWindow && state.dataSurface->SurfWinOriginalClass(SurfNum) != SurfaceClass::TDD_Diffuser) {
     302   146751590 :                     if (SurfIterations == 0 && NOT_SHADED(state.dataSurface->SurfWinShadingFlag(SurfNum))) {
     303             :                         // If the window is bare this TS and it is the first time through we use the previous TS glass
     304             :                         // temperature whether or not the window was shaded in the previous TS. If the window was shaded
     305             :                         // the previous time step this temperature is a better starting value than the shade temperature.
     306    66543514 :                         SurfaceTempRad[ZoneSurfNum] = surface_window.ThetaFace(2 * construct.TotGlassLayers) - DataGlobalConstants::KelvinConv;
     307    66543514 :                         SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
     308             :                         // For windows with an interior shade or blind an effective inside surface temp
     309             :                         // and emiss is used here that is a weighted combination of shade/blind and glass temp and emiss.
     310    80208076 :                     } else if (ANY_INTERIOR_SHADE_BLIND(state.dataSurface->SurfWinShadingFlag(SurfNum))) {
     311      255807 :                         SurfaceTempRad[ZoneSurfNum] = state.dataSurface->SurfWinEffInsSurfTemp(SurfNum);
     312      255807 :                         SurfaceEmiss[ZoneSurfNum] = state.dataHeatBalSurf->SurfAbsThermalInt(SurfNum);
     313             :                     } else {
     314    79952269 :                         SurfaceTempRad[ZoneSurfNum] = SurfaceTemp(SurfNum);
     315    79952269 :                         SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
     316             :                     }
     317             :                 } else {
     318   948855322 :                     SurfaceTempRad[ZoneSurfNum] = SurfaceTemp(SurfNum);
     319   948855322 :                     SurfaceEmiss[ZoneSurfNum] = construct.InsideAbsorpThermal;
     320             :                 }
     321  1095891993 :                 if (state.dataHeatBalIntRadExchg->CarrollMethod) {
     322     9384094 :                     CarrollMRTNumerator += SurfaceTempRad[ZoneSurfNum] * zone_info.Fp[ZoneSurfNum] * zone_info.Area[ZoneSurfNum];
     323     9384094 :                     CarrollMRTDenominator += zone_info.Fp[ZoneSurfNum] * zone_info.Area[ZoneSurfNum];
     324             :                 }
     325  1095891993 :                 SurfaceTempInKto4th[ZoneSurfNum] = pow_4(SurfaceTempRad[ZoneSurfNum] + DataGlobalConstants::KelvinConv);
     326             :             }
     327             : 
     328   126811018 :             if (state.dataHeatBalIntRadExchg->CarrollMethod) {
     329     1128467 :                 if (CarrollMRTDenominator > 0.0) {
     330     1128467 :                     CarrollMRTInKTo4th = pow_4(CarrollMRTNumerator / CarrollMRTDenominator + DataGlobalConstants::KelvinConv);
     331             :                 } else {
     332             :                     // Likely only one surface in this enclosure
     333           0 :                     CarrollMRTInKTo4th = 293.15; // arbitrary value, IR will be zero
     334             :                 }
     335             :             }
     336             : 
     337             :             // These are the money loops
     338   126811018 :             if (state.dataHeatBalIntRadExchg->CarrollMethod) {
     339    10512561 :                 for (size_type RecZoneSurfNum = 0; RecZoneSurfNum < s_zone_Surfaces; ++RecZoneSurfNum) {
     340     9384094 :                     int const RecSurfNum = zone_SurfacePtr[RecZoneSurfNum];
     341     9384094 :                     int const ConstrNumRec = state.dataSurface->Surface(RecSurfNum).Construction;
     342     9384094 :                     auto const &rec_construct(state.dataConstruction->Construct(ConstrNumRec));
     343     9384094 :                     auto &netLWRadToRecSurf(NetLWRadToSurf(RecSurfNum));
     344     9384094 :                     if (rec_construct.TypeIsWindow) {
     345             :                         //                        auto& rec_surface_window(SurfaceWindow(RecSurfNum));
     346      712716 :                         Real64 CarrollMRTInKTo4thWin = CarrollMRTInKTo4th; // arbitrary value, IR will be zero
     347      712716 :                         Real64 CarrollMRTNumeratorWin(0.0);
     348      712716 :                         Real64 CarrollMRTDenominatorWin(0.0);
     349     6414444 :                         for (size_type SendZoneSurfNum = 0; SendZoneSurfNum < s_zone_Surfaces; ++SendZoneSurfNum) {
     350     5701728 :                             if (SendZoneSurfNum != RecZoneSurfNum) {
     351     4989012 :                                 CarrollMRTNumeratorWin +=
     352     4989012 :                                     SurfaceTempRad[SendZoneSurfNum] * zone_info.Fp[SendZoneSurfNum] * zone_info.Area[SendZoneSurfNum];
     353     4989012 :                                 CarrollMRTDenominatorWin += zone_info.Fp[SendZoneSurfNum] * zone_info.Area[SendZoneSurfNum];
     354             :                             }
     355             :                         }
     356      712716 :                         if (CarrollMRTDenominatorWin > 0.0) {
     357      712716 :                             CarrollMRTInKTo4thWin = pow_4(CarrollMRTNumeratorWin / CarrollMRTDenominatorWin + DataGlobalConstants::KelvinConv);
     358             :                         }
     359      712716 :                         state.dataSurface->SurfWinIRfromParentZone(RecSurfNum) +=
     360      712716 :                             (zone_info.Fp[RecZoneSurfNum] * CarrollMRTInKTo4thWin) / SurfaceEmiss[RecZoneSurfNum];
     361             :                     }
     362     9384094 :                     netLWRadToRecSurf += zone_info.Fp[RecZoneSurfNum] * (CarrollMRTInKTo4th - SurfaceTempInKto4th[RecZoneSurfNum]);
     363             :                 }
     364             :             } else {
     365  1212190450 :                 for (size_type RecZoneSurfNum = 0; RecZoneSurfNum < s_zone_Surfaces; ++RecZoneSurfNum) {
     366  1086507899 :                     int const RecSurfNum = zone_SurfacePtr[RecZoneSurfNum];
     367  1086507899 :                     int const ConstrNumRec = state.dataSurface->Surface(RecSurfNum).Construction;
     368  1086507899 :                     auto const &rec_construct(state.dataConstruction->Construct(ConstrNumRec));
     369  1086507899 :                     auto &netLWRadToRecSurf(NetLWRadToSurf(RecSurfNum));
     370             : 
     371             :                     // Calculate net long-wave radiation for opaque surfaces and incident
     372             :                     // long-wave radiation for windows.
     373  1086507899 :                     if (rec_construct.TypeIsWindow) {      // Window
     374             :                                                            //                        auto& rec_surface_window(SurfaceWindow(RecSurfNum));
     375   146352701 :                         Real64 scriptF_acc(0.0);           // Local accumulator
     376   146352701 :                         Real64 netLWRadToRecSurf_cor(0.0); // Correction
     377   146352701 :                         Real64 IRfromParentZone_acc(0.0);  // Local accumulator
     378  2103033848 :                         for (size_type SendZoneSurfNum = 0; SendZoneSurfNum < s_zone_Surfaces; ++SendZoneSurfNum) {
     379  1956681147 :                             size_type lSR = RecZoneSurfNum * s_zone_Surfaces + SendZoneSurfNum;
     380  1956681147 :                             Real64 const scriptF(zone_ScriptF[lSR]); // [ lSR ] == ( SendZoneSurfNum+1, RecZoneSurfNum+1 )
     381  1956681147 :                             Real64 const scriptF_temp_ink_4th(scriptF * SurfaceTempInKto4th[SendZoneSurfNum]);
     382             :                             // Calculate interior LW incident on window rather than net LW for use in window layer heat balance calculation.
     383  1956681147 :                             IRfromParentZone_acc += scriptF_temp_ink_4th;
     384             : 
     385  1956681147 :                             if (RecZoneSurfNum != SendZoneSurfNum) {
     386  1810328446 :                                 scriptF_acc += scriptF;
     387             :                             } else {
     388   146352701 :                                 netLWRadToRecSurf_cor = scriptF_temp_ink_4th;
     389             :                             }
     390             : 
     391             :                             // Per BG -- this should never happened.  (CR6346,CR6550 caused this to be put in.  Now removed. LKL 1/2013)
     392             :                             //          IF (SurfaceWindow(RecSurfNum)%IRfromParentZone < 0.0) THEN
     393             :                             //            CALL ShowRecurringWarningErrorAtEnd(state, 'CalcInteriorRadExchange: Window_IRFromParentZone negative,
     394             :                             //            Window="'// &
     395             :                             //                TRIM(Surface(RecSurfNum)%Name)//'"',  &
     396             :                             //                SurfaceWindow(RecSurfNum)%IRErrCount)
     397             :                             //            CALL ShowRecurringContinueErrorAtEnd(state, '..occurs in Zone="'//TRIM(Surface(RecSurfNum)%ZoneName)//  &
     398             :                             //                '", reset to 0.0 for remaining calculations.',SurfaceWindow(RecSurfNum)%IRErrCountC)
     399             :                             //            SurfaceWindow(RecSurfNum)%IRfromParentZone=0.0
     400             :                             //          ENDIF
     401             :                         }
     402   146352701 :                         netLWRadToRecSurf += IRfromParentZone_acc - netLWRadToRecSurf_cor - (scriptF_acc * SurfaceTempInKto4th[RecZoneSurfNum]);
     403   146352701 :                         state.dataSurface->SurfWinIRfromParentZone(RecSurfNum) += IRfromParentZone_acc / SurfaceEmiss[RecZoneSurfNum];
     404             :                     } else {
     405   940155198 :                         Real64 netLWRadToRecSurf_acc(0.0); // Local accumulator
     406   940155198 :                         zone_ScriptF[RecZoneSurfNum * s_zone_Surfaces + RecZoneSurfNum] = 0;
     407  9923194780 :                         for (size_type SendZoneSurfNum = 0; SendZoneSurfNum < s_zone_Surfaces; ++SendZoneSurfNum) {
     408  8983039582 :                             size_type lSR = RecZoneSurfNum * s_zone_Surfaces + SendZoneSurfNum;
     409  8983039582 :                             netLWRadToRecSurf_acc +=
     410 17966079164 :                                 zone_ScriptF[lSR] * (SurfaceTempInKto4th[SendZoneSurfNum] -
     411  8983039582 :                                                      SurfaceTempInKto4th[RecZoneSurfNum]); // [ lSR ] == ( SendZoneSurfNum+1, RecZoneSurfNum+1 )
     412             :                         }
     413   940155198 :                         netLWRadToRecSurf += netLWRadToRecSurf_acc;
     414             :                     }
     415             :                 }
     416             :             }
     417             :         }
     418             : 
     419             :         // Automatic Surface Multipliers: Update values of surfaces not simulated
     420    19060283 :         if (state.dataSurface->UseRepresentativeSurfaceCalculations) {
     421     2007960 :             for (auto SurfNum : state.dataSurface->AllHTSurfaceList) {
     422     2004498 :                 auto &RepSurfNum = state.dataSurface->Surface(SurfNum).RepresentativeCalcSurfNum;
     423     2004498 :                 if (SurfNum != RepSurfNum) {
     424      453522 :                     state.dataSurface->SurfWinIRfromParentZone(SurfNum) = state.dataSurface->SurfWinIRfromParentZone(RepSurfNum);
     425      453522 :                     NetLWRadToSurf(SurfNum) = NetLWRadToSurf(RepSurfNum);
     426             :                 }
     427             :             }
     428             :         }
     429             :     }
     430             : 
     431      181260 :     void UpdateMovableInsulationFlag(EnergyPlusData &state, bool &MovableInsulationChange, int const SurfNum)
     432             :     {
     433             : 
     434             :         // SUBROUTINE INFORMATION:
     435             :         //       AUTHOR         Rick Strand
     436             :         //       DATE WRITTEN   July 2016
     437             : 
     438             :         // PURPOSE OF THIS SUBROUTINE:
     439             :         // To determine if any changes in interior movable insulation have happened.
     440             :         // If there have been changes due to a schedule change AND a change in properties,
     441             :         // then the matrices which are used to calculate interior radiation must be recalculated.
     442             : 
     443      181260 :         MovableInsulationChange = false;
     444             : 
     445      181260 :         if (state.dataHeatBalSurf->SurfMovInsulIntPresent(SurfNum) != state.dataHeatBalSurf->SurfMovInsulIntPresentPrevTS(SurfNum)) {
     446           0 :             auto const &thissurf(state.dataSurface->Surface(SurfNum));
     447           0 :             Real64 AbsorpDiff = std::abs(state.dataConstruction->Construct(thissurf.Construction).InsideAbsorpThermal -
     448           0 :                                          state.dataMaterial->Material(state.dataSurface->SurfMaterialMovInsulInt(SurfNum)).AbsorpThermal);
     449           0 :             if (AbsorpDiff > 0.01) {
     450           0 :                 MovableInsulationChange = true;
     451             :             }
     452             :         }
     453      181260 :     }
     454             : 
     455         771 :     void InitInteriorRadExchange(EnergyPlusData &state)
     456             :     {
     457             : 
     458             :         // SUBROUTINE INFORMATION:
     459             :         //       AUTHOR         Rick Strand
     460             :         //       DATE WRITTEN   September 2000
     461             :         //       MODIFIED       na
     462             :         //       RE-ENGINEERED  na
     463             : 
     464             :         // PURPOSE OF THIS SUBROUTINE:
     465             :         // Initializes the various parameters for Hottel's ScriptF method for
     466             :         // the grey interchange between surfaces in an enclosure.
     467             : 
     468             :         // Using/Aliasing
     469             : 
     470             :         using General::ScanForReports;
     471             : 
     472             :         // SUBROUTINE PARAMETER DEFINITIONS:
     473             : 
     474             :         // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
     475             :         bool NoUserInputF; // Logical flag signifying no input F's for zone
     476         771 :         bool ErrorsFound(false);
     477             :         Real64 CheckValue1;
     478             :         Real64 CheckValue2;
     479             :         Real64 FinalCheckValue;
     480        1542 :         Array2D<Real64> SaveApproximateViewFactors; // Save for View Factor reporting
     481             :         Real64 RowSum;
     482             :         Real64 FixedRowSum;
     483             :         int NumIterations;
     484        1542 :         std::string Option1; // view factor report option
     485             : 
     486         771 :         auto &ViewFactorReport(state.dataHeatBalIntRadExchg->ViewFactorReport);
     487             : 
     488         771 :         ScanForReports(state, "ViewFactorInfo", ViewFactorReport, _, Option1);
     489             : 
     490         771 :         if (ViewFactorReport) { // Print heading
     491           6 :             print(state.files.eio, "{}\n", "! <Surface View Factor and Grey Interchange Information>");
     492           6 :             print(state.files.eio, "{}\n", "! <View Factor - Zone/Enclosure Information>,Zone/Enclosure Name,Number of Surfaces");
     493           6 :             print(state.files.eio,
     494             :                   "{}\n",
     495           6 :                   "! <View Factor - Surface Information>,Surface Name,Surface Class,Area {m2},Azimuth,Tilt,Thermal Emissivity,#Sides,Vertices");
     496           6 :             print(state.files.eio, "{}\n", "! <View Factor / Grey Interchange Type>,Surface Name(s)");
     497           6 :             print(state.files.eio, "{}\n", "! <View Factor>,Surface Name,Surface Class,Row Sum,View Factors for each Surface");
     498             :         }
     499             : 
     500         771 :         state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs = 0;
     501        5582 :         for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfRadiantEnclosures; ++enclosureNum) {
     502        4811 :             auto &thisEnclosure(state.dataViewFactor->EnclRadInfo(enclosureNum));
     503        4811 :             if (enclosureNum == 1) {
     504         741 :                 if (state.dataGlobal->DisplayAdvancedReportVariables) {
     505          20 :                     print(state.files.eio,
     506             :                           "{}\n",
     507             :                           "! <Surface View Factor Check Values>,Zone/Enclosure Name,Original Check Value,Calculated Fixed Check Value,Final Check "
     508          20 :                           "Value,Number of Iterations,Fixed RowSum Convergence,Used RowSum Convergence");
     509             :                 }
     510             :             }
     511        4811 :             int numEnclosureSurfaces = 0;
     512        9629 :             for (int spaceNum : thisEnclosure.spaceNums) {
     513             :                 // Note that Space.Surfaces only includes HT surfs, see SurfaceGeometry::CreateMissingSpaces
     514             :                 // But it also includes air boundary surfaces which need to be excluded here
     515       47038 :                 for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
     516       42220 :                     if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
     517             :                     // Automatic Surface Multipliers: Only include representative surfaces
     518       42206 :                     if (surfNum == state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
     519       42075 :                         ++numEnclosureSurfaces;
     520             :                     }
     521             :                 }
     522             :             }
     523        4811 :             thisEnclosure.NumOfSurfaces = numEnclosureSurfaces;
     524        4811 :             state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs =
     525        4811 :                 max(state.dataHeatBalIntRadExchg->MaxNumOfRadEnclosureSurfs, numEnclosureSurfaces);
     526        4811 :             if (numEnclosureSurfaces < 1) ShowFatalError(state, "No surfaces in an enclosure in InitInteriorRadExchange");
     527             : 
     528             :             // Allocate the parts of the derived type
     529        4811 :             thisEnclosure.F.dimension(numEnclosureSurfaces, numEnclosureSurfaces, 0.0);
     530        4811 :             thisEnclosure.ScriptF.dimension(numEnclosureSurfaces, numEnclosureSurfaces, 0.0);
     531        4811 :             thisEnclosure.Area.dimension(numEnclosureSurfaces, 0.0);
     532        4811 :             thisEnclosure.Emissivity.dimension(numEnclosureSurfaces, 0.0);
     533        4811 :             thisEnclosure.Azimuth.dimension(numEnclosureSurfaces, 0.0);
     534        4811 :             thisEnclosure.Tilt.dimension(numEnclosureSurfaces, 0.0);
     535        4811 :             thisEnclosure.Fp.dimension(numEnclosureSurfaces, 1.0);
     536        4811 :             thisEnclosure.FMRT.dimension(numEnclosureSurfaces, 0.0);
     537        4811 :             thisEnclosure.SurfacePtr.dimension(numEnclosureSurfaces, 0);
     538             : 
     539             :             // Initialize the enclosure surface arrays
     540        4811 :             int enclosureSurfNum = 0;
     541        9629 :             for (int const spaceNum : thisEnclosure.spaceNums) {
     542        4818 :                 int priorZoneTotEnclSurfs = enclosureSurfNum;
     543       47038 :                 for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
     544       42220 :                     if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
     545             :                     // Automatic Surface Multipliers: Only include representative surfaces
     546       42206 :                     if (surfNum == state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
     547       42075 :                         ++enclosureSurfNum;
     548       42075 :                         thisEnclosure.SurfacePtr(enclosureSurfNum) = surfNum;
     549       42075 :                         thisEnclosure.Area(enclosureSurfNum) = state.dataSurface->Surface(surfNum).Area;
     550       42075 :                         thisEnclosure.Emissivity(enclosureSurfNum) =
     551       42075 :                             state.dataConstruction->Construct(state.dataSurface->Surface(surfNum).Construction).InsideAbsorpThermal;
     552       42075 :                         thisEnclosure.Azimuth(enclosureSurfNum) = state.dataSurface->Surface(surfNum).Azimuth;
     553       42075 :                         thisEnclosure.Tilt(enclosureSurfNum) = state.dataSurface->Surface(surfNum).Tilt;
     554             :                     }
     555             :                 }
     556             : 
     557       47038 :                 for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
     558       42220 :                     if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
     559             :                     // Map non-representative surfaces
     560       42206 :                     if (surfNum != state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
     561             :                         // Automatic Surface Multipliers: search for corresponding enclosure surface number
     562        1642 :                         for (int enclSNum = priorZoneTotEnclSurfs + 1; enclSNum <= enclosureSurfNum; ++enclSNum) {
     563        1511 :                             if (thisEnclosure.SurfacePtr(enclSNum) == state.dataSurface->Surface(surfNum).RepresentativeCalcSurfNum) {
     564             :                                 // enclosureSurfMap[allSurfNum] = enclSNum;
     565             :                                 //  Increase the area for the combined enclosure surface
     566         131 :                                 thisEnclosure.Area(enclSNum) += state.dataSurface->Surface(surfNum).Area;
     567             :                             }
     568             :                         }
     569             :                     }
     570             :                     // Store SurfaceReportNums to maintain original reporting order
     571      239207 :                     for (int enclSNum = priorZoneTotEnclSurfs + 1; enclSNum <= enclosureSurfNum; ++enclSNum) {
     572      239076 :                         if (thisEnclosure.SurfacePtr(enclSNum) == state.dataSurface->AllSurfaceListReportOrder[surfNum - 1]) {
     573       42075 :                             thisEnclosure.SurfaceReportNums.push_back(enclSNum);
     574       42075 :                             break;
     575             :                         }
     576             :                     }
     577             :                 }
     578             :             }
     579             : 
     580        4813 :             if (thisEnclosure.NumOfSurfaces == 1) {
     581             :                 // If there is only one surface in a zone, then there is no radiant exchange
     582           2 :                 thisEnclosure.F = 0.0;
     583           2 :                 thisEnclosure.ScriptF = 0.0;
     584           2 :                 thisEnclosure.Fp = 0.0;
     585           2 :                 thisEnclosure.FMRT = 0.0;
     586           2 :                 if (state.dataGlobal->DisplayAdvancedReportVariables)
     587           1 :                     print(state.files.eio, "Surface View Factor Check Values,{},0,0,0,-1,0,0\n", thisEnclosure.Name);
     588           2 :                 continue; // Go to the next enclosure in the loop
     589             :             }
     590             : 
     591             :             // Process View Factors
     592        4809 :             if (state.dataHeatBalIntRadExchg->CarrollMethod) {
     593             : 
     594             :                 // User View Factors cannot be used with Carroll method.
     595          19 :                 if (state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "ZoneProperty:UserViewFactors:BySurfaceName")) {
     596           0 :                     ShowWarningError(state, "ZoneProperty:UserViewFactors:BySurfaceName objects have been defined, however View");
     597           0 :                     ShowContinueError(state, "  Factors are not used when Zone Radiant Exchange Algorithm is set to CarrollMRT.");
     598             :                 }
     599          19 :                 CalcFMRT(state, thisEnclosure.NumOfSurfaces, thisEnclosure.Area, thisEnclosure.FMRT);
     600          19 :                 CalcFp(thisEnclosure.NumOfSurfaces, thisEnclosure.Emissivity, thisEnclosure.FMRT, thisEnclosure.Fp);
     601             :             } else {
     602             :                 //  Get user supplied view factors if available in idf.
     603             : 
     604        4790 :                 NoUserInputF = true;
     605             : 
     606        9580 :                 std::string cCurrentModuleObject = "ZoneProperty:UserViewFactors:BySurfaceName";
     607        4790 :                 int NumZonesWithUserFbyS = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, cCurrentModuleObject);
     608        4790 :                 if (NumZonesWithUserFbyS > 0) {
     609             : 
     610           2 :                     GetInputViewFactorsbyName(state,
     611             :                                               thisEnclosure.Name,
     612             :                                               thisEnclosure.NumOfSurfaces,
     613             :                                               thisEnclosure.F,
     614             :                                               thisEnclosure.SurfacePtr,
     615             :                                               NoUserInputF,
     616             :                                               ErrorsFound); // Obtains user input view factors from input file
     617             :                 }
     618             : 
     619        4790 :                 if (NoUserInputF) {
     620             : 
     621             :                     // Calculate the view factors and make sure they satisfy reciprocity
     622        4789 :                     CalcApproximateViewFactors(state,
     623             :                                                thisEnclosure.NumOfSurfaces,
     624             :                                                thisEnclosure.Area,
     625             :                                                thisEnclosure.Azimuth,
     626             :                                                thisEnclosure.Tilt,
     627             :                                                thisEnclosure.F,
     628             :                                                thisEnclosure.SurfacePtr);
     629             :                 }
     630             : 
     631        4790 :                 if (ViewFactorReport) { // Allocate and save user or approximate view factors for reporting.
     632          28 :                     SaveApproximateViewFactors.allocate(thisEnclosure.NumOfSurfaces, thisEnclosure.NumOfSurfaces);
     633          28 :                     SaveApproximateViewFactors = thisEnclosure.F;
     634             :                 }
     635             : 
     636        4790 :                 bool anyIntMassInZone = DoesZoneHaveInternalMass(state, thisEnclosure.NumOfSurfaces, thisEnclosure.SurfacePtr);
     637        4790 :                 FixViewFactors(state,
     638             :                                thisEnclosure.NumOfSurfaces,
     639             :                                thisEnclosure.Area,
     640             :                                thisEnclosure.F,
     641             :                                thisEnclosure.Name,
     642             :                                thisEnclosure.spaceNums,
     643             :                                CheckValue1,
     644             :                                CheckValue2,
     645             :                                FinalCheckValue,
     646             :                                NumIterations,
     647             :                                FixedRowSum,
     648             :                                anyIntMassInZone);
     649             : 
     650             :                 // Calculate the script F factors
     651        4790 :                 CalcScriptF(state, thisEnclosure.NumOfSurfaces, thisEnclosure.Area, thisEnclosure.F, thisEnclosure.Emissivity, thisEnclosure.ScriptF);
     652        4790 :                 if (ViewFactorReport) { // Write to SurfInfo File
     653             :                     // Zone Surface Information Output
     654          56 :                     print(
     655          28 :                         state.files.eio, "Surface View Factor - Zone/Enclosure Information,{},{}\n", thisEnclosure.Name, thisEnclosure.NumOfSurfaces);
     656             : 
     657         266 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     658         714 :                         print(state.files.eio,
     659             :                               "Surface View Factor - Surface Information,{},{},{:.4R},{:.4R},{:.4R},{:.4R},{}",
     660         238 :                               state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
     661         476 :                               cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Class),
     662             :                               thisEnclosure.Area(SurfNum),
     663             :                               thisEnclosure.Azimuth(SurfNum),
     664             :                               thisEnclosure.Tilt(SurfNum),
     665             :                               thisEnclosure.Emissivity(SurfNum),
     666         476 :                               state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides);
     667        1190 :                         for (int Vindex = 1; Vindex <= state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides; ++Vindex) {
     668         952 :                             auto &Vertex = state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Vertex(Vindex);
     669         952 :                             print(state.files.eio, ",{:.4R},{:.4R},{:.4R}", Vertex.x, Vertex.y, Vertex.z);
     670             :                         }
     671         238 :                         print(state.files.eio, "\n");
     672             :                     }
     673             : 
     674          28 :                     print(state.files.eio, "Approximate or User Input ViewFactors,To Surface,Surface Class,RowSum");
     675         266 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     676         238 :                         print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
     677             :                     }
     678          28 :                     print(state.files.eio, "\n");
     679             : 
     680         266 :                     for (int Findex : thisEnclosure.SurfaceReportNums) {
     681         238 :                         RowSum = sum(SaveApproximateViewFactors(_, Findex));
     682         476 :                         print(state.files.eio,
     683             :                               "{},{},{},{:.4R}",
     684             :                               "View Factor",
     685         238 :                               state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
     686         476 :                               cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
     687         238 :                               RowSum);
     688        2400 :                         for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     689        2162 :                             print(state.files.eio, ",{:.4R}", SaveApproximateViewFactors(SurfNum, Findex));
     690             :                         }
     691         238 :                         print(state.files.eio, "\n");
     692             :                     }
     693             : 
     694          28 :                     print(state.files.eio, "Final ViewFactors,To Surface,Surface Class,RowSum");
     695         266 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     696         238 :                         print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
     697             :                     }
     698          28 :                     print(state.files.eio, "\n");
     699             : 
     700         266 :                     for (int Findex : thisEnclosure.SurfaceReportNums) {
     701         238 :                         RowSum = sum(thisEnclosure.F(_, Findex));
     702         476 :                         print(state.files.eio,
     703             :                               "{},{},{},{:.4R}",
     704             :                               "View Factor",
     705         238 :                               state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
     706         476 :                               cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
     707         238 :                               RowSum);
     708        2400 :                         for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     709        2162 :                             print(state.files.eio, ",{:.4R}", thisEnclosure.F(SurfNum, Findex));
     710             :                         }
     711         238 :                         print(state.files.eio, "\n");
     712             :                     }
     713             : 
     714          28 :                     if (Option1 == "IDF") {
     715             :                         // TODO Both "original" and "final" print the same output. This is likely a bug
     716             :                         // (discovered while updating output to {fmt}
     717             :                         // see:
     718             :                         // https://github.com/NREL/EnergyPlusArchive/commit/1c08247853c297dce59f3f53cde47ccfa67720c0#diff-124964a7e9b73ce494c1952ab1acdeeb
     719           8 :                         print(state.files.debug, "{}\n", "!======== original input factors ===========================");
     720           8 :                         print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
     721          88 :                         for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     722         956 :                             for (int Findex : thisEnclosure.SurfaceReportNums) {
     723        3504 :                                 print(state.files.debug,
     724             :                                       "  {},{},{:.6R}",
     725         876 :                                       state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
     726         876 :                                       state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
     727         876 :                                       thisEnclosure.F(Findex, SurfNum));
     728         876 :                                 if (!(SurfNum == thisEnclosure.NumOfSurfaces && Findex == thisEnclosure.NumOfSurfaces)) {
     729         868 :                                     print(state.files.debug, ",\n");
     730             :                                 } else {
     731           8 :                                     print(state.files.debug, ";\n");
     732             :                                 }
     733             :                             }
     734             :                         }
     735           8 :                         print(state.files.debug, "{}\n", "!============= end of data ======================");
     736             : 
     737           8 :                         print(state.files.debug, "{}\n", "!============ final view factors =======================");
     738           8 :                         print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
     739          88 :                         for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     740         956 :                             for (int Findex : thisEnclosure.SurfaceReportNums) {
     741        3504 :                                 print(state.files.debug,
     742             :                                       "  {},{},{:.6R}",
     743         876 :                                       state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
     744         876 :                                       state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
     745         876 :                                       thisEnclosure.F(Findex, SurfNum));
     746         876 :                                 if (!(SurfNum == thisEnclosure.SurfaceReportNums.back() && Findex == thisEnclosure.SurfaceReportNums.back())) {
     747         868 :                                     print(state.files.debug, ",\n");
     748             :                                 } else {
     749           8 :                                     print(state.files.debug, ";\n");
     750             :                                 }
     751             :                             }
     752             :                         }
     753           8 :                         print(state.files.debug, "{}\n", "!============= end of data ======================");
     754             :                     }
     755             : 
     756          28 :                     print(state.files.eio, "Script F Factors,X Surface");
     757         266 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     758         238 :                         print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
     759             :                     }
     760          28 :                     print(state.files.eio, "\n");
     761         266 :                     for (int Findex : thisEnclosure.SurfaceReportNums) {
     762         238 :                         print(state.files.eio, "{},{}", "Script F Factor", state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name);
     763        2400 :                         for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     764        2162 :                             print(state.files.eio, ",{:.4R}", thisEnclosure.ScriptF(Findex, SurfNum));
     765             :                         }
     766         238 :                         print(state.files.eio, "\n");
     767             :                     }
     768             : 
     769             :                     // Deallocate saved approximate/user view factors
     770          28 :                     SaveApproximateViewFactors.deallocate();
     771             :                 }
     772             : 
     773        4790 :                 RowSum = 0.0;
     774       46705 :                 for (int Findex : thisEnclosure.SurfaceReportNums) {
     775       41915 :                     RowSum += sum(thisEnclosure.F(_, Findex));
     776             :                 }
     777        4790 :                 RowSum = std::abs(RowSum - thisEnclosure.NumOfSurfaces);
     778        4790 :                 FixedRowSum = std::abs(FixedRowSum - thisEnclosure.NumOfSurfaces);
     779        4790 :                 if (state.dataGlobal->DisplayAdvancedReportVariables) {
     780         112 :                     print(state.files.eio,
     781             :                           "Surface View Factor Check Values,{},{:.6R},{:.6R},{:.6R},{},{:.6R},{:.6R}\n",
     782             :                           thisEnclosure.Name,
     783             :                           CheckValue1,
     784             :                           CheckValue2,
     785             :                           FinalCheckValue,
     786             :                           NumIterations,
     787             :                           FixedRowSum,
     788          56 :                           RowSum);
     789             :                 }
     790             :             }
     791             :         }
     792             : 
     793         771 :         if (ErrorsFound) {
     794           0 :             ShowFatalError(state, "InitInteriorRadExchange: Errors found during initialization of radiant exchange.  Program terminated.");
     795             :         }
     796         771 :     }
     797             : 
     798         771 :     void InitSolarViewFactors(EnergyPlusData &state)
     799             :     {
     800             : 
     801             :         // Initializes view factors for diffuse solar distribution between surfaces in an enclosure.
     802             : 
     803        1542 :         Array2D<Real64> SaveApproximateViewFactors; // Save for View Factor reporting
     804        1542 :         std::string Option1;                        // view factor report option
     805             : 
     806         771 :         bool ErrorsFound = false;
     807         771 :         bool ViewFactorReport = false;
     808         771 :         General::ScanForReports(state, "ViewFactorInfo", ViewFactorReport, _, Option1);
     809             : 
     810         771 :         if (ViewFactorReport) { // Print heading
     811           6 :             print(state.files.eio, "{}\n", "! <Solar View Factor Information>");
     812           6 :             print(state.files.eio, "{}\n", "! <Solar View Factor - Zone/Enclosure Information>,Zone/Enclosure Name,Number of Surfaces");
     813           6 :             print(state.files.eio,
     814             :                   "{}\n",
     815           6 :                   "! <Solar View Factor - Surface Information>,Surface Name,Surface Class,Area {m2},Azimuth,Tilt,Solar Absorbtance,#Sides,Vertices");
     816           6 :             print(state.files.eio, "{}\n", "! <Solar View Factor / Interchange Type>,Surface Name(s)");
     817           6 :             print(state.files.eio, "{}\n", "! <Solar View Factor>,Surface Name,Surface Class,Row Sum,View Factors for each Surface");
     818             :         }
     819             : 
     820        1542 :         std::string cCurrentModuleObject = "ZoneProperty:UserViewFactors:BySurfaceName";
     821         771 :         int NumZonesWithUserFbyS = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, cCurrentModuleObject);
     822         771 :         if (NumZonesWithUserFbyS > 0) AlignInputViewFactors(state, cCurrentModuleObject, ErrorsFound);
     823             : 
     824        5582 :         for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfSolarEnclosures; ++enclosureNum) {
     825        4811 :             auto &thisEnclosure(state.dataViewFactor->EnclSolInfo(enclosureNum));
     826        4811 :             if (enclosureNum == 1) {
     827         741 :                 if (state.dataGlobal->DisplayAdvancedReportVariables)
     828          20 :                     print(state.files.eio,
     829             :                           "{}\n",
     830             :                           "! <Solar View Factor Check Values>,Zone/Enclosure Name,Original Check Value,Calculated Fixed Check "
     831             :                           "Value,Final Check Value,Number of Iterations,Fixed RowSum Convergence,Used RowSum "
     832          20 :                           "Convergence");
     833             :             }
     834        4811 :             int numEnclosureSurfaces = 0;
     835        9629 :             for (int spaceNum : thisEnclosure.spaceNums) {
     836             :                 // Note that Space.Surfaces only includes HT surfs, see SurfaceGeometry::CreateMissingSpaces
     837             :                 // But it also includes air boundary surfaces which need to be excluded here
     838       47038 :                 for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
     839       42220 :                     if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
     840       42206 :                     ++numEnclosureSurfaces;
     841             :                 }
     842             :             }
     843        4811 :             thisEnclosure.NumOfSurfaces = numEnclosureSurfaces;
     844        4811 :             if (numEnclosureSurfaces < 1) ShowFatalError(state, "No surfaces in an enclosure in InitSolarViewFactors");
     845             : 
     846             :             // Allocate the parts of the derived type
     847        4811 :             thisEnclosure.F.dimension(numEnclosureSurfaces, numEnclosureSurfaces, 0.0);
     848        4811 :             thisEnclosure.Area.dimension(numEnclosureSurfaces, 0.0);
     849        4811 :             thisEnclosure.SolAbsorptance.dimension(numEnclosureSurfaces, 0.0);
     850        4811 :             thisEnclosure.Azimuth.dimension(numEnclosureSurfaces, 0.0);
     851        4811 :             thisEnclosure.Tilt.dimension(numEnclosureSurfaces, 0.0);
     852        4811 :             thisEnclosure.SurfacePtr.dimension(numEnclosureSurfaces, 0);
     853             : 
     854             :             // Initialize the surface pointer array
     855        4811 :             int enclosureSurfNum = 0;
     856        9629 :             for (int const spaceNum : thisEnclosure.spaceNums) {
     857        4818 :                 int priorZoneTotEnclSurfs = enclosureSurfNum;
     858       47038 :                 for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
     859       42220 :                     if (state.dataSurface->Surface(surfNum).IsAirBoundarySurf) continue;
     860       42206 :                     ++enclosureSurfNum;
     861       42206 :                     thisEnclosure.SurfacePtr(enclosureSurfNum) = surfNum;
     862             :                     // Store pointers back to here
     863       42206 :                     state.dataSurface->Surface(surfNum).SolarEnclSurfIndex = enclosureSurfNum;
     864       42206 :                     state.dataSurface->Surface(surfNum).SolarEnclIndex = enclosureNum;
     865       42206 :                     state.dataSurface->Surface(surfNum).RadEnclIndex = enclosureNum; // Radiant and Solar enclosures are parallel for now
     866       90388 :                     if ((state.dataConstruction->Construct(state.dataSurface->Surface(surfNum).Construction).TypeIsWindow) &&
     867       42220 :                         (state.dataSurface->Surface(surfNum).ExtBoundCond > 0) &&
     868          14 :                         (state.dataSurface->Surface(surfNum).BaseSurf != surfNum)) { // Interzone window present
     869          14 :                         thisEnclosure.HasInterZoneWindow = true;
     870             :                     }
     871             :                 }
     872             :                 // Store SurfaceReportNums to maintain original reporting order
     873       47038 :                 for (int surfNum : state.dataHeatBal->space(spaceNum).surfaces) {
     874      241498 :                     for (int enclSNum = priorZoneTotEnclSurfs + 1; enclSNum <= enclosureSurfNum; ++enclSNum) {
     875      241484 :                         if (thisEnclosure.SurfacePtr(enclSNum) == state.dataSurface->AllSurfaceListReportOrder[surfNum - 1]) {
     876       42206 :                             thisEnclosure.SurfaceReportNums.push_back(enclSNum);
     877       42206 :                             break;
     878             :                         }
     879             :                     }
     880             :                 }
     881             :             }
     882             :             // Initialize the area and related arrays
     883       47017 :             for (int enclSurfNum = 1; enclSurfNum <= thisEnclosure.NumOfSurfaces; ++enclSurfNum) {
     884       42206 :                 int const SurfNum = thisEnclosure.SurfacePtr(enclSurfNum);
     885       42206 :                 thisEnclosure.Area(enclSurfNum) = state.dataSurface->Surface(SurfNum).Area;
     886       42206 :                 thisEnclosure.SolAbsorptance(enclSurfNum) =
     887       42206 :                     state.dataConstruction->Construct(state.dataSurface->Surface(SurfNum).Construction).InsideAbsorpSolar;
     888       42206 :                 thisEnclosure.Azimuth(enclSurfNum) = state.dataSurface->Surface(SurfNum).Azimuth;
     889       42206 :                 thisEnclosure.Tilt(enclSurfNum) = state.dataSurface->Surface(SurfNum).Tilt;
     890             :             }
     891             : 
     892        4813 :             if (thisEnclosure.NumOfSurfaces == 1) {
     893             :                 // If there is only one surface in a zone, then there is no solar distribution
     894           2 :                 if (state.dataGlobal->DisplayAdvancedReportVariables)
     895           1 :                     print(state.files.eio, "Solar View Factor Check Values,{},0,0,0,-1,0,0\n", thisEnclosure.Name);
     896           2 :                 continue; // Go to the next enclosure in the loop
     897             :             }
     898             : 
     899             :             //  Get user supplied view factors if available in idf.
     900             : 
     901        4809 :             bool NoUserInputF = true;
     902             : 
     903        4809 :             if (NumZonesWithUserFbyS > 0) {
     904             : 
     905           2 :                 GetInputViewFactorsbyName(state,
     906             :                                           thisEnclosure.Name,
     907             :                                           thisEnclosure.NumOfSurfaces,
     908             :                                           thisEnclosure.F,
     909             :                                           thisEnclosure.SurfacePtr,
     910             :                                           NoUserInputF,
     911             :                                           ErrorsFound); // Obtains user input view factors from input file
     912             :             }
     913             : 
     914        4809 :             if (NoUserInputF) {
     915             : 
     916             :                 // Calculate the view factors and make sure they satisfy reciprocity
     917        4808 :                 CalcApproximateViewFactors(state,
     918             :                                            thisEnclosure.NumOfSurfaces,
     919             :                                            thisEnclosure.Area,
     920             :                                            thisEnclosure.Azimuth,
     921             :                                            thisEnclosure.Tilt,
     922             :                                            thisEnclosure.F,
     923             :                                            thisEnclosure.SurfacePtr);
     924             :             }
     925             : 
     926        4809 :             if (ViewFactorReport) { // Allocate and save user or approximate view factors for reporting.
     927          28 :                 SaveApproximateViewFactors.allocate(thisEnclosure.NumOfSurfaces, thisEnclosure.NumOfSurfaces);
     928          28 :                 SaveApproximateViewFactors = thisEnclosure.F;
     929             :             }
     930             : 
     931        4809 :             Real64 CheckValue1 = 0.0;
     932        4809 :             Real64 CheckValue2 = 0.0;
     933        4809 :             Real64 FinalCheckValue = 0.0;
     934        4809 :             Real64 FixedRowSum = 0.0;
     935        4809 :             int NumIterations = 0;
     936             : 
     937        4809 :             bool anyIntMassInZone = DoesZoneHaveInternalMass(state, thisEnclosure.NumOfSurfaces, thisEnclosure.SurfacePtr);
     938        4809 :             FixViewFactors(state,
     939             :                            thisEnclosure.NumOfSurfaces,
     940             :                            thisEnclosure.Area,
     941             :                            thisEnclosure.F,
     942             :                            thisEnclosure.Name,
     943             :                            thisEnclosure.spaceNums,
     944             :                            CheckValue1,
     945             :                            CheckValue2,
     946             :                            FinalCheckValue,
     947             :                            NumIterations,
     948             :                            FixedRowSum,
     949             :                            anyIntMassInZone);
     950             : 
     951        4809 :             if (ViewFactorReport) { // Write to SurfInfo File
     952             :                 // Zone Surface Information Output
     953          28 :                 print(state.files.eio, "Solar View Factor - Zone/Enclosure Information,{},{}\n", thisEnclosure.Name, thisEnclosure.NumOfSurfaces);
     954             : 
     955         266 :                 for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     956         714 :                     print(state.files.eio,
     957             :                           "Solar View Factor - Surface Information,{},{},{:.4R},{:.4R},{:.4R},{:.4R},{}",
     958         238 :                           state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
     959         476 :                           cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Class),
     960             :                           thisEnclosure.Area(SurfNum),
     961             :                           thisEnclosure.Azimuth(SurfNum),
     962             :                           thisEnclosure.Tilt(SurfNum),
     963             :                           thisEnclosure.SolAbsorptance(SurfNum),
     964         476 :                           state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides);
     965             : 
     966        1190 :                     for (int Vindex = 1; Vindex <= state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Sides; ++Vindex) {
     967         952 :                         auto &Vertex = state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Vertex(Vindex);
     968         952 :                         print(state.files.eio, ",{:.4R},{:.4R},{:.4R}", Vertex.x, Vertex.y, Vertex.z);
     969             :                     }
     970         238 :                     print(state.files.eio, "\n");
     971             :                 }
     972             : 
     973          28 :                 print(state.files.eio, "Approximate or User Input Solar ViewFactors,To Surface,Surface Class,RowSum");
     974         266 :                 for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     975         238 :                     print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
     976             :                 }
     977          28 :                 print(state.files.eio, "\n");
     978             : 
     979         266 :                 for (int Findex : thisEnclosure.SurfaceReportNums) {
     980         238 :                     Real64 RowSum = sum(SaveApproximateViewFactors(_, Findex));
     981         476 :                     print(state.files.eio,
     982             :                           "Solar View Factor,{},{},{:.4R}",
     983         238 :                           state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
     984         476 :                           cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
     985         238 :                           RowSum);
     986        2400 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     987        2162 :                         print(state.files.eio, ",{:.4R}", SaveApproximateViewFactors(SurfNum, Findex));
     988             :                     }
     989         238 :                     print(state.files.eio, "\n");
     990             :                 }
     991             :             }
     992             : 
     993        4809 :             if (ViewFactorReport) {
     994          28 :                 print(state.files.eio, "Final Solar ViewFactors,To Surface,Surface Class,RowSum");
     995         266 :                 for (int SurfNum : thisEnclosure.SurfaceReportNums) {
     996         238 :                     print(state.files.eio, ",{}", state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name);
     997             :                 }
     998          28 :                 print(state.files.eio, "\n");
     999             : 
    1000         266 :                 for (int Findex : thisEnclosure.SurfaceReportNums) {
    1001         238 :                     Real64 RowSum = sum(thisEnclosure.F(_, Findex));
    1002         476 :                     print(state.files.eio,
    1003             :                           "{},{},{},{:.4R}",
    1004             :                           "Solar View Factor",
    1005         238 :                           state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
    1006         476 :                           cSurfaceClass(state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Class),
    1007         238 :                           RowSum);
    1008        2400 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
    1009        2162 :                         print(state.files.eio, ",{:.4R}", thisEnclosure.F(SurfNum, Findex));
    1010             :                     }
    1011         238 :                     print(state.files.eio, "\n");
    1012             :                 }
    1013             : 
    1014          28 :                 if (Option1 == "IDF") {
    1015             :                     // TODO Both "original" and "final" print the same output. This is likely a bug
    1016             :                     // see:
    1017             :                     // https://github.com/NREL/EnergyPlusArchive/commit/1c08247853c297dce59f3f53cde47ccfa67720c0#diff-124964a7e9b73ce494c1952ab1acdeeb
    1018           8 :                     print(state.files.debug, "{}\n", "!======== original input factors ===========================");
    1019           8 :                     print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
    1020          88 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
    1021         956 :                         for (int Findex : thisEnclosure.SurfaceReportNums) {
    1022        3504 :                             print(state.files.debug,
    1023             :                                   "  {},{},{:.6R}",
    1024         876 :                                   state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
    1025         876 :                                   state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
    1026         876 :                                   thisEnclosure.F(Findex, SurfNum));
    1027         876 :                             if (!(SurfNum == thisEnclosure.NumOfSurfaces && Findex == thisEnclosure.NumOfSurfaces)) {
    1028         868 :                                 print(state.files.debug, ",\n");
    1029             :                             } else {
    1030           8 :                                 print(state.files.debug, ";\n");
    1031             :                             }
    1032             :                         }
    1033             :                     }
    1034           8 :                     print(state.files.debug, "{}\n", "!============= end of data ======================");
    1035             : 
    1036           8 :                     print(state.files.debug, "{}\n", "!============ final view factors =======================");
    1037           8 :                     print(state.files.debug, "ZoneProperty:UserViewFactors:BySurfaceName,{},\n", thisEnclosure.Name);
    1038          88 :                     for (int SurfNum : thisEnclosure.SurfaceReportNums) {
    1039         956 :                         for (int Findex : thisEnclosure.SurfaceReportNums) {
    1040        3504 :                             print(state.files.debug,
    1041             :                                   "  {},{},{:.6R}",
    1042         876 :                                   state.dataSurface->Surface(thisEnclosure.SurfacePtr(SurfNum)).Name,
    1043         876 :                                   state.dataSurface->Surface(thisEnclosure.SurfacePtr(Findex)).Name,
    1044         876 :                                   thisEnclosure.F(Findex, SurfNum));
    1045         876 :                             if (!(SurfNum == thisEnclosure.NumOfSurfaces && Findex == thisEnclosure.NumOfSurfaces)) {
    1046         868 :                                 print(state.files.debug, ",\n");
    1047             :                             } else {
    1048           8 :                                 print(state.files.debug, ";\n");
    1049             :                             }
    1050             :                         }
    1051             :                     }
    1052           8 :                     print(state.files.debug, "{}\n", "!============= end of data ======================");
    1053             :                 }
    1054             :             }
    1055             : 
    1056        4809 :             if (ViewFactorReport) { // Deallocate saved approximate/user view factors
    1057          28 :                 SaveApproximateViewFactors.deallocate();
    1058             :             }
    1059             : 
    1060        4809 :             Real64 RowSum = 0.0;
    1061       47013 :             for (int Findex : thisEnclosure.SurfaceReportNums) {
    1062       42204 :                 RowSum += sum(thisEnclosure.F(_, Findex));
    1063             :             }
    1064        4809 :             RowSum = std::abs(RowSum - thisEnclosure.NumOfSurfaces);
    1065        4809 :             FixedRowSum = std::abs(FixedRowSum - thisEnclosure.NumOfSurfaces);
    1066        4809 :             if (state.dataGlobal->DisplayAdvancedReportVariables) {
    1067         112 :                 print(state.files.eio,
    1068             :                       "Solar View Factor Check Values,{},{:.6R},{:.6R},{:.6R},{},{:.6R},{:.6R}\n",
    1069             :                       thisEnclosure.Name,
    1070             :                       CheckValue1,
    1071             :                       CheckValue2,
    1072             :                       FinalCheckValue,
    1073             :                       NumIterations,
    1074             :                       FixedRowSum,
    1075          56 :                       RowSum);
    1076             :             }
    1077             :         }
    1078             : 
    1079         771 :         if (ErrorsFound) {
    1080           0 :             ShowFatalError(state, "InitSolarViewFactors: Errors found during initialization of diffuse solar distribution.  Program terminated.");
    1081             :         }
    1082         771 :     }
    1083             : 
    1084           0 :     void GetInputViewFactors(EnergyPlusData &state,
    1085             :                              std::string const &ZoneName,              // Needed to check for user input view factors.
    1086             :                              int const N,                              // NUMBER OF SURFACES
    1087             :                              Array2A<Real64> F,                        // USER INPUT DIRECT VIEW FACTOR MATRIX (N X N)
    1088             :                              [[maybe_unused]] const Array1D_int &SPtr, // pointer to actual surface number
    1089             :                              bool &NoUserInputF,                       // Flag signifying no input F's for this
    1090             :                              bool &ErrorsFound                         // True when errors are found in number of fields vs max args
    1091             :     )
    1092             :     {
    1093             : 
    1094             :         // SUBROUTINE INFORMATION:
    1095             :         //       AUTHOR         Curt Pedersen
    1096             :         //       DATE WRITTEN   September 2005
    1097             :         //       MODIFIED       Linda Lawrie;September 2010
    1098             :         //       RE-ENGINEERED  na
    1099             : 
    1100             :         // PURPOSE OF THIS SUBROUTINE:
    1101             :         // This routine gets the user view factor info.
    1102             : 
    1103             :         // Using/Aliasing
    1104             : 
    1105             :         // Argument array dimensioning
    1106           0 :         F.dim(N, N);
    1107             :         // EP_SIZE_CHECK(SPtr, N);
    1108             : 
    1109             :         // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    1110             :         //  INTEGER   :: NumZonesWithUserF
    1111             :         int UserFZoneIndex;
    1112             :         int NumAlphas;
    1113             :         int NumNums;
    1114             :         int IOStat;
    1115             :         int index;
    1116             :         int inx1;
    1117             :         int inx2;
    1118             : 
    1119           0 :         NoUserInputF = true;
    1120           0 :         UserFZoneIndex = state.dataInputProcessing->inputProcessor->getObjectItemNum(state, "ZoneProperty:UserViewFactors", ZoneName);
    1121           0 :         auto &cCurrentModuleObject = state.dataIPShortCut->cCurrentModuleObject;
    1122           0 :         if (UserFZoneIndex > 0) {
    1123           0 :             NoUserInputF = false;
    1124             : 
    1125           0 :             state.dataInputProcessing->inputProcessor->getObjectItem(state,
    1126             :                                                                      "ZoneProperty:UserViewFactors",
    1127             :                                                                      UserFZoneIndex,
    1128           0 :                                                                      state.dataIPShortCut->cAlphaArgs,
    1129             :                                                                      NumAlphas,
    1130           0 :                                                                      state.dataIPShortCut->rNumericArgs,
    1131             :                                                                      NumNums,
    1132             :                                                                      IOStat,
    1133           0 :                                                                      state.dataIPShortCut->lNumericFieldBlanks,
    1134           0 :                                                                      state.dataIPShortCut->lAlphaFieldBlanks,
    1135           0 :                                                                      state.dataIPShortCut->cAlphaFieldNames,
    1136           0 :                                                                      state.dataIPShortCut->cNumericFieldNames);
    1137             : 
    1138           0 :             if (NumNums < 3 * pow_2(N)) {
    1139           0 :                 ShowSevereError(state, "GetInputViewFactors: " + cCurrentModuleObject + "=\"" + ZoneName + "\", not enough values.");
    1140           0 :                 ShowContinueError(state, format("...Number of input values [{}] is less than the required number=[{}].", NumNums, 3 * pow_2(N)));
    1141           0 :                 ErrorsFound = true;
    1142           0 :                 NumNums = 0;
    1143             :             }
    1144           0 :             F = 0.0;
    1145           0 :             for (index = 1; index <= NumNums; index += 3) {
    1146           0 :                 inx1 = state.dataIPShortCut->rNumericArgs(index);
    1147           0 :                 inx2 = state.dataIPShortCut->rNumericArgs(index + 1);
    1148           0 :                 F(inx2, inx1) = state.dataIPShortCut->rNumericArgs(index + 2);
    1149             :             }
    1150             :         }
    1151           0 :     }
    1152             : 
    1153           1 :     void AlignInputViewFactors(EnergyPlusData &state,
    1154             :                                std::string const &cCurrentModuleObject, // Object type
    1155             :                                bool &ErrorsFound                        // True when errors are found
    1156             :     )
    1157             :     {
    1158           2 :         auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(cCurrentModuleObject);
    1159           1 :         auto &instancesValue = instances.value();
    1160           2 :         for (auto instance = instancesValue.begin(); instance != instancesValue.end(); ++instance) {
    1161           1 :             auto const &fields = instance.value();
    1162             :             std::string const thisSpaceOrSpaceListName =
    1163           1 :                 UtilityRoutines::MakeUPPERCase(fields.at("zone_or_zonelist_or_space_or_spacelist_name").get<std::string>());
    1164             :             // do not mark object as used here - let GetInputViewFactorsbyName do that
    1165             : 
    1166             :             // Look for matching radiant enclosure name (solar enclosures have the same names)
    1167           1 :             bool enclMatchFound = false;
    1168           1 :             for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfRadiantEnclosures; ++enclosureNum) {
    1169           1 :                 auto &thisEnclosure(state.dataViewFactor->EnclRadInfo(enclosureNum));
    1170           1 :                 if (UtilityRoutines::SameString(thisSpaceOrSpaceListName, thisEnclosure.Name)) {
    1171             :                     // View factor space name matches enclosure name
    1172           1 :                     enclMatchFound = true;
    1173           1 :                     break;
    1174             :                 }
    1175             :             }
    1176           1 :             if (enclMatchFound) continue; // We're done with this instance
    1177             :             // Find matching SpaceList name
    1178           0 :             int spaceListNum = UtilityRoutines::FindItemInList(thisSpaceOrSpaceListName, state.dataHeatBal->spaceList);
    1179           0 :             int zoneListNum = 0;
    1180           0 :             int inputZoneNum = 0;
    1181           0 :             size_t listSize = 0;
    1182           0 :             if (spaceListNum > 0) {
    1183           0 :                 listSize = int(state.dataHeatBal->spaceList(spaceListNum).spaces.size());
    1184             :             } else {
    1185             :                 // Check Zone and ZoneList for backwards compatibility
    1186           0 :                 zoneListNum = UtilityRoutines::FindItemInList(thisSpaceOrSpaceListName, state.dataHeatBal->ZoneList);
    1187           0 :                 if (zoneListNum > 0) {
    1188           0 :                     for (int const zoneNum : state.dataHeatBal->ZoneList(zoneListNum).Zone) {
    1189           0 :                         listSize += state.dataHeatBal->Zone(zoneNum).spaceIndexes.size();
    1190             :                     }
    1191             :                 } else {
    1192           0 :                     inputZoneNum = UtilityRoutines::FindItemInList(thisSpaceOrSpaceListName, state.dataHeatBal->Zone);
    1193           0 :                     if (inputZoneNum > 0) {
    1194           0 :                         listSize = state.dataHeatBal->Zone(inputZoneNum).spaceIndexes.size();
    1195             :                     }
    1196             :                 }
    1197             :             }
    1198             : 
    1199           0 :             if (listSize > 0) {
    1200             :                 // Look for radiant enclosure with same list of spaces (solar enclosures are the same, so this is sufficient)
    1201           0 :                 for (int enclosureNum = 1; enclosureNum <= state.dataViewFactor->NumOfRadiantEnclosures; ++enclosureNum) {
    1202           0 :                     auto &thisEnclosure(state.dataViewFactor->EnclRadInfo(enclosureNum));
    1203           0 :                     bool anySpaceNotFound = false;
    1204             :                     // If the number of enclosure spaces is not the same as the number of spacelist space, go to the next enclosure
    1205           0 :                     if (listSize != thisEnclosure.spaceNums.size()) continue;
    1206           0 :                     if (spaceListNum > 0) {
    1207           0 :                         for (int sListSpaceNum : state.dataHeatBal->spaceList(spaceListNum).spaces) {
    1208             :                             // Search for matching spaces
    1209           0 :                             bool thisSpaceFound = false;
    1210           0 :                             for (int enclSpaceNum : thisEnclosure.spaceNums) {
    1211           0 :                                 if (enclSpaceNum == sListSpaceNum) {
    1212           0 :                                     thisSpaceFound = true;
    1213           0 :                                     break;
    1214             :                                 }
    1215             :                             }
    1216           0 :                             if (!thisSpaceFound) {
    1217           0 :                                 anySpaceNotFound = true;
    1218           0 :                                 break;
    1219             :                             }
    1220             :                         }
    1221           0 :                     } else if (zoneListNum > 0) {
    1222           0 :                         for (int zListZoneNum : state.dataHeatBal->ZoneList(zoneListNum).Zone) {
    1223           0 :                             for (int spaceNum : state.dataHeatBal->Zone(zListZoneNum).spaceIndexes) {
    1224             :                                 // Search for matching spaces
    1225           0 :                                 bool thisSpaceFound = false;
    1226           0 :                                 for (int enclSpaceNum : thisEnclosure.spaceNums) {
    1227           0 :                                     if (enclSpaceNum == spaceNum) {
    1228           0 :                                         thisSpaceFound = true;
    1229           0 :                                         break;
    1230             :                                     }
    1231             :                                 }
    1232           0 :                                 if (!thisSpaceFound) {
    1233           0 :                                     anySpaceNotFound = true;
    1234           0 :                                     break;
    1235             :                                 }
    1236             :                             }
    1237           0 :                             if (anySpaceNotFound) {
    1238           0 :                                 break;
    1239             :                             }
    1240             :                         }
    1241           0 :                     } else if (inputZoneNum > 0) {
    1242           0 :                         for (int spaceNum : state.dataHeatBal->Zone(inputZoneNum).spaceIndexes) {
    1243             :                             // Search for matching spaces
    1244           0 :                             bool thisSpaceFound = false;
    1245           0 :                             for (int enclSpaceNum : thisEnclosure.spaceNums) {
    1246           0 :                                 if (enclSpaceNum == spaceNum) {
    1247           0 :                                     thisSpaceFound = true;
    1248           0 :                                     break;
    1249             :                                 }
    1250             :                             }
    1251           0 :                             if (!thisSpaceFound) {
    1252           0 :                                 anySpaceNotFound = true;
    1253           0 :                                 break;
    1254             :                             }
    1255             :                         }
    1256             :                     }
    1257           0 :                     if (anySpaceNotFound) {
    1258           0 :                         continue; // On to the next enclosure
    1259             :                     } else {
    1260           0 :                         enclMatchFound = true;
    1261             :                         // If matching SpaceList or ZoneList or Zone found, set the radiant and solar enclosure names to match
    1262           0 :                         thisEnclosure.Name = thisSpaceOrSpaceListName;
    1263           0 :                         state.dataViewFactor->EnclSolInfo(enclosureNum).Name = thisSpaceOrSpaceListName;
    1264           0 :                         break; // We're done with radiant enclosures
    1265             :                     }
    1266             :                 }
    1267             :             }
    1268           0 :             if (!enclMatchFound) {
    1269           0 :                 if (spaceListNum > 0) {
    1270           0 :                     ShowSevereError(
    1271             :                         state,
    1272           0 :                         "AlignInputViewFactors: " + cCurrentModuleObject + "=\"" + thisSpaceOrSpaceListName +
    1273             :                             "\" found a matching SpaceList, but did not find a matching radiant or solar enclosure with the same spaces.");
    1274           0 :                     ErrorsFound = true;
    1275             : 
    1276           0 :                 } else if (zoneListNum > 0) {
    1277           0 :                     ShowSevereError(state,
    1278           0 :                                     "AlignInputViewFactors: " + cCurrentModuleObject + "=\"" + thisSpaceOrSpaceListName +
    1279             :                                         "\" found a matching ZoneList, but did not find a matching radiant or solar enclosure with the same spaces.");
    1280           0 :                     ErrorsFound = true;
    1281             : 
    1282             :                 } else {
    1283           0 :                     ShowSevereError(state,
    1284           0 :                                     "AlignInputViewFactors: " + cCurrentModuleObject + "=\"" + thisSpaceOrSpaceListName +
    1285             :                                         "\" did not find a matching radiant or solar enclosure name.");
    1286           0 :                     ErrorsFound = true;
    1287             :                 }
    1288             :             }
    1289             :         }
    1290           1 :     }
    1291             : 
    1292           4 :     void GetInputViewFactorsbyName(EnergyPlusData &state,
    1293             :                                    std::string const &EnclosureName, // Needed to check for user input view factors.
    1294             :                                    int const N,                      // NUMBER OF SURFACES
    1295             :                                    Array2A<Real64> F,                // USER INPUT DIRECT VIEW FACTOR MATRIX (N X N)
    1296             :                                    const Array1D_int &SPtr,          // pointer to actual surface number
    1297             :                                    bool &NoUserInputF,               // Flag signifying no input F's for this
    1298             :                                    bool &ErrorsFound                 // True when errors are found in number of fields vs max args
    1299             :     )
    1300             :     {
    1301             : 
    1302             :         // SUBROUTINE INFORMATION:
    1303             :         //       AUTHOR         Curt Pedersen
    1304             :         //       DATE WRITTEN   September 2005
    1305             :         //       MODIFIED       Linda Lawrie;September 2010
    1306             : 
    1307             :         // PURPOSE OF THIS SUBROUTINE:
    1308             :         // This routine gets the user view factor info for an enclosure which could be a zone or a group of zones
    1309             : 
    1310             :         // Using/Aliasing
    1311             : 
    1312             :         // Argument array dimensioning
    1313           4 :         F.dim(N, N);
    1314           4 :         EP_SIZE_CHECK(SPtr, N);
    1315             : 
    1316             :         // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    1317             :         int UserFZoneIndex;
    1318             :         int NumAlphas;
    1319             :         int NumNums;
    1320             :         int IOStat;
    1321             :         int index;
    1322             :         int numinx1;
    1323             :         int inx1;
    1324             :         int inx2;
    1325           8 :         Array1D_string enclosureSurfaceNames;
    1326           4 :         auto &cCurrentModuleObject = state.dataIPShortCut->cCurrentModuleObject;
    1327           4 :         NoUserInputF = true;
    1328           4 :         UserFZoneIndex = state.dataInputProcessing->inputProcessor->getObjectItemNum(
    1329             :             state, "ZoneProperty:UserViewFactors:BySurfaceName", "zone_or_zonelist_or_space_or_spacelist_name", EnclosureName);
    1330             : 
    1331           4 :         if (UserFZoneIndex > 0) {
    1332           2 :             enclosureSurfaceNames.allocate(N);
    1333          22 :             for (index = 1; index <= N; ++index) {
    1334          20 :                 enclosureSurfaceNames(index) = state.dataSurface->Surface(SPtr(index)).Name;
    1335             :             }
    1336           2 :             NoUserInputF = false;
    1337             : 
    1338          16 :             state.dataInputProcessing->inputProcessor->getObjectItem(state,
    1339             :                                                                      "ZoneProperty:UserViewFactors:BySurfaceName",
    1340             :                                                                      UserFZoneIndex,
    1341           2 :                                                                      state.dataIPShortCut->cAlphaArgs,
    1342             :                                                                      NumAlphas,
    1343           2 :                                                                      state.dataIPShortCut->rNumericArgs,
    1344             :                                                                      NumNums,
    1345             :                                                                      IOStat,
    1346           2 :                                                                      state.dataIPShortCut->lNumericFieldBlanks,
    1347           2 :                                                                      state.dataIPShortCut->lAlphaFieldBlanks,
    1348           2 :                                                                      state.dataIPShortCut->cAlphaFieldNames,
    1349           4 :                                                                      state.dataIPShortCut->cNumericFieldNames);
    1350             : 
    1351           2 :             F = 0.0;
    1352           2 :             numinx1 = 0;
    1353           2 :             if (NumNums < pow_2(N)) {
    1354           0 :                 ShowWarningError(state, "GetInputViewFactors: " + cCurrentModuleObject + "=\"" + EnclosureName + "\", not enough values.");
    1355           0 :                 ShowContinueError(state,
    1356           0 :                                   format("...Number of input values [{}] is less than the required number=[{}] Missing surface pairs will have a "
    1357             :                                          "zero view factor.",
    1358             :                                          NumNums,
    1359           0 :                                          pow_2(N)));
    1360             :             }
    1361             : 
    1362         202 :             for (index = 2; index <= NumAlphas; index += 2) {
    1363         200 :                 inx1 = UtilityRoutines::FindItemInList(state.dataIPShortCut->cAlphaArgs(index), enclosureSurfaceNames, N);
    1364         200 :                 inx2 = UtilityRoutines::FindItemInList(state.dataIPShortCut->cAlphaArgs(index + 1), enclosureSurfaceNames, N);
    1365         200 :                 if (inx1 == 0) {
    1366           0 :                     ShowSevereError(state, "GetInputViewFactors: " + cCurrentModuleObject + "=\"" + EnclosureName + "\", invalid surface name.");
    1367           0 :                     ShowContinueError(state, "...Surface name=\"" + state.dataIPShortCut->cAlphaArgs(index) + "\", not in this zone or enclosure.");
    1368           0 :                     ErrorsFound = true;
    1369             :                 }
    1370         200 :                 if (inx2 == 0) {
    1371           0 :                     ShowSevereError(state, "GetInputViewFactors: " + cCurrentModuleObject + "=\"" + EnclosureName + "\", invalid surface name.");
    1372           0 :                     ShowContinueError(state,
    1373           0 :                                       "...Surface name=\"" + state.dataIPShortCut->cAlphaArgs(index + 2) + "\", not in this zone or enclosure.");
    1374           0 :                     ErrorsFound = true;
    1375             :                 }
    1376         200 :                 ++numinx1;
    1377         200 :                 if (inx1 > 0 && inx2 > 0) F(inx2, inx1) = state.dataIPShortCut->rNumericArgs(numinx1);
    1378             :             }
    1379             : 
    1380           2 :             enclosureSurfaceNames.deallocate();
    1381             :         }
    1382           4 :     }
    1383             : 
    1384        9597 :     void CalcApproximateViewFactors(EnergyPlusData &state,
    1385             :                                     int const N,                    // NUMBER OF SURFACES
    1386             :                                     const Array1D<Real64> &A,       // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
    1387             :                                     const Array1D<Real64> &Azimuth, // Facing angle of the surface (in degrees)
    1388             :                                     const Array1D<Real64> &Tilt,    // Tilt angle of the surface (in degrees)
    1389             :                                     Array2A<Real64> F,              // APPROXIMATE DIRECT VIEW FACTOR MATRIX (N X N)
    1390             :                                     const Array1D_int &SPtr         // pointer to REAL(r64) surface number (for error message)
    1391             :     )
    1392             :     {
    1393             : 
    1394             :         // SUBROUTINE INFORMATION:
    1395             :         //       AUTHOR         Curt Pedersen
    1396             :         //       DATE WRITTEN   July 2000
    1397             :         //       MODIFIED       March 2001 (RKS) to disallow surfaces facing the same direction to interact radiatively
    1398             :         //                      May 2002 (COP) to include INTMASS, FLOOR, ROOF and CEILING.
    1399             :         //       RE-ENGINEERED  September 2000 (RKS for EnergyPlus)
    1400             : 
    1401             :         // PURPOSE OF THIS SUBROUTINE:
    1402             :         // This subroutine approximates view factors using an area weighting.
    1403             :         // This is improved by one degree by not allowing surfaces facing the same
    1404             :         // direction to "see" each other.
    1405             : 
    1406             :         // METHODOLOGY EMPLOYED:
    1407             :         // Each surface sees some area of other surfaces within the zone.  The view
    1408             :         // factors from the surface to the other seen surfaces are defined by their
    1409             :         // area over the summed area of seen surfaces.  Surfaces facing the same angle
    1410             :         // are assumed to not be able to see each other.
    1411             :         //  Modified May 2002 to cover poorly defined surface orientation.  Now all thermal masses, roofs and
    1412             :         //  ceilings are "seen" by other surfaces. Floors are seen by all other surfaces, but
    1413             :         //  not by other floors.
    1414             : 
    1415             :         // REFERENCES:
    1416             :         // na
    1417             : 
    1418             :         // USE STATEMENTS:
    1419             :         // na
    1420             : 
    1421             :         // Argument array dimensioning
    1422        9597 :         EP_SIZE_CHECK(A, N);
    1423        9597 :         EP_SIZE_CHECK(Azimuth, N);
    1424        9597 :         EP_SIZE_CHECK(Tilt, N);
    1425        9597 :         F.dim(N, N);
    1426        9597 :         EP_SIZE_CHECK(SPtr, N);
    1427             : 
    1428             :         // Locals
    1429             :         // SUBROUTINE ARGUMENTS:
    1430             : 
    1431             :         // SUBROUTINE PARAMETER DEFINITIONS:
    1432        9597 :         Real64 constexpr SameAngleLimit(10.0); // If the difference in the azimuth angles are above this value (degrees),
    1433             :         // then the surfaces are assumed to be facing different directions.
    1434             : 
    1435             :         // INTERFACE BLOCK SPECIFICATIONS
    1436             :         // na
    1437             : 
    1438             :         // DERIVED TYPE DEFINITIONS
    1439             :         // na
    1440             : 
    1441             :         // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    1442             : 
    1443             :         int i; // DO loop counters for surfaces in the zone
    1444             :         int j;
    1445       19194 :         Array1D<Real64> ZoneArea; // Sum of the area of all zone surfaces seen
    1446             : 
    1447             :         // Calculate the sum of the areas seen by all zone surfaces
    1448        9597 :         ZoneArea.dimension(N, 0.0);
    1449       93696 :         for (i = 1; i <= N; ++i) {
    1450      957020 :             for (j = 1; j <= N; ++j) {
    1451             :                 //  No surface sees itself and no floor sees another floor
    1452      983281 :                 if (i == j || (state.dataSurface->Surface(SPtr(j)).Class == SurfaceClass::Floor &&
    1453      110360 :                                state.dataSurface->Surface(SPtr(i)).Class == SurfaceClass::Floor))
    1454      102359 :                     continue;
    1455             :                 // All surface types see internal mass
    1456             :                 // All surface types see floors
    1457             :                 // Floors always see ceilings/roofs
    1458             :                 // All other surfaces whose tilt or facing angle differences are greater than 10 degrees see each other
    1459     2272003 :                 if ((state.dataSurface->Surface(SPtr(j)).Class == SurfaceClass::IntMass) ||
    1460     1422491 :                     (state.dataSurface->Surface(SPtr(i)).Class == SurfaceClass::IntMass) ||
    1461     1296044 :                     (state.dataSurface->Surface(SPtr(j)).Class == SurfaceClass::Floor) ||
    1462     1121684 :                     (state.dataSurface->Surface(SPtr(i)).Class == SurfaceClass::Floor) ||
    1463     1805066 :                     (std::abs(Azimuth(i) - Azimuth(j)) > SameAngleLimit && std::abs(Azimuth(i) - Azimuth(j)) < 360.0 - SameAngleLimit) ||
    1464      137938 :                     (std::abs(Tilt(i) - Tilt(j)) > SameAngleLimit)) {
    1465      670618 :                     ZoneArea(i) += A(j);
    1466             :                 }
    1467             :             }
    1468       84099 :             if (ZoneArea(i) <= 0.0) {
    1469           0 :                 ShowWarningError(state, "CalcApproximateViewFactors: Zero area for all other zone surfaces.");
    1470           0 :                 ShowContinueError(state,
    1471           0 :                                   "Happens for Surface=\"" + state.dataSurface->Surface(SPtr(i)).Name +
    1472           0 :                                       "\" in Zone=" + state.dataHeatBal->Zone(state.dataSurface->Surface(SPtr(i)).Zone).Name);
    1473             :             }
    1474             :         }
    1475             : 
    1476             :         // Set up the approximate view factors.  First these are initialized to all zero.
    1477             :         // This will clear out any junk leftover from whenever.  Then, for each zone
    1478             :         // surface, set the view factor from that surface to other surfaces as the
    1479             :         // area of the other surface divided by the sum of the area of all zone surfaces
    1480             :         // that the original surface can actually see (calculated above).  This will
    1481             :         // allow that the sum of all view factors from the original surface to all other
    1482             :         // surfaces will equal unity.  F(I,J)=0 if I=J or if the surfaces face the same
    1483             :         // direction.
    1484        9597 :         F = 0.0;
    1485       93696 :         for (i = 1; i <= N; ++i) {
    1486      957020 :             for (j = 1; j <= N; ++j) {
    1487             :                 //  No surface sees itself and no floor sees another floor
    1488      983281 :                 if (i == j || (state.dataSurface->Surface(SPtr(j)).Class == SurfaceClass::Floor &&
    1489      110360 :                                state.dataSurface->Surface(SPtr(i)).Class == SurfaceClass::Floor))
    1490      102359 :                     continue;
    1491             :                 // All surface types see internal mass
    1492             :                 // All surface types see floors
    1493             :                 // Floors always see ceilings/roofs
    1494             :                 // All other surfaces whose tilt or facing angle differences are greater than 10 degrees see each other
    1495     2272003 :                 if ((state.dataSurface->Surface(SPtr(j)).Class == SurfaceClass::IntMass) ||
    1496     1422491 :                     (state.dataSurface->Surface(SPtr(i)).Class == SurfaceClass::IntMass) ||
    1497     1296044 :                     (state.dataSurface->Surface(SPtr(j)).Class == SurfaceClass::Floor) ||
    1498     1121684 :                     (state.dataSurface->Surface(SPtr(i)).Class == SurfaceClass::Floor) ||
    1499     1805066 :                     (std::abs(Azimuth(i) - Azimuth(j)) > SameAngleLimit && std::abs(Azimuth(i) - Azimuth(j)) < 360.0 - SameAngleLimit) ||
    1500      137938 :                     (std::abs(Tilt(i) - Tilt(j)) > SameAngleLimit)) {
    1501             :                     // avoid a divide by zero if there are no other surfaces in the zone that can be seen
    1502      670618 :                     if (ZoneArea(i) > 0.0) F(j, i) = A(j) / (ZoneArea(i));
    1503             :                 }
    1504             :             }
    1505             :         }
    1506             : 
    1507        9597 :         ZoneArea.deallocate();
    1508        9597 :     }
    1509             : 
    1510        9599 :     void FixViewFactors(EnergyPlusData &state,
    1511             :                         int const N,                      // NUMBER OF SURFACES
    1512             :                         const Array1D<Real64> &A,         // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
    1513             :                         Array2A<Real64> F,                // APPROXIMATE DIRECT VIEW FACTOR MATRIX (N X N)
    1514             :                         std::string &enclName,            // Name of Enclosure being fixed
    1515             :                         std::vector<int> const spaceNums, // Zones which are part of this enclosure
    1516             :                         Real64 &OriginalCheckValue,       // check of SUM(F) - N
    1517             :                         Real64 &FixedCheckValue,          // check after fixed of SUM(F) - N
    1518             :                         Real64 &FinalCheckValue,          // the one to go with
    1519             :                         int &NumIterations,               // number of iterations to fixed
    1520             :                         Real64 &RowSum,                   // RowSum of Fixed
    1521             :                         bool const anyIntMassInZone       // are there any internal mass surfaces in the zone
    1522             :     )
    1523             :     {
    1524             : 
    1525             :         // SUBROUTINE INFORMATION:
    1526             :         //       AUTHOR         Curt Pedersen
    1527             :         //       DATE WRITTEN   July 2000
    1528             :         //       MODIFIED       September 2000 (RKS for EnergyPlus)
    1529             :         //                      April 2005,COP added capability to handle a
    1530             :         //                      surface larger than sum of all others (nonenclosure)
    1531             :         //                      by using a Fii view factor for that surface. Process is
    1532             :         //                      now much more robust and stable.
    1533             :         //       RE-ENGINEERED  na
    1534             : 
    1535             :         // PURPOSE OF THIS SUBROUTINE:
    1536             :         // This subroutine fixes approximate view factors and enforces reciprocity
    1537             :         // and completeness.
    1538             : 
    1539             :         // METHODOLOGY EMPLOYED:
    1540             :         // A(i)*F(i,j)=A(j)*F(j,i); F(i,i)=0.; SUM(F(i,j)=1.0, j=1,N)
    1541             :         // Subroutine takes approximate view factors and enforces reciprocity by
    1542             :         // averaging AiFij and AjFji.  Then it determines a set of row coefficients
    1543             :         // which can be multipled by each AF product to force the sum of AiFij for
    1544             :         // each row to equal Ai, and applies them. Completeness is checked, and if
    1545             :         // not satisfied, the AF averaging and row modifications are repeated until
    1546             :         // completeness is within a preselected small deviation from 1.0
    1547             :         // The routine also checks the number of surfaces and if N<=3, just enforces reciprocity.
    1548             : 
    1549             :         // REFERENCES:
    1550             :         // na
    1551             : 
    1552             :         // Using/Aliasing
    1553             : 
    1554             :         // Argument array dimensioning
    1555        9599 :         EP_SIZE_CHECK(A, N);
    1556        9599 :         F.dim(N, N);
    1557             : 
    1558             :         // Locals
    1559             :         // SUBROUTINE ARGUMENTS:
    1560             : 
    1561             :         // SUBROUTINE PARAMETER DEFINITIONS:
    1562        9599 :         Real64 constexpr PrimaryConvergence(0.001);
    1563        9599 :         Real64 constexpr DifferenceConvergence(0.00001);
    1564             : 
    1565             :         // INTERFACE BLOCK SPECIFICATIONS
    1566             :         // na
    1567             : 
    1568             :         // DERIVED TYPE DEFINITIONS
    1569             :         // na
    1570             : 
    1571             :         // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    1572             :         Real64 LargestArea;
    1573             :         Real64 ConvrgNew;
    1574             :         Real64 ConvrgOld;
    1575             :         Real64 Accelerator;            // RowCoefficient multipler to accelerate convergence
    1576             :         Real64 CheckConvergeTolerance; // check value for actual warning
    1577             : 
    1578             :         bool Converged;
    1579             :         int i;
    1580             :         int j;
    1581             :         bool severeErrorPresent;
    1582             :         // OriginalCheckValue is the first pass at a completeness check.  Even if this is zero,
    1583             :         // there is no guarantee that reciprocity is satisfied.  As a result, the rest of this
    1584             :         // routine is needed to correct any issues even if the user defined view factors
    1585             :         // satisfy completeness (OriginalCheckValue = 0).
    1586        9599 :         OriginalCheckValue = std::abs(sum(F) - N);
    1587             : 
    1588             :         //  Allocate and zero arrays
    1589       19196 :         Array2D<Real64> FixedAF(F); // store for largest area check
    1590             : 
    1591        9599 :         Accelerator = 1.0;
    1592        9599 :         ConvrgOld = 10.0;
    1593        9599 :         LargestArea = maxval(A);
    1594        9599 :         severeErrorPresent = false;
    1595             :         // Check for Strange Geometry
    1596             :         // When one surface has an area that exceeds the sum of all other surface areas in a zone,
    1597             :         // essentially the situation is a non-complete enclosure.  As a result, the view factors
    1598             :         // calculated using the standard procedure below will not converge and may result in invalid
    1599             :         // view factors where either reciprocity or completeness is not satisfied.  However, when
    1600             :         // the largest surface is just slightly smaller than the rest of the surface areas in the
    1601             :         // zone, it has been shown that there can still be problems.  The correction below can
    1602             :         // be helpful in avoiding these problems.  So, the criteria below (with the 0.99 term added
    1603             :         // into the comparison in the next line of code) intends to capture more cases that are
    1604             :         // "unbalanced" in surface area distribution so other strange cases can take advantage of
    1605             :         // this correction.  The use of 0.99 is simply to provide some reasonable boundary numerically
    1606             :         // and does not have some derived theoretical basis.
    1607        9599 :         if (LargestArea > 0.99 * (sum(A) - LargestArea) && (N > 3)) {
    1608         214 :             for (i = 1; i <= N; ++i) {
    1609         214 :                 if (LargestArea != A(i)) continue;
    1610          26 :                 state.dataHeatBalIntRadExchg->LargestSurf = i;
    1611          26 :                 break;
    1612             :             }
    1613          26 :             FixedAF(state.dataHeatBalIntRadExchg->LargestSurf, state.dataHeatBalIntRadExchg->LargestSurf) =
    1614          26 :                 min(0.9, 1.2 * LargestArea / sum(A)); // Give self view to big surface
    1615             :         }
    1616             : 
    1617             :         //  Set up AF matrix.
    1618       19196 :         Array2D<Real64> AF(N, N); // = (AREA * DIRECT VIEW FACTOR) MATRIX
    1619       93718 :         for (i = 1; i <= N; ++i) {
    1620      957240 :             for (j = 1; j <= N; ++j) {
    1621      873121 :                 AF(j, i) = FixedAF(j, i) * A(i);
    1622             :             }
    1623             :         }
    1624             : 
    1625             :         //  Enforce reciprocity by averaging AiFij and AjFji
    1626        9599 :         FixedAF = 0.5 * (AF + transpose(AF)); // Performance Slow way to average with transpose (heap use)
    1627             : 
    1628        9599 :         AF.deallocate();
    1629             : 
    1630       19196 :         Array2D<Real64> FixedF(N, N); // CORRECTED MATRIX OF VIEW FACTORS (N X N)
    1631             : 
    1632        9599 :         NumIterations = 0;
    1633        9599 :         RowSum = 0.0;
    1634             :         //  Check for physically unreasonable enclosures.
    1635             : 
    1636        9599 :         if (N <= 3) {
    1637           8 :             for (i = 1; i <= N; ++i) {
    1638          24 :                 for (j = 1; j <= N; ++j) {
    1639          18 :                     FixedF(j, i) = FixedAF(j, i) / A(i);
    1640             :                 }
    1641             :             }
    1642             : 
    1643           2 :             ShowWarningError(state, "Surfaces in Zone/Enclosure=\"" + enclName + "\" do not define an enclosure.");
    1644           2 :             ShowContinueError(state, "Number of surfaces <= 3, view factors are set to force reciprocity but may not fulfill completeness.");
    1645           2 :             ShowContinueError(state, "Reciprocity means that radiant exchange between two surfaces will match and not lead to an energy loss.");
    1646           2 :             ShowContinueError(state,
    1647             :                               "Completeness means that all of the view factors between a surface and the other surfaces in a zone add up to unity.");
    1648           2 :             ShowContinueError(state,
    1649             :                               "So, when there are three or less surfaces in a zone, EnergyPlus will make sure there are no losses of energy but");
    1650           2 :             ShowContinueError(
    1651             :                 state, "it will not exchange the full amount of radiation with the rest of the zone as it would if there was a completed enclosure.");
    1652             : 
    1653           2 :             RowSum = sum(FixedF);
    1654           2 :             if (RowSum > (N + 0.01)) {
    1655             :                 // Reciprocity enforced but there is more radiation than possible somewhere since the sum of one of the rows
    1656             :                 // is now greater than unity.  This should not be allowed as it can cause issues with the heat balance.
    1657             :                 // Correct this by finding the largest row summation and dividing all of the elements in the F matrix by
    1658             :                 // this max summation.  This will provide a cap on radiation so that no row has a sum greater than unity
    1659             :                 // and will still maintain reciprocity.
    1660           0 :                 Array1D<Real64> sumFixedF;
    1661             :                 Real64 MaxFixedFRowSum;
    1662           0 :                 sumFixedF.allocate(N);
    1663           0 :                 sumFixedF = 0.0;
    1664           0 :                 for (i = 1; i <= N; ++i) {
    1665           0 :                     for (j = 1; j <= N; ++j) {
    1666           0 :                         sumFixedF(i) += FixedF(i, j);
    1667             :                     }
    1668           0 :                     if (i == 1) {
    1669           0 :                         MaxFixedFRowSum = sumFixedF(i);
    1670             :                     } else {
    1671           0 :                         if (sumFixedF(i) > MaxFixedFRowSum) MaxFixedFRowSum = sumFixedF(i);
    1672             :                     }
    1673             :                 }
    1674           0 :                 sumFixedF.deallocate();
    1675           0 :                 if (MaxFixedFRowSum < 1.0) {
    1676           0 :                     ShowFatalError(state, " FixViewFactors: Three surface or less zone failing ViewFactorFix correction which should never happen.");
    1677             :                 } else {
    1678           0 :                     FixedF *= (1.0 / MaxFixedFRowSum);
    1679             :                 }
    1680           0 :                 RowSum = sum(FixedF); // needs to be recalculated
    1681             :             }
    1682           2 :             FinalCheckValue = FixedCheckValue = std::abs(RowSum - N);
    1683           2 :             F = FixedF;
    1684           4 :             for (int spaceNum : spaceNums) {
    1685           2 :                 state.dataHeatBal->Zone(state.dataHeatBal->space(spaceNum).zoneNum).EnforcedReciprocity = true;
    1686             :             }
    1687           2 :             return; // Do not iterate, stop with reciprocity satisfied.
    1688             : 
    1689             :         } //  N <= 3 Case
    1690             : 
    1691             :         //  Regular fix cases
    1692       19194 :         Array1D<Real64> RowCoefficient(N);
    1693        9597 :         Converged = false;
    1694      410025 :         while (!Converged) {
    1695      200214 :             ++NumIterations;
    1696     2076121 :             for (i = 1; i <= N; ++i) {
    1697             :                 // Determine row coefficients which will enforce closure.
    1698     1875907 :                 Real64 const sum_FixedAF_i(sum(FixedAF(_, i)));
    1699     1875907 :                 if (std::abs(sum_FixedAF_i) > 1.0e-10) {
    1700     1875907 :                     RowCoefficient(i) = A(i) / sum_FixedAF_i;
    1701             :                 } else {
    1702           0 :                     RowCoefficient(i) = 1.0;
    1703             :                 }
    1704     1875907 :                 FixedAF(_, i) *= RowCoefficient(i);
    1705             :             }
    1706             : 
    1707             :             //  Enforce reciprocity by averaging AiFij and AjFji
    1708      200214 :             FixedAF = 0.5 * (FixedAF + transpose(FixedAF));
    1709             : 
    1710             :             //  Form FixedF matrix
    1711     2076121 :             for (i = 1; i <= N; ++i) {
    1712    23606802 :                 for (j = 1; j <= N; ++j) {
    1713    21730895 :                     FixedF(j, i) = FixedAF(j, i) / A(i);
    1714    21730895 :                     if (std::abs(FixedF(j, i)) < 1.e-10) {
    1715     5671051 :                         FixedF(j, i) = 0.0;
    1716     5671051 :                         FixedAF(j, i) = 0.0;
    1717             :                     }
    1718             :                 }
    1719             :             }
    1720             : 
    1721      200214 :             ConvrgNew = std::abs(sum(FixedF) - N);
    1722      200214 :             if (std::abs(ConvrgOld - ConvrgNew) < DifferenceConvergence || ConvrgNew <= PrimaryConvergence) { //  Change in sum of Fs must be small.
    1723        9597 :                 Converged = true;
    1724             :             }
    1725      200214 :             ConvrgOld = ConvrgNew;
    1726      200214 :             if (NumIterations > 400) { //  If everything goes bad,enforce reciprocity and go home.
    1727             :                 //  Enforce reciprocity by averaging AiFij and AjFji
    1728           0 :                 FixedAF = 0.5 * (FixedAF + transpose(FixedAF));
    1729             : 
    1730             :                 //  Form FixedF matrix
    1731           0 :                 for (i = 1; i <= N; ++i) {
    1732           0 :                     for (j = 1; j <= N; ++j) {
    1733           0 :                         FixedF(j, i) = FixedAF(j, i) / A(i);
    1734             :                     }
    1735             :                 }
    1736           0 :                 Real64 const sum_FixedF(sum(FixedF));
    1737           0 :                 FinalCheckValue = FixedCheckValue = CheckConvergeTolerance = std::abs(sum_FixedF - N);
    1738           0 :                 RowSum = sum_FixedF;
    1739           0 :                 if (CheckConvergeTolerance > 0.005) {
    1740           0 :                     if (CheckConvergeTolerance > 0.1) {
    1741           0 :                         ShowSevereError(state,
    1742             :                                         "FixViewFactors: View factors convergence has failed "
    1743           0 :                                         "and will lead to heat balance errors in zone=\"" +
    1744           0 :                                             enclName + "\".");
    1745           0 :                         severeErrorPresent = true;
    1746             :                     }
    1747           0 :                     ShowWarningError(state,
    1748             :                                      "FixViewFactors: View factors not complete. Check "
    1749           0 :                                      "for bad surface descriptions or unenclosed zone=\"" +
    1750           0 :                                          enclName + "\".");
    1751           0 :                     ShowContinueError(state,
    1752           0 :                                       format("Enforced reciprocity has tolerance (ideal is "
    1753             :                                              "0)=[{:.6R}], Row Sum (ideal is {})=[{:.2R}].",
    1754             :                                              CheckConvergeTolerance,
    1755             :                                              N,
    1756           0 :                                              RowSum));
    1757           0 :                     ShowContinueError(state,
    1758             :                                       "If zone is unusual or tolerance is on the order of 0.001, view "
    1759             :                                       "factors might be OK but results should be checked carefully.");
    1760           0 :                     if (anyIntMassInZone) {
    1761           0 :                         ShowContinueError(state,
    1762             :                                           "For zones with internal mass like this one, this"
    1763             :                                           "can happen when the internal mass has an area that"
    1764             :                                           "is much larger than the other surfaces in the zone.");
    1765           0 :                         ShowContinueError(state,
    1766             :                                           "If a single thermal mass element exists in this zone"
    1767             :                                           "that has an area that is larger than the sum of the"
    1768             :                                           "rest of the surface areas, consider breaking it up"
    1769             :                                           "into two or more separate internal mass elements.");
    1770             :                     }
    1771             :                 }
    1772           0 :                 if (std::abs(FixedCheckValue) < std::abs(OriginalCheckValue)) {
    1773           0 :                     F = FixedF;
    1774           0 :                     FinalCheckValue = FixedCheckValue;
    1775             :                 }
    1776           0 :                 return;
    1777             :             }
    1778             :         }
    1779        9597 :         FixedCheckValue = ConvrgNew;
    1780        9597 :         if (FixedCheckValue < OriginalCheckValue) {
    1781           0 :             F = FixedF;
    1782           0 :             FinalCheckValue = FixedCheckValue;
    1783             :         } else {
    1784        9597 :             FinalCheckValue = OriginalCheckValue;
    1785        9597 :             RowSum = sum(FixedF);
    1786        9597 :             if (std::abs(RowSum - N) < PrimaryConvergence) {
    1787        9597 :                 F = FixedF;
    1788        9597 :                 FinalCheckValue = FixedCheckValue;
    1789             :             } else {
    1790           0 :                 ShowWarningError(
    1791           0 :                     state, "FixViewFactors: View factors not complete. Check for bad surface descriptions or unenclosed zone=\"" + enclName + "\".");
    1792             :             }
    1793             :         }
    1794        9597 :         if (severeErrorPresent) {
    1795           0 :             ShowFatalError(state,
    1796             :                            "FixViewFactors: View factor calculations significantly out "
    1797             :                            "of tolerance.  See above messages for more information.");
    1798             :         }
    1799             :     }
    1800             : 
    1801        9599 :     bool DoesZoneHaveInternalMass(EnergyPlusData &state, int const numZoneSurfaces, const Array1D_int &surfPointer)
    1802             :     {
    1803       81545 :         for (int i = 1; i <= numZoneSurfaces; ++i) {
    1804       76310 :             if (state.dataSurface->Surface(surfPointer(i)).Class == SurfaceClass::IntMass) return true;
    1805             :         }
    1806        5235 :         return false;
    1807             :     }
    1808             : 
    1809       90537 :     void CalcScriptF(EnergyPlusData &state,
    1810             :                      int const N,              // Number of surfaces
    1811             :                      Array1D<Real64> const &A, // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
    1812             :                      Array2<Real64> const &F,  // DIRECT VIEW FACTOR MATRIX (N X N)
    1813             :                      Array1D<Real64> &EMISS,   // VECTOR OF SURFACE EMISSIVITIES
    1814             :                      Array2<Real64> &ScriptF   // MATRIX OF SCRIPT F FACTORS (N X N) //Tuned Transposed
    1815             :     )
    1816             :     {
    1817             : 
    1818             :         // SUBROUTINE INFORMATION:
    1819             :         //       AUTHOR         Curt Pedersen
    1820             :         //       DATE WRITTEN   1980
    1821             :         //       MODIFIED       July 2000 (COP for the ASHRAE Loads Toolkit)
    1822             :         //       RE-ENGINEERED  September 2000 (RKS for EnergyPlus)
    1823             :         //       RE-ENGINEERED  June 2014 (Stuart Mentzer): Performance tuned
    1824             : 
    1825             :         // PURPOSE OF THIS SUBROUTINE:
    1826             :         // Determines Hottel's ScriptF coefficients which account for the total
    1827             :         // grey interchange between surfaces in an enclosure.
    1828             : 
    1829             :         // METHODOLOGY EMPLOYED:
    1830             :         // See reference
    1831             : 
    1832             :         // REFERENCES:
    1833             :         // Hottel, H. C. and A. F. Sarofim, Radiative Transfer, Ch 3, McGraw Hill, 1967.
    1834             : 
    1835             :         // USE STATEMENTS:
    1836             :         // na
    1837             : 
    1838             :         // Locals
    1839             :         // SUBROUTINE ARGUMENTS:
    1840             :         // --Must satisfy reciprocity and completeness:
    1841             :         //  A(i)*F(i,j)=A(j)*F(j,i); F(i,i)=0.; SUM(F(i,j)=1.0, j=1,N)
    1842             : 
    1843             :         // SUBROUTINE PARAMETER DEFINITIONS:
    1844       90537 :         Real64 constexpr MaxEmissLimit(0.99999); // Limit the emissivity internally/avoid a divide by zero error
    1845             : 
    1846             :         // INTERFACE BLOCK SPECIFICATIONS
    1847             :         // na
    1848             : 
    1849             :         // DERIVED TYPE DEFINITIONS
    1850             :         // na
    1851             : 
    1852             :         // Validate argument array dimensions
    1853       90537 :         assert(N >= 0); // Do we need to allow for N==0?
    1854       90537 :         assert((A.l() == 1) && (A.u() == N));
    1855       90537 :         assert((F.l1() == 1) && (F.u1() == N));
    1856       90537 :         assert((F.l2() == 1) && (F.u2() == N));
    1857       90537 :         assert((EMISS.l() == 1) && (EMISS.u() == N));
    1858       90537 :         assert(equal_dimensions(F, ScriptF));
    1859             : 
    1860             :         // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
    1861             : 
    1862             : #ifdef EP_Count_Calls
    1863             :         ++state.dataTimingsData->NumCalcScriptF_Calls;
    1864             : #endif
    1865             : 
    1866             :         // Load Cmatrix with AF (AREA * DIRECT VIEW FACTOR) matrix
    1867      181074 :         Array2D<Real64> Cmatrix(N, N);        // = (AF - EMISS/REFLECTANCE) matrix (but plays other roles)
    1868       90537 :         assert(equal_dimensions(Cmatrix, F)); // For linear indexing
    1869       90537 :         Array2D<Real64>::size_type l(0u);
    1870      838133 :         for (int j = 1; j <= N; ++j) {
    1871     7860912 :             for (int i = 1; i <= N; ++i, ++l) {
    1872     7113316 :                 Cmatrix[l] = A(i) * F[l]; // [ l ] == ( i, j )
    1873             :             }
    1874             :         }
    1875             : 
    1876             :         // Load Cmatrix with (AF - EMISS/REFLECTANCE) matrix
    1877      181074 :         Array1D<Real64> Excite(N); // Excitation vector = A*EMISS/REFLECTANCE
    1878       90537 :         l = 0u;
    1879      838133 :         for (int i = 1; i <= N; ++i, l += N + 1) {
    1880      747596 :             Real64 EMISS_i(EMISS(i));
    1881      747596 :             if (EMISS_i > MaxEmissLimit) { // Check/limit EMISS for this surface to avoid divide by zero below
    1882           0 :                 EMISS_i = EMISS(i) = MaxEmissLimit;
    1883           0 :                 ShowWarningError(state, "A thermal emissivity above 0.99999 was detected. This is not allowed. Value was reset to 0.99999");
    1884             :             }
    1885      747596 :             Real64 const EMISS_i_fac(A(i) / (1.0 - EMISS_i));
    1886      747596 :             Excite(i) = -EMISS_i * EMISS_i_fac; // Set up matrix columns for partial radiosity calculation
    1887      747596 :             Cmatrix[l] -= EMISS_i_fac;          // Coefficient matrix for partial radiosity calculation // [ l ] == ( i, i )
    1888             :         }
    1889             : 
    1890      181074 :         Array2D<Real64> Cinverse(N, N);       // Inverse of Cmatrix
    1891       90537 :         CalcMatrixInverse(Cmatrix, Cinverse); // SOLVE THE LINEAR SYSTEM
    1892       90537 :         Cmatrix.clear();                      // Release memory ASAP
    1893             : 
    1894             :         // Scale Cinverse colums by excitation to get partial radiosity matrix
    1895       90537 :         l = 0u;
    1896      838133 :         for (int j = 1; j <= N; ++j) {
    1897      747596 :             Real64 const e_j(Excite(j));
    1898     7860912 :             for (int i = 1; i <= N; ++i, ++l) {
    1899     7113316 :                 Cinverse[l] *= e_j; // [ l ] == ( i, j )
    1900             :             }
    1901             :         }
    1902       90537 :         Excite.clear(); // Release memory ASAP
    1903             : 
    1904             :         // Form Script F matrix transposed
    1905       90537 :         assert(equal_dimensions(Cinverse, ScriptF)); // For linear indexing
    1906       90537 :         Array2D<Real64>::size_type m(0u);
    1907      838133 :         for (int i = 1; i <= N; ++i) { // Inefficient order for cache but can reuse multiplier so faster choice depends on N
    1908      747596 :             Real64 const EMISS_i(EMISS(i));
    1909      747596 :             Real64 const EMISS_fac(EMISS_i / (1.0 - EMISS_i));
    1910      747596 :             l = static_cast<Array2D<Real64>::size_type>(i - 1);
    1911     7860912 :             for (int j = 1; j <= N; ++j, l += N, ++m) {
    1912     7113316 :                 if (i == j) {
    1913             :                     //        ScriptF(I,J) = EMISS(I)/(1.0d0-EMISS(I))*(Jmatrix(I,J)-Delta*EMISS(I)), where Delta=1
    1914      747596 :                     ScriptF[m] = EMISS_fac * (Cinverse[l] - EMISS_i); // [ l ] = ( i, j ), [ m ] == ( j, i )
    1915             :                 } else {
    1916             :                     //        ScriptF(I,J) = EMISS(I)/(1.0d0-EMISS(I))*(Jmatrix(I,J)-Delta*EMISS(I)), where Delta=0
    1917     6365720 :                     ScriptF[m] = EMISS_fac * Cinverse[l]; // [ l ] == ( i, j ), [ m ] == ( j, i )
    1918             :                 }
    1919             :             }
    1920             :         }
    1921       90537 :     }
    1922             : 
    1923       90537 :     void CalcMatrixInverse(Array2<Real64> &A, // Matrix: Gets reduced to L\U form
    1924             :                            Array2<Real64> &I  // Returned as inverse matrix
    1925             :     )
    1926             :     {
    1927             :         // SUBROUTINE INFORMATION:
    1928             :         //       AUTHOR         Jakob Asmundsson
    1929             :         //       DATE WRITTEN   January 1999
    1930             :         //       MODIFIED       September 2000 (RKS for EnergyPlus)
    1931             :         //       RE-ENGINEERED  June 2014 (Stuart Mentzer): Performance/memory tuning rewrite
    1932             : 
    1933             :         // PURPOSE OF THIS SUBROUTINE:
    1934             :         // To find the inverse of Matrix, using partial pivoting.
    1935             : 
    1936             :         // METHODOLOGY EMPLOYED:
    1937             :         // Inverse is found using partial pivoting and Gauss elimination
    1938             : 
    1939             :         // REFERENCES:
    1940             :         // Any Linear Algebra book
    1941             : 
    1942             :         // Validation
    1943       90537 :         assert(A.square());
    1944       90537 :         assert(A.I1() == A.I2());
    1945       90537 :         assert(equal_dimensions(A, I));
    1946             : 
    1947             :         // Initialization
    1948       90537 :         int const l(A.l1());
    1949       90537 :         int const u(A.u1());
    1950       90537 :         int const n(u - l + 1);
    1951       90537 :         I.to_identity(); // I starts out as identity
    1952             : 
    1953             :         // Could do row scaling here to improve condition and then check min pivot isn't too small
    1954             : 
    1955             :         // Compute in-place LU decomposition of [A|I] with row pivoting
    1956      838133 :         for (int i = l; i <= u; ++i) {
    1957             : 
    1958             :             // Find pivot row in column i below diagonal
    1959      747596 :             int iPiv = i;
    1960      747596 :             Real64 aPiv(std::abs(A(i, i)));
    1961      747596 :             auto ik(A.index(i, i + 1));
    1962     3930456 :             for (int k = i + 1; k <= u; ++k, ++ik) {
    1963     3182860 :                 Real64 const aAki(std::abs(A[ik])); // [ ik ] == ( i, k )
    1964     3182860 :                 if (aAki > aPiv) {
    1965           0 :                     iPiv = k;
    1966           0 :                     aPiv = aAki;
    1967             :                 }
    1968             :             }
    1969      747596 :             assert(aPiv != 0.0); //? Is zero pivot possible for some user inputs? If so if test/handler needed
    1970             : 
    1971             :             // Swap row i with pivot row
    1972      747596 :             if (iPiv != i) {
    1973           0 :                 auto ji(A.index(l, i));    // [ ji ] == ( j, i )
    1974           0 :                 auto pj(A.index(l, iPiv)); // [ pj ] == ( j, iPiv )
    1975           0 :                 for (int j = l; j <= u; ++j, ji += n, pj += n) {
    1976           0 :                     Real64 const Aij(A[ji]);
    1977           0 :                     A[ji] = A[pj];
    1978           0 :                     A[pj] = Aij;
    1979           0 :                     Real64 const Iij(I[ji]);
    1980           0 :                     I[ji] = I[pj];
    1981           0 :                     I[pj] = Iij;
    1982             :                 }
    1983             :             }
    1984             : 
    1985             :             // Put multipliers in column i and reduce block below A(i,i)
    1986      747596 :             Real64 const Aii_inv(1.0 / A(i, i));
    1987     3930456 :             for (int k = i + 1; k <= u; ++k) {
    1988     3182860 :                 Real64 const multiplier(A(i, k) * Aii_inv);
    1989     3182860 :                 A(i, k) = multiplier;
    1990     3182860 :                 if (multiplier != 0.0) {
    1991     3094162 :                     auto ji(A.index(i + 1, i)); // [ ji ] == ( j, i )
    1992     3094162 :                     auto jk(A.index(i + 1, k)); // [ jk ] == ( j, k )
    1993    31498602 :                     for (int j = i + 1; j <= u; ++j, ji += n, jk += n) {
    1994    28404440 :                         A[jk] -= multiplier * A[ji];
    1995             :                     }
    1996     3094162 :                     ji = A.index(l, i);
    1997     3094162 :                     jk = A.index(l, k);
    1998    47602654 :                     for (int j = l; j <= u; ++j, ji += n, jk += n) {
    1999    44508492 :                         Real64 const Iij(I[ji]);
    2000    44508492 :                         if (Iij != 0.0) {
    2001    16057486 :                             I[jk] -= multiplier * Iij;
    2002             :                         }
    2003             :                     }
    2004             :                 }
    2005             :             }
    2006             :         }
    2007             : 
    2008             :         // Perform back-substitution on [U|I] to put inverse in I
    2009      838133 :         for (int k = u; k >= l; --k) {
    2010      747596 :             Real64 const Akk_inv(1.0 / A(k, k));
    2011      747596 :             auto jk(A.index(l, k)); // [ jk ] == ( j, k )
    2012     7860912 :             for (int j = l; j <= u; ++j, jk += n) {
    2013     7113316 :                 I[jk] *= Akk_inv;
    2014             :             }
    2015      747596 :             auto ik(A.index(k, l));             // [ ik ] == ( i, k )
    2016     3930456 :             for (int i = l; i < k; ++i, ++ik) { // Eliminate kth column entries from I in rows above k
    2017     3182860 :                 Real64 const Aik(A[ik]);
    2018     3182860 :                 auto ji(A.index(l, i)); // [ ji ] == ( j, i )
    2019     3182860 :                 auto jk(A.index(l, k)); // [ jk ] == ( k, j )
    2020    48637656 :                 for (int j = l; j <= u; ++j, ji += n, jk += n) {
    2021    45454796 :                     I[ji] -= Aik * I[jk];
    2022             :                 }
    2023             :             }
    2024             :         }
    2025       90537 :     }
    2026             : 
    2027          19 :     void CalcFMRT(EnergyPlusData &state,
    2028             :                   int const N,              // Number of surfaces
    2029             :                   Array1D<Real64> const &A, // AREA VECTOR- ASSUMED,BE N ELEMENTS LONG
    2030             :                   Array1D<Real64> &FMRT     // VECTOR OF MEAN RADIANT TEMPERATURE "VIEW FACTORS"
    2031             :     )
    2032             :     {
    2033          19 :         double sumAF = 0.0;
    2034         177 :         for (int iS = 0; iS < N; iS++) {
    2035         158 :             FMRT[iS] = 1.0;
    2036         158 :             sumAF += A[iS];
    2037             :         }
    2038             : 
    2039          19 :         constexpr int maxIt = 100;
    2040          19 :         constexpr double tol = 0.0001;
    2041             :         double fChange, fLast;
    2042          19 :         double sumAFNew = sumAF;
    2043         687 :         for (unsigned i = 0; i < maxIt; i++) {
    2044         687 :             fChange = 0.;
    2045         687 :             bool errorsFound(false);
    2046         687 :             sumAF = sumAFNew;
    2047         687 :             sumAFNew = 0.0;
    2048        6726 :             for (int iS = 0; iS < N; iS++) {
    2049        6039 :                 fLast = FMRT[iS];
    2050        6039 :                 FMRT[iS] = 1. / (1. - A[iS] * FMRT[iS] / (sumAF));
    2051        6039 :                 if (FMRT[iS] > 100.) {
    2052           0 :                     errorsFound = true;
    2053           0 :                     ShowSevereError(state, "Geometry not compatible with Carroll MRT Zone Radiant Exchange method.");
    2054           0 :                     break;
    2055             :                 }
    2056        6039 :                 fChange += fabs(FMRT[iS] - fLast);
    2057        6039 :                 sumAFNew += A[iS] * FMRT[iS];
    2058             :             }
    2059             : 
    2060         687 :             if (errorsFound || fChange / N < tol) {
    2061             :                 break;
    2062             :             }
    2063         668 :             if (i >= maxIt) {
    2064           0 :                 errorsFound = true;
    2065           0 :                 ShowSevereError(state, "Carroll MRT Zone Radiant Exchange method unable to converge on \"view factor\" calculation.");
    2066             :             }
    2067         668 :             if (errorsFound) {
    2068           0 :                 ShowFatalError(state, "CalcFMRT: Errors found while calculating mean radiant temperature view factors.  Program terminated.");
    2069             :             }
    2070             :         }
    2071          19 :         return;
    2072             :     }
    2073             : 
    2074         475 :     void CalcFp(int const N,            // Number of surfaces
    2075             :                 Array1D<Real64> &EMISS, // VECTOR OF SURFACE EMISSIVITIES
    2076             :                 Array1D<Real64> &FMRT,  // VECTOR OF MEAN RADIANT TEMPERATURE "VIEW FACTORS"
    2077             :                 Array1D<Real64> &Fp     // VECTOR OF OPPENHEIM RESISTANCE VALUES
    2078             :     )
    2079             :     {
    2080         475 :         Real64 SB = DataGlobalConstants::StefanBoltzmann;
    2081        4425 :         for (int iS = 0; iS < N; iS++) {
    2082        3950 :             Fp[iS] = SB * EMISS[iS] / (EMISS[iS] / FMRT[iS] + 1. - EMISS[iS]); // actually sigma *
    2083             :         }
    2084         475 :         return;
    2085             :     }
    2086             : 
    2087         113 :     int GetRadiantSystemSurface(EnergyPlusData &state,
    2088             :                                 std::string const &cCurrentModuleObject, // Calling Object type
    2089             :                                 std::string const &RadSysName,           // Calling Object name
    2090             :                                 int const RadSysZoneNum,                 // Radiant system zone number
    2091             :                                 std::string const &SurfaceName,          // Referenced surface name
    2092             :                                 bool &ErrorsFound                        // True when errors are found
    2093             :     )
    2094             :     {
    2095             :         static constexpr std::string_view routineName("GetRadiantSystemSurface: "); // include trailing blank space
    2096             : 
    2097             :         // For radiant zone equipment, find the referenced surface and check if it is in the same zone or radiant enclosure
    2098         113 :         int const surfNum = UtilityRoutines::FindItemInList(SurfaceName, state.dataSurface->Surface);
    2099             : 
    2100             :         // Trap for surfaces that do not exist
    2101         113 :         if (surfNum == 0) {
    2102           0 :             ShowSevereError(state, std::string{routineName} + "Invalid Surface name = " + SurfaceName);
    2103           0 :             ShowContinueError(state, "Occurs for " + cCurrentModuleObject + " = " + RadSysName);
    2104           0 :             ErrorsFound = true;
    2105           0 :             return surfNum;
    2106             :         }
    2107             : 
    2108         113 :         if (RadSysZoneNum == 0) {
    2109           0 :             ShowSevereError(state, std::string{routineName} + "Invalid Zone number passed by " + cCurrentModuleObject + " = " + RadSysName);
    2110           0 :             ErrorsFound = true;
    2111           0 :             return surfNum;
    2112             :         }
    2113             : 
    2114             :         // Check if the surface and equipment are in the same zone
    2115         113 :         int const surfZoneNum = state.dataSurface->Surface(surfNum).Zone;
    2116         113 :         if (RadSysZoneNum == 0) {
    2117             :             // This should never happen - but it does in some simple unit tests that are designed to throw errors
    2118           0 :             ShowSevereError(
    2119           0 :                 state, std::string{routineName} + "Somehow the radiant system zone number is zero for" + cCurrentModuleObject + " = " + RadSysName);
    2120           0 :             ErrorsFound = true;
    2121         113 :         } else if (surfZoneNum == 0) {
    2122             :             // This should never happen
    2123           0 :             ShowSevereError(state,
    2124           0 :                             std::string{routineName} + "Somehow  the surface zone number is zero for" + cCurrentModuleObject + " = " + RadSysName +
    2125             :                                 " and Surface = " + SurfaceName); // LCOV_EXCL_LINE
    2126             :             ErrorsFound = true;                                   // LCOV_EXCL_LINE
    2127         113 :         } else if (surfZoneNum != RadSysZoneNum) {
    2128           0 :             ShowSevereError(state, std::string(routineName) + "Surface = " + SurfaceName + " is not in the same zone  as the radiant equipment.");
    2129           0 :             ShowContinueError(state, "Surface zone or enclosure = " + state.dataHeatBal->Zone(surfZoneNum).Name);
    2130           0 :             ShowContinueError(state, "Radiant equipment zone or enclosure = " + state.dataHeatBal->Zone(RadSysZoneNum).Name);
    2131           0 :             ShowContinueError(state, "Occurs for " + cCurrentModuleObject + " = " + RadSysName);
    2132           0 :             ErrorsFound = true;
    2133             :         }
    2134         113 :         return surfNum;
    2135             :     }
    2136             : 
    2137             : } // namespace HeatBalanceIntRadExchange
    2138             : 
    2139        2313 : } // namespace EnergyPlus

Generated by: LCOV version 1.13