' Dick van Dijk, July 2013
  ' Dick van Dijk, March/April 2014
  ' EPB-research@dickvandijk.nl
  ' Major contributions by Pim van den Helm
  '
  ' Further refinements
  ' June 2014 - April 2015
  ' Revision: February 2017
  ' Dick van Dijk
  ' Continued revision: March 2018
  '  reverse numbering of nodes (PLI) in elements (ELI):
  '    was from inside to outside
  '    changed into: from outside to inside
  '  Identifier for these changes: #Pl#
  '
Option Explicit
 '
 ' Global variables:
 Dim MatA(1 To 51, 1 To 51) As Double
 Dim VecB(1 To 51, 1 To 3) As Double
 Dim colB_H, colB_C, colB_act As Integer ' Used to indicate which col is H and which is C
 '                                          and which is actually used in case of recalc
 Dim nrHCmodes As Integer        ' Used to indicate the HC calculation mode
 Dim Theta_old(1 To 51) As Double     ' Temp.vector at previous timestep
 Dim XTest, XTest1, XTest2, XTemp As Variant
 Dim ElnMax, PlnMax As Integer  'absolute max nr of elements and number of layers
 Dim Eln As Integer     'Number of actual building elements
 Dim Pln(1 To 20) As Integer     'Number of actual layers in elements, per building element,
 Dim PlnSum(1 To 20) As Integer     'Nr of layers summed over all lower numbered elements
 Dim Msg As String
 '
 

Sub FillMatrix()
  ' FillMatrix Macro
  ' Actually: for time being one big routine, including time loop
  ' Completes the matrix of temperatures with time independent and time dependent data
  ' (see Sheet Explanation)
  '
  '
 Dim RunDate0, RunDate1, RunTime0, RunTime1 As Variant
 
 Dim R00, REl00, CEl00, RTstep00, RTstepClim As Integer          'Specific starting positions on sheets
 Dim ri, ci, ck As Integer       'Counter for rows and columns matrix
 Dim Rn, Cn As Integer       'Matrix dimension
 Dim Eli, Elk, Pli, Plk As Integer     'Counter for building elements, layers in elements
 Dim Tstepi, Tstep1, Tstepn, Rstep As Integer    'time steps, row
  '
  '
 Dim i As Integer
 Dim j As Integer
 Dim n As Integer
 Dim m As Integer
 Dim cBi As Integer
 '
 '
 Dim Input_t As Variant
 Dim Output_t As Variant
 Dim Month, Week, Day, Hour As Integer 'Not all used yet (Month is used for gorund floor)
  '
 Dim f_int_c As Single
 Dim f_sol_c As Single
 Dim f_HC_c As Single
 Dim f_H_c As Single
 Dim f_C_c As Single
 Dim Delta_Theta_er As Single
  '
 Dim H_tb As Variant ' Thermal bridges
 Dim Name_eli(1 To 20) As String     ' but cannot read a row without loop?
 Dim Type_eli(1 To 20) As String
 Dim TypeSub_eli(1 To 20) As String
 Dim TypeClass_eli(1 To 20) As String
 Dim U_eli, R_c_eli As Variant
 Dim kappa_m_eli As Variant  ' Not used for VB calc, but maybe for reporting
 Dim A_eli, A_eli_tot As Variant
 Dim a_sol_eli As Variant    ' For monthly energy balance
 Dim g_w_eli As Variant      ' Maybe for monthly energy balance
 Dim Or_eli, Or_clima As Variant ' Not used for VB calc, but for testing if available in ClimData sheet
 '
 Dim h_ci_eli As Variant
 Dim Ah_ci As Variant        ' Weighted sum A_eli * h_ci
 Dim h_ri_eli, h_ri_eli_mn As Variant
 Dim h_ce_eli As Variant
 Dim h_re_eli As Variant
 Dim F_sk_eli As Variant
 Dim DeltaTheta_er As Variant
 Dim h_pli_eli As Variant
 Dim kappa_pli_eli As Variant
 Dim a_sol_pli_eli As Variant
 Dim tau_sol_eli As Variant
 Dim g_sol_eli As Variant  'For any bldng element; for monthly energy balance; don't confuse with g_w_eli
 '                          (but.. actually the same; except for DYNAMIC elements?)
 Dim Phi_sol_zi As Double 'For direct solar into zone
 '
 Dim Dtime As Single
 Dim C_int, A_use As Single
 
 Dim Theta_int_air(1 To 3), Theta_int_r_mn(1 To 3), Theta_int_op(1 To 3) As Double
 Dim R_gr_ve, Theta_gr_ve(1 To 12) As Double
 Dim Theta_H_set, Theta_C_set, Theta_op_set, Theta_op_act As Double
 Dim Phi_HC_nd, Phi_HC_nd_calc(1 To 3), Phi_HC_nd_act As Double
 Dim Phi_H_nd_max, Phi_C_nd_max, Phi_H_nd_max_act, Phi_C_nd_max_act As Single
 
  '
  'Get time-independent input data
  ' Some of the following variables (see comments at declarations)
  ' are not used for VB calc, but could be used for reporting
  '
  ' Attention: Switching between worksheets makes it slow
  ' (but of course only within time and/or iteration loop)!
  
  ' Enable recalculaion of all open worksheets during initialization before time loop starts
  '
  Sheets("INFO").EnableCalculation = True
  Sheets("Explanation").EnableCalculation = True
  
  Sheets("Method_input").EnableCalculation = True
  Sheets("Input_GrFl").EnableCalculation = True
  Sheets("Calc_GrFl").EnableCalculation = True
  
  Sheets("Method_calculation").EnableCalculation = True
  Sheets("Method_output").EnableCalculation = True
  
  Sheets("Input_0").EnableCalculation = True
  Sheets("Input_p").EnableCalculation = True
  Sheets("Input_t").EnableCalculation = True
  
  Sheets("Calc_H_m").EnableCalculation = True
  Sheets("Calc_C_m").EnableCalculation = True
  
  Sheets("Graphs").EnableCalculation = True
  Sheets("Clipboard").EnableCalculation = True
  
  Sheets("Output_m").EnableCalculation = True
  Sheets("Output_t").EnableCalculation = True
  
  Sheets("ClimDat_m").EnableCalculation = True
  Sheets("ClimData").EnableCalculation = True
  
  Sheets("VCS_special").EnableCalculation = True
 
  Application.Calculation = xlCalculationAutomatic
  
  '
  ' Record calculation starting date and time
  RunDate0 = Date
  RunTime0 = Time
  Worksheets("Output_m").Range("b1").Value = "Calculation started at: " & RunDate0 & " " & RunTime0
  Worksheets("Output_m").Range("b2").Value = "MACRO is running; please be patient (Full year calculation takes typically less than 10 minutes)"
  Sheets("Output_m").Select         ' To show this sheet at start of calculation
  'Cells(1, 1).Select
  Range("A1").Select
  Application.ScreenUpdating = False 'to keep this sheet on screen
'
'
'
'  Fixed Values:
    ElnMax = 20
    PlnMax = 5
'
  With Sheets("Input_0")
    Delta_Theta_er = .Cells(6, 4).Value
    A_use = .Cells(12, 4).Value
    C_int = .Cells(13, 8).Value
    f_int_c = .Cells(14, 4).Value
    f_sol_c = .Cells(15, 4).Value
    f_H_c = .Cells(16, 4).Value
    f_C_c = .Cells(17, 4).Value
    '
    Phi_H_nd_max = .Cells(18, 4).Value
    Phi_C_nd_max = -.Cells(19, 4).Value 'NOTE: Minus sign introduced!!!
    '
    ' First row for input data per element
    '
    H_tb = .Cells(21, 4).Value
    REl00 = 29 'Third character is a L!
    
    Eln = 0
    For Eli = 1 To ElnMax ' Could this be done without loop? But now we need loop for nr of elements
      ri = REl00
      Name_eli(Eli) = .Cells(ri, 4 + Eli).Value
      ri = ri + 1
      TypeSub_eli(Eli) = .Cells(ri, 4 + Eli).Value
      If TypeSub_eli(Eli) <> "OP" And TypeSub_eli(Eli) <> "W" And _
         TypeSub_eli(Eli) <> "GR" And TypeSub_eli(Eli) <> "AD" And _
         TypeSub_eli(Eli) <> "NA" Then
            ' Stop, because number of nodes exceeds max
            Msg = MsgBox("Subtype of building element " & Eli & " in sheet Input_0 not recognized; should be OP, W, GR, AD or NA. Calculation stopped")
            Exit Sub
      End If
      If TypeSub_eli(Eli) = "OP" Then Type_eli(Eli) = "EXT" ' Consistent with Input_0 sheet (EXT=1)
      If TypeSub_eli(Eli) = "W" Then Type_eli(Eli) = "EXT"   ' Consistent with Input_0 sheet (EXT=1)
      If TypeSub_eli(Eli) = "GR" Then Type_eli(Eli) = "GR"  ' Maybe not necessary, but doesn't harm
      If TypeSub_eli(Eli) = "AD" Then Type_eli(Eli) = "AD"  ' Added March 19, 2018. Until now: no effect, because "AD" had nothing to do
      If TypeSub_eli(Eli) = "NA" Then Type_eli(Eli) = "NA"
      
      
      ri = ri + 1
      TypeClass_eli(Eli) = .Cells(ri, 4 + Eli).Value
      ' If TypeSub is not "W" (window), then the impact of the TypeClass
      ' is already taken care of in the Input_0 sheet
      If TypeSub_eli(Eli) = "W" Then
      'set default in case of wrong input:
        If TypeClass_eli(Eli) <> "F" And TypeClass_eli(Eli) <> "V1" And TypeClass_eli(Eli) <> "V2" Then TypeClass_eli(Eli) = "F"
      End If
    ' Check for last actual element (Eln):
      If (Type_eli(Eli) <> "NA") Then
        Eln = Eln + 1
      Else
        GoTo 20
      End If
    Next Eli

20  CEl00 = 5  ' l is not one but letter el
    ri = ri + 1
    ' This U-value is only input at the sheet
    ri = ri + 1
     ' This R-value is only input at the sheet
    ri = ri + 1
    kappa_m_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    ' kappa_m_eli is also not used, but no harm to get the values
    ri = ri + 1
    A_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    A_eli_tot = .Cells(ri, CEl00 + ElnMax).Value
    ri = ri + 1
    a_sol_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    ' a_sol_eli is  not used, but no harm to get the values
    ri = ri + 1
    g_w_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    ' g_w_eli is  not used, but no harm to get the values
    ri = ri + 1
    Or_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    ' Or_eli is not needed in this macro calc, but used to test if present in ClimaData sheet
    ri = ri + 1
    'Empty row
    ri = ri + 1
    h_ci_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    Ah_ci = .Cells(ri, CEl00 + ElnMax).Value
    
    ri = ri + 1
    h_ri_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    h_ri_eli_mn = .Cells(ri, CEl00 + ElnMax).Value
    ri = ri + 1
    h_ce_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    ri = ri + 1
    h_re_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    ri = ri + 1
    F_sk_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    'Solar shading obstacles select
    'not needed in this Macro, so skip:
    ri = ri + 5
    'Monthly input data ground floor:
    ri = ri + 2
    For Month = 1 To 12
    ri = ri + 1
       Theta_gr_ve(Month) = .Cells(ri, 3).Value
       If Month = 6 Then R_gr_ve = .Cells(ri, 8).Value 'Month 6 is just for the position of the cell with the value
    Next Month
    ri = ri + 2
    U_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address) ' Used for monthly balance
    ri = ri + 1
    R_c_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address) ' Maybe not used in macro (yet)?
    ri = ri + 1
    h_pli_eli = .Range(Cells(ri, CEl00).Address, Cells(ri + 3, CEl00 + Eln - 1).Address)
    ri = ri + 4
    kappa_pli_eli = .Range(Cells(ri, CEl00).Address, Cells(ri + 4, CEl00 + Eln - 1).Address)
    ri = ri + 5
    a_sol_pli_eli = .Range(Cells(ri, CEl00).Address, Cells(ri + 4, CEl00 + Eln - 1).Address)
    ri = ri + 5
    tau_sol_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    ri = ri + 1
    For Eli = 1 To ElnMax ' Could this be done without loop? But now we need loop for nr of elements
       Pln(Eli) = .Cells(ri, 4 + Eli).Value
    Next Eli
    '
    'For monthly energy balance:
    ri = ri + 5
    g_sol_eli = .Range(Cells(ri, CEl00).Address, Cells(ri, CEl00 + Eln - 1).Address)
    '
  End With         ' Sheets("Input_0")
  '
  ' Other derived variables:
  '
  ' Initial start: assume three modes needed (no H or C, H, C)
  ' (But will be overruled at start of timestep)
  nrHCmodes = 3
  colB_act = 1 ' Initial value of column used for actual calculation
  '
  '
  ' Looping structure to look at array.
  '      For i = 1 To UBound(MatA, 1)
  '      For j = 1 To UBound(MatA, 2)
  '        MsgBox "MatA(" & Format(i, "0") & "," & Format(j, "0") & ")=" & MatA(i, j)
  '      Next j
  '      Next i
  '
  ' INITIALIZATION OF FIXED VALUES:
  Rn = 1
    For Eli = 1 To Eln
      PlnSum(Eli) = Rn - 1
      Rn = Rn + Pln(Eli)
  Next Eli
  Cn = Rn          'Square matrix
  If Rn > 51 Then
  ' Stop, because number of nodes exceeds max
    Msg = MsgBox("Number of nodes exceeds max = 51. Calculation stopped")
    Exit Sub
  End If
  
  '
  'Position of cell preceding upper left corner of the matrix at the sheet
  R00 = 5          ' Needed for reading and writing (other) output to sheets
  '
  ' Position of cell preceding timeseries on the sheet(s)
  RTstep00 = 6
  RTstepClim = 14 'extra preceding rows in CLimData compared to Input_t and Output_t
  '
  'Test if all requested OrTilt angles are available in ClimaData sheet:
  'Max 20 angles, each 2 columns:
  With Sheets("ClimData")
    Or_clima = .Range(Cells(16, 21).Address, Cells(16, (21 - 1) + (2 * 20)).Address)
    For Eli = 1 To Eln
     XTest = 0
     For ci = 1 To 40
        If (Or_clima(1, ci) = Or_eli(1, Eli)) Then
          XTest = 1 'Means: match found
        End If
     Next ci
       If XTest = 0 Then
          MsgBox ("Error: requested orientation/tilt with id = " & Or_eli(1, Eli) & " not found in ClimData sheet")
          Stop
       End If
    Next Eli
  End With
  '
  ' Set initial values
  Dtime = 0
  Month = 0
  Week = 0
  'TEST:
  XTest = 0
  XTest1 = 0
  XTest2 = 0
  '
  ' Start to set initial values for the node temperatures
  ' Initial values temperatures previous timestep
  ' FOR TIME BEING: fixed value
  For ri = 1 To Rn
   Theta_old(ri) = 20#
   For cBi = 1 To 3
  ' Initial values temperatures each node
  ' (Only needed in case some conditions are temperature dependent
  '  while building up the heat balance)
     VecB(ri, cBi) = 20#
   Next cBi
  Next ri                  ' Ri = 1 To Rn
  '
  '
  'Reset possible output previous calculation
  'With Sheets("Output_t") '
  ' ATTENTION!!:
  ' The following loop (ignore the outdated cell numbers) would take > 2 minutes !!!
  ' For Tstepi = 1 To 8760 ' full year
  '  .Cells(RTstep00 + Tstepi, 16).Value = 0# ' HC needs
  '  .Cells(RTstep00 + Tstepi, 17).Value = 0# 'Indoor air temp.
  ' Next Tstepi
  'End With
  'But the following way takes less than a second!!:
  Sheets("Output_t").Select  'Funny that Select is needed!
  'Note: Dec. added before Jan., hence 24*31 hours extra
  Sheets("Output_t").Range(Cells(RTstep00 + 1, 29), Cells(RTstep00 + 8760 + 24 * 31, 33)).Value = 0#
  '
  ' Length of time period
  '
  Tstep1 = Sheets("Input_0").Cells(7, 10).Value
  Tstepn = Sheets("Input_0").Cells(8, 10).Value
  '

 ' Disable recalculation of all open workbooks while time loop runs
 ' CHECK DYNAMIC!  This assumes there is no need for recalc due to writing at sheets during time loop.
 ' If events are built in triggered by specific conditions, then a recalc may be needed?!
' Only Output_t matters so far, because of output added each time step
  
  
 ' !?! If the following list is deleted, the calculation time appears 5 times higher...
 
  Sheets("INFO").EnableCalculation = False
  Sheets("Explanation").EnableCalculation = False
  
  Sheets("Method_input").EnableCalculation = False
  Sheets("Input_GrFl").EnableCalculation = False
  Sheets("Calc_GrFl").EnableCalculation = False
  
  Sheets("Method_calculation").EnableCalculation = False
  Sheets("Method_output").EnableCalculation = False
    
  Sheets("Input_0").EnableCalculation = False
  Sheets("Input_p").EnableCalculation = False
  Sheets("Input_t").EnableCalculation = False
  
  Sheets("Calc_H_m").EnableCalculation = False
  Sheets("Calc_C_m").EnableCalculation = False
  
  Sheets("Graphs").EnableCalculation = False
  Sheets("Clipboard").EnableCalculation = False
  
  Sheets("Output_m").EnableCalculation = False
  Sheets("Output_t").EnableCalculation = False
  
  Sheets("ClimDat_m").EnableCalculation = False
  Sheets("ClimData").EnableCalculation = False
  
  Sheets("VCS_special").EnableCalculation = False
  '
  'So apparently (see previous comment) the following code is not equivalent to the list above
  
  Application.Calculation = xlCalculationManual
  Application.ScreenUpdating = False

  'Redimensioning arrays: FUTURE)
  ' ReDim X_xxx(1 To Eln) As Variant
  '
  ' START TIME LOOP - START TIME LOOP - START TIME LOOP - START TIME LOOP - START TIME LOOP -
  '
 
  For Tstepi = Tstep1 To Tstepn        'Start time loop
  '
  ' INITIALIZATION PER TIME STEP:
  '
  '
    ' First: get the relevant time-dependent input data:
    '
    '
    Rstep = Tstepi + RTstep00       'Row number at sheet for input and output data per timestep
    ' Note: from Input_t several columns are used at this moment; others are transfered to Output_t and read from there.
    ' To be CHECKed at later date: this requires recalc of Output_t per timestep;
    ' does recalc of Output_t take much calc. time?
    'CHECK DYNAMIC! Recalc shall not be suppressed if we want possibly variable data from the sheet
    'DYNAMIC VCS:
    'CHECK: could be refined: read VCS = 0 or 1 => only if =1.
    ' On the other hand: this calculation may be needed for any type of dynamic element!
    'Update Output_t previous timestep:
    Sheets("Output_t").Range(Cells(Rstep - 1, 1).Address, Cells(Rstep - 1, 52).Address).Calculate
    'For row before first row, row1-1, there is nothing to calculate:
    ' but no problem, because it is not row=0
    Sheets("VCS_special").Range(Cells(Rstep - 1, 1).Address, Cells(Rstep, 8).Address).Calculate
    Sheets("Input_t").Range(Cells(Rstep, 1).Address, Cells(Rstep, 30).Address).Calculate
    'Vent. flow is converted to Hve in Output_t, so also update sheet Output_t:
    Sheets("Output_t").Range(Cells(Rstep, 1).Address, Cells(Rstep, 52).Address).Calculate
   'Read input data:
    Input_t = Sheets("Input_t").Range(Cells(Rstep, 1).Address, Cells(Rstep, 24).Address)
    'Some dynamic input data have been adapted in Output_t, so:
    Output_t = Sheets("Output_t").Range(Cells(Rstep, 1).Address, Cells(Rstep, 27).Address) 'if too many: does not matter
    '
    Theta_H_set = Input_t(1, 13)
    Theta_C_set = Input_t(1, 14)
    '
    '
    Dtime = Input_t(1, 2)
    Month = Input_t(1, 3)
    If Week < Input_t(1, 4) Then
      ' message?
    End If
    Week = Input_t(1, 4)
   
    '
    ' Second: determine time variable conditions:
    ' Store temperature vector previous timestep
    ' Previous time step (with or without iteration) always ends with solution for colB=1,
    ' except if maxH or maxC used -> colB_act!
    For ri = 1 To Rn
      Theta_old(ri) = VecB(ri, colB_act)
    Next ri                  ' Ri = 1 To Rn
    '
    If Theta_H_set < -995 Then
      Phi_H_nd_max_act = 0
    Else
      Phi_H_nd_max_act = Phi_H_nd_max
    End If
    If Theta_C_set > 995 Then
      Phi_C_nd_max_act = 0
    Else
      Phi_C_nd_max_act = Phi_C_nd_max
    End If
    ' Which mode(s) to calculate:
    ' First pass of this time step: calculate three modes:
    ' As basis for interpolation (called "step 2" in the standard?)
    ' nrHCmodes and mod_HC are used
    '   cBi is the column number of the vector: (1, 2 or 3)
    ' Normally 3 modes:
    '   1: no H, no C
    '   2: H colB_H=2
    '   3: C colB_C=3
    '
    ' But if Phi_H_nd_max_act = 0 or Phi_C_nd_max_act = 0:
    ' the number of modes can be reduced:
    ' If Phi_C_nd_max_act = 0, two modes, namely:
    '   1: no H, no C
    '   2: H -> colB_H=2
    ' If Phi_H_nd_max_act = 0, two modes, namely:
    '   1: no H, no C
    '   2: C -> colB_C=2
    ' If Phi_H_nd_max_act = 0 and Phi_C_nd_max_act = 0, single mode, namely:
    '   1: no H, no C
    '
    ' In case of recalculation with actual H or C power: single mode, namely:
    '   1: actual H or C
    '
    '
    Phi_HC_nd_calc(1) = 0
    If Phi_H_nd_max_act = 0 And Phi_C_nd_max_act = 0 Then
    ' No H and no C
      nrHCmodes = 1
    ElseIf Phi_C_nd_max_act = 0 Then
    ' Only H
      colB_H = 2
      nrHCmodes = 2
      Phi_HC_nd_calc(2) = Phi_H_nd_max_act
    ElseIf Phi_H_nd_max_act = 0 Then
    ' Only C
      colB_C = 2
      nrHCmodes = 2
      Phi_HC_nd_calc(2) = Phi_C_nd_max_act ' Note: C power already put in column 2,
    Else
    ' Both H and C
      nrHCmodes = 3
      colB_H = 2
      colB_C = 3
      Phi_HC_nd_calc(2) = Phi_H_nd_max_act
      Phi_HC_nd_calc(3) = Phi_C_nd_max_act
    End If
 
    '
'   HERE STARTS RECALCULATION IN CASE OF ITERATION
'   Note: if specific input data changed as part of the iteration, these have to be re-read here:
    ' Reset all matrix and vector values to zero before (re-)calculation
10    For ri = 1 To Rn
      For cBi = 1 To 3 'Safest to reset all three ...
       VecB(ri, cBi) = 0#
      Next cBi
      For ci = 1 To Cn
        MatA(ri, ci) = 0#
      Next ci                ' Ci = 1 To Cn
    Next ri                  ' Ri = 1 To Rn
    '
    'Read DYNAMIC properties ==SPECIAL!!:
     For Eli = 1 To Eln
       If TypeSub_eli(Eli) = "W" And TypeClass_eli(Eli) <> "F" Then
        If TypeClass_eli(Eli) = "V1" Then
          h_pli_eli(1, Eli) = Input_t(1, 21) '#Pl#: no change needed: only one layer for windows
         tau_sol_eli(1, Eli) = Input_t(1, 22) '#Pl#: no change needed: tau_sol is not per node
        End If
        If TypeClass_eli(Eli) = "V2" Then
          h_pli_eli(1, Eli) = Input_t(1, 23)
          tau_sol_eli(1, Eli) = Input_t(1, 24)
        End If
      End If
     Next Eli
      Phi_sol_zi = 0 'DYNAMIC: New Sept. 2014: now in macro and not directly on sheet, because tau_sol may vary
    For Eli = 1 To Eln
       '  Calculate direct solar gains:
       '  Phi_sol_zi = SUM over eli: (tau_sol_eli * A_eli * I_sol_eli)
        ' Only then? Yes: = OPaque or Window (not GR.floor or AD or ...)
        If Type_eli(Eli) = "EXT" Then
          'External element:
           Phi_sol_zi = Phi_sol_zi + tau_sol_eli(1, Eli) * A_eli(1, Eli) * Output_t(1, 7 + Eli)
        End If
    Next Eli
      'CHECK DYNAMIC: Do we need this as new column output in Output_t ?
      '
    ' FILL
    ' Start to (re-)fill matrix matA and vector VecB
    ' First row: zone air node balance
    ri = 1
    '
    'BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB():
    '  Theta_(int;air) node
    '  H_tb (is constant)
    '  Output_t(1, 6)  =  H_ve(t)
    '  Input_t(1, 19)   =  Theta_sup(t)
    '  Output_t(1, 7)  =  Phi_(int;zi)(t)
    '  OVERRULED BECAUSE OF DYNAMIC ELEMENTS Output_t(1, 28) =  Phi_(sol;zi)(t)
    '  Output_t(1, 30) =  Phi_(HC;nd;zi)(t) 'but this is an output
    '
    XTemp = H_tb * Sheets("ClimData").Cells(Rstep + RTstepClim, 13).Value _
    + Output_t(1, 6) * Input_t(1, 19) _
    + f_int_c * Output_t(1, 7) _
    + f_sol_c * Phi_sol_zi _
    + (C_int / Dtime) * Theta_old(ri)
    '
    For cBi = 1 To nrHCmodes
      If Phi_HC_nd_calc(cBi) > 0 Then
        f_HC_c = f_H_c
      Else
        f_HC_c = f_C_c
      End If
      VecB(ri, cBi) = VecB(ri, cBi) + XTemp + f_HC_c * Phi_HC_nd_calc(cBi)
    Next cBi
    '
    '
    'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA():
    ci = 1
    ' Theta_(int;air) node:
    '  H_tb (is constant = thermal bridges)
    '  Output_t(1, 6) =  H_ve(t)
    MatA(ri, ci) = MatA(ri, ci) + (C_int / Dtime) + Ah_ci + H_tb + Output_t(1, 6)
    
    ' Other columns:
    For Eli = 1 To Eln
      ' Still first row, next columns, if linked to indoor air node:
      ' #Pl# For each internal surface node Pli=Pln, of each element Eli
      ' Link to Theta_pli;eli) node:
      '#Pl# Pli = 1
      Pli = Pln(Eli)
      ci = 1 + PlnSum(Eli) + Pli
      MatA(ri, ci) = MatA(ri, ci) - A_eli(1, Eli) * h_ci_eli(1, Eli)
    Next Eli
    ' DONE: First row (balance equation for indoor air node)
    '
    ' Next rows: Balance equations for all other nodes: all nodes Pli for each element Eli
    For Eli = 1 To Eln

      For Pli = 1 To Pln(Eli)
        ri = ri + 1          'Row: Equation for node Pli, Eli
        ' BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB():
        ' each node:
        ' Output_t(1, 7 + Eli) = I_sol(or, tilt)(t)
        XTemp = a_sol_pli_eli(Pli, Eli) * Output_t(1, 7 + Eli) _
          + (kappa_pli_eli(Pli, Eli) / Dtime) * Theta_old(ri)
         For cBi = 1 To nrHCmodes
           VecB(ri, cBi) = VecB(ri, cBi) + XTemp
         Next cBi
          
        '#Pl# If Pli = 1 Then
        If Pli = Pln(Eli) Then
           ' #Pl# For each element: last node is indoor (zone) facing surface:
           ' balance in W/m2!
           '  Output_t(1, 7)  =  Phi_(int;zi)(t)
           '  OVERRULED BECAUSE OF DYNAMIC ELEMENTS: Output_t(1, 28) =  Phi_(sol;zi)(t)
           '  Output_t(1, 30) =  Phi_(HC;nd;zi)(t), but this is output
           '
          
           XTemp = ((1 - f_int_c) * Output_t(1, 7) + (1 - f_sol_c) * Phi_sol_zi)
           For cBi = 1 To nrHCmodes
             If Phi_HC_nd_calc(cBi) > 0 Then
               f_HC_c = f_H_c
             Else
               f_HC_c = f_C_c
             End If
             VecB(ri, cBi) = VecB(ri, cBi) _
                + (XTemp + (1 - f_HC_c) * Phi_HC_nd_calc(cBi)) / A_eli_tot
           Next cBi

        '#Pl# ElseIf Pli = Pln(Eli) Then
        ElseIf Pli = 1 Then
          ' #Pl# For each element: first node is external (environment) facing surface:
          
          If Type_eli(Eli) = "EXT" Then
          'External element:
            'For external elements: Sheets("ClimData").Cells(Rstep, 6).Value = Theta_(e;air)(t)
            'Here the sky radiation has been added too
            XTemp = (h_ce_eli(1, Eli) + h_re_eli(1, Eli)) * Sheets("ClimData").Cells(Rstep + RTstepClim, 13).Value _
              - F_sk_eli(1, Eli) * h_re_eli(1, Eli) * Delta_Theta_er
            For cBi = 1 To nrHCmodes
               VecB(ri, cBi) = VecB(ri, cBi) + XTemp
            Next cBi
          ElseIf Type_eli(Eli) = "AD" Then
          ' Adiabatic element:
          ' Nothing to do
          '
          ElseIf Type_eli(Eli) = "GR" Then
           ' Ground floor element:
           XTemp = (1 / R_gr_ve) * Theta_gr_ve(Month)
             For cBi = 1 To nrHCmodes
               VecB(ri, cBi) = VecB(ri, cBi) + XTemp
             Next cBi
          End If
          
        Else
         ' For each element: node inside element (if any):
         ' Nothing to add
        End If
        '
        '
        '
        ' AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA():
        ' For each row ri > 1 (each Eli, Pli):
        '
        'Diagonal node:
        ci = 1 + PlnSum(Eli) + Pli
        MatA(ri, ci) = MatA(ri, ci) + kappa_pli_eli(Pli, Eli) / Dtime
        '
        '#Pl# If Pli = 1 Then
        If Pli = Pln(Eli) Then 'The diagonal node is an indoor facing surface
          ' for the diagonal node ri, ci)
          MatA(ri, ci) = MatA(ri, ci) + h_ci_eli(1, Eli) + h_ri_eli_mn
          ' for other nodes:
          MatA(ri, 1) = MatA(ri, 1) - h_ci_eli(1, Eli) 'Link to indoor air node (ci=1)
          For Elk = 1 To Eln
            ' Long wave (thermal) radiative link to other indoor (zone) facing surfaces
            '#Pl# Plk = 1
            Plk = Pln(Elk)
            ck = 1 + PlnSum(Elk) + Plk
            MatA(ri, ck) = MatA(ri, ck) _
            - (A_eli(1, Elk) / A_eli_tot) * h_ri_eli(1, Elk)
          Next Elk
          
        '#Pl# ElseIf Pli = Pln(Eli) Then
        ElseIf Pli = 1 Then 'The diagonal node is an external surface
          ' For the diagonal node:
          If Type_eli(Eli) = "EXT" Then
            ' External element:
            MatA(ri, ci) = MatA(ri, ci) + h_ce_eli(1, Eli) + h_re_eli(1, Eli)
          ElseIf Type_eli(Eli) = "AD" Then
            ' Adiabatic element:
            ' Nothing to do
            '
          ElseIf Type_eli(Eli) = "GR" Then
            ' Ground floor element:
            MatA(ri, ci) = MatA(ri, ci) + 1 / R_gr_ve
          End If
        Else
         ' For node inside element (if any):
         ' Nothing to add here, but see below
        End If
        
        If Pli > 1 Then '#Pl#: no need to change: only needed to avoid nodes less than 1 or higher than Pln
          MatA(ri, ci) = MatA(ri, ci) + h_pli_eli(Pli - 1, Eli) ' diagonal node ri, ci
          MatA(ri, ci - 1) = MatA(ri, ci - 1) - h_pli_eli(Pli - 1, Eli) ' lower node
        End If
        If (Pli < Pln(Eli)) Then
          MatA(ri, ci) = MatA(ri, ci) + h_pli_eli(Pli, Eli) ' diagonal node ri, ci
          MatA(ri, ci + 1) = MatA(ri, ci + 1) - h_pli_eli(Pli, Eli) 'higher node
        End If
         
      Next Pli               ' Pli = 1 To Pln(Eli)
      
    Next Eli                 ' Eli = 1 To Eln
    '
    ' End matrix matA and vector vecB filled for this time step
    '
    '
    '#Pl# Produce output VecB and MatA:
    If Tstepi = Tstepn Then
       'Output for Vector B:
       ri = 1
       ci = 1
       Sheets("Test_MatVec").Cells(1, 1).Value = "vecB"
       Sheets("Test_MatVec").Cells(1, ci + 2).Value = "air"
       Sheets("Test_MatVec").Cells(2, 1).Value = "[Output for last timestep= " & Tstepn & "]"
       Sheets("Test_MatVec").Cells(3, ci + 2).Value = VecB(ri, 1) ' "1)" is for the HC mode=1
      
       For Eli = 1 To Eln
        For Pli = 1 To Pln(Eli)
            ri = ri + 1
            ci = ci + 1
            Sheets("Test_MatVec").Cells(1, ci + 2).Value = Eli
            Sheets("Test_MatVec").Cells(2, ci + 2).Value = Pli
            Sheets("Test_MatVec").Cells(3, ci + 2).Value = VecB(ri, 1) ' "1)" is for the HC mode=1
         Next Pli
       Next Eli
       
       'Output for Matrix A:
       Sheets("Test_MatVec").Cells(4, 1).Value = "matA"
       Sheets("Test_MatVec").Cells(5, 1).Value = "air"
       ri = 1
       ci = 1
       Sheets("Test_MatVec").Cells(ri + 4, ci + 2).Value = MatA(ri, ci)
       For Elk = 1 To Eln
          For Plk = 1 To Pln(Elk)
            ci = ci + 1
            Sheets("Test_MatVec").Cells(ri + 4, ci + 2).Value = MatA(ri, ci)
          Next Plk
         Next Elk
       For Eli = 1 To Eln
        For Pli = 1 To Pln(Eli)
          ri = ri + 1
          ci = 1
          Sheets("Test_MatVec").Cells(ri + 4, ci + 2).Value = MatA(ri, ci)
          Sheets("Test_MatVec").Cells(ri + 4, 1).Value = Eli
          Sheets("Test_MatVec").Cells(ri + 4, 2).Value = Pli
          For Elk = 1 To Eln
          For Plk = 1 To Pln(Elk)
            ci = ci + 1
            Sheets("Test_MatVec").Cells(ri + 4, ci + 2).Value = MatA(ri, ci)
          Next Plk
         Next Elk
        Next Pli
       Next Eli
    End If
     '#Pl# END TEST + END TEST + END TEST + END TEST + END TEST + END TEST + END TEST
    '
    ' Matrix solving is  done in GaussJordan
    '
    '
  
    ' Matrix solved with function GaussJordan
    XTest = MatA(1, 1)
    Call GaussJordan
    '
    ' Produce output:
    ' Communicating with the sheets here is unavoidable, unless the data are stored for the whole time series
    ' CHECK: Work around?: writing to csv file is quicker? Then at the end: read from CSV and write to sheets
    '
    ' VecB returned by GaussJ is the Theta vector
    '
    '
    '#Pl# Produce output VecB (=Theta):
    If Tstepi = Tstepn Then
       'Output for Vector B:
       ri = 1
       ci = 1
       Sheets("Test_MatVec").Cells(1 + Cn + 6, 1).Value = "Theta"
       Sheets("Test_MatVec").Cells(1 + Cn + 6, ci + 2).Value = "air"
       Sheets("Test_MatVec").Cells(3 + Cn + 6, ci + 2).Value = VecB(ri, 1) ' "1)" is for the HC mode=1
      
       For Eli = 1 To Eln
        For Pli = 1 To Pln(Eli)
            ri = ri + 1
            ci = ci + 1
            Sheets("Test_MatVec").Cells(1 + Cn + 6, ci + 2).Value = Eli
            Sheets("Test_MatVec").Cells(2 + Cn + 6, ci + 2).Value = Pli
            Sheets("Test_MatVec").Cells(3 + Cn + 6, ci + 2).Value = VecB(ri, 1) ' "1)" is for the HC mode=1
         Next Pli
       Next Eli
    End If
    '
    '
    '
    ' Calculate operative temperature:
    For cBi = 1 To nrHCmodes

        Theta_int_air(cBi) = VecB(1, cBi)
        Theta_int_r_mn(cBi) = 0#
        For Eli = 1 To Eln
           '#Pl# ri = 1 + PlnSum(Eli) + 1 'for each Pli=1
           ri = 1 + PlnSum(Eli) + Pln(Eli) 'for each Pli=Pln
           Theta_int_r_mn(cBi) = Theta_int_r_mn(cBi) + A_eli(1, Eli) * VecB(ri, cBi)
        Next Eli
        Theta_int_r_mn(cBi) = Theta_int_r_mn(cBi) / A_eli_tot
        Theta_int_op(cBi) = 0.5 * Theta_int_air(cBi) + 0.5 * Theta_int_r_mn(cBi)
    Next cBi
    '
    If nrHCmodes > 1 Then
      '
      If Theta_int_op(1) < Theta_H_set Then
          Theta_op_set = Theta_H_set
          ' If there is H and nrcColsB > 1, then H is always the second mode, but to be sure: use colB_H:
          Phi_HC_nd_act = Phi_H_nd_max _
           * (Theta_op_set - Theta_int_op(1)) / (Theta_int_op(colB_H) - Theta_int_op(1))
          '
          If Phi_HC_nd_act > Phi_H_nd_max Then
          ' DONE!!: max power used; output for this time step is the output of mode 2
           Phi_HC_nd_act = Phi_H_nd_max
           Theta_op_act = Theta_int_op(colB_H)
           colB_act = colB_H ' This means temperature vector as output for this time step is colB=2
          Else
           Phi_HC_nd_calc(1) = Phi_HC_nd_act
           Theta_op_act = Theta_op_set
           colB_act = 1 ' But for recalculation not used
           nrHCmodes = 1
           GoTo 10 ' Recalculate
          End If
      ElseIf Theta_int_op(1) > Theta_C_set Then
          Theta_op_set = Theta_C_set
          ' If there is C and nrcColsB > 1, then C is the second or the third colB:
          Phi_HC_nd_act = Phi_C_nd_max _
               * (Theta_op_set - Theta_int_op(1)) / (Theta_int_op(colB_C) - Theta_int_op(1))
          '
          If Phi_HC_nd_act < Phi_C_nd_max Then 'ATTENTION: "<" because negative values!
          ' DONE!!: max power used; output for this time step is the output of mode 2
           Phi_HC_nd_act = Phi_C_nd_max
           Theta_op_act = Theta_int_op(colB_C)
           colB_act = colB_C ' This means temperature vector as output for this time step is colB=colB_C
          Else
           Phi_HC_nd_calc(1) = Phi_HC_nd_act
           Theta_op_act = Theta_op_set
           colB_act = 1 ' But for recalculation not used
           nrHCmodes = 1
           GoTo 10
          End If
       Else
        'if neither: No need for H power nor for C power
        ' Occurs only in case of nrHCmodes > 1, if
        ' at first pass (nrHCmodes=3) temperature between H and C setpoint
        Phi_HC_nd_act = 0
        Theta_op_act = Theta_int_op(1)
        colB_act = 1 ' This means temperature vector as output for this time step is colB=1
      End If
    Else   ' ("If nrHCmodes > 1")
      ' Means: nrHCmodes = 1, so either no H and C at all, or already recalculated H or C
      ' DONE for this time step!
        Phi_HC_nd_act = Phi_HC_nd_calc(1)
        Theta_op_act = Theta_int_op(1)
        colB_act = 1
    End If ' ("If nrHCmodes > 1")
    '
    'Postprocessing this timestep:
       
    ' Write that this hour has been calculated
    Sheets("Output_t").Cells(Rstep, 29).Value = 1
    ' Write actual HC power to output sheet:
    Sheets("Output_t").Cells(Rstep, 30).Value = Phi_HC_nd_act
    '
    ' Write actual op. temperature to output sheet:
    Sheets("Output_t").Cells(Rstep, 31).Value = Theta_op_act
    Sheets("Output_t").Cells(Rstep, 32).Value = Theta_int_air(colB_act)
    Sheets("Output_t").Cells(Rstep, 33).Value = Theta_int_r_mn(colB_act)
   'DYNAMIC:
    Sheets("Output_t").Range(Cells(Rstep, 29).Address, Cells(Rstep, 33).Address).Calculate
   
   '
   '
   'CHECK: is this a good moment to recalculate the sheets?
   ' No! Unless put back to manual at start of next timestep...
   ' If DYNAMIC: some recalc may be needed
 ' DYNAMIC:
 'recalc placed at beginning of timestep, bevause there it picks up the values of the previous time step

  ' Tstepi = Tstep1 To Tstepn
  ' END THIS TIMESTEP - END THIS TIMESTEP - END THIS TIMESTEP
  Next Tstepi
  ' END TIME LOOP - END TIME LOOP - END TIME LOOP - END TIME LOOP
  '
  ' Now is the moment to recalculate
  Sheets("INFO").EnableCalculation = True
  Sheets("Explanation").EnableCalculation = True
  
  Sheets("Method_input").EnableCalculation = True
  Sheets("Input_GrFl").EnableCalculation = True
  Sheets("Calc_GrFl").EnableCalculation = True
  
  Sheets("Method_calculation").EnableCalculation = True
  Sheets("Method_output").EnableCalculation = True
  
  Sheets("Input_0").EnableCalculation = True
  Sheets("Input_p").EnableCalculation = True
  Sheets("Input_t").EnableCalculation = True
  
  Sheets("Calc_H_m").EnableCalculation = True
  Sheets("Calc_C_m").EnableCalculation = True
  
  Sheets("Graphs").EnableCalculation = True
  Sheets("Clipboard").EnableCalculation = True
  
  Sheets("Output_m").EnableCalculation = True
  Sheets("Output_t").EnableCalculation = True
  
  Sheets("ClimDat_m").EnableCalculation = True
  Sheets("ClimData").EnableCalculation = True
  
  Sheets("VCS_special").EnableCalculation = True
 
  '
  '
  'Record calculation ending date and time
  RunDate1 = Date
  RunTime1 = Time
  Worksheets("Output_m").Range("b2").Value = "Calculation ended at: " & RunDate1 & " " & RunTime1
  
'
'
  Sheets("Output_m").Select         ' To show this sheet after calculation finished
  'Cells(1, 1).Select
  Range("A1").Select
  
  Application.Calculation = xlCalculationAutomatic   ' CHECK: is this a good moment to recalculate the sheets?
  Application.ScreenUpdating = True ' To update the screen of the computer
  '
  '
End Sub            ' FillMatrix()