-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRadiation_Functions.py
2253 lines (1851 loc) · 100 KB
/
Radiation_Functions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import os
import numpy
import math
import copy
'''
Water Functions:
Developed by Mohsen Moradi and Amir A. Aliabadi
Atmospheric Innovations Research (AIR) Laboratory, University of Guelph, Guelph, Canada
Last update: May 2021
Originally developed by Naika Meili
'''
class RadiationFunctions(object):
def TotalLWRabsorbed(self,TemperatureC,geometry,MeteoData,FractionsGround,PropOpticalGround,PropOpticalWall,
PropOpticalTree,ParTree,ViewFactor):
"""
------
INPUT:
TemperatureC: Temperature of the exterior surfaces [K]
geometry: Normalized geometric parameters
MeteoData: Forcing variables
FractionsGround:
PropOpticalGround: Optical properties of the ground (albedo and emissivity)
PropOpticalWall: Optical properties of the wall (albedo and emissivity)
PropOpticalTree: Optical properties of the trees (albedo and emissivity)
ParTree: Trees parameter
ViewFactor:
-------
OUTPUT:
LWRin_t: Incoming longwave radiations [W m^-2]
LWRout_t: Outgoing longwave radiations [W m^-2]
LWRabs_t: Absorbed longwave radiations [W m^-2]
LWREB_t: Energy conservation of longwave radiations [W m^-2]
"""
Tgrimp = TemperatureC[0]
Tgbare = TemperatureC[1]
Tgveg = TemperatureC[2]
Twsun = TemperatureC[3]
Twshade = TemperatureC[4]
Ttree = TemperatureC[5]
if ParTree.trees == 1:
# Longwave radiation absorbed without trees
LWRin_nT, LWRout_nT, LWRabs_nT, LWREB_nT = self.LWRabsorbedNoTree(geometry.hcanyon,geometry.wcanyon,MeteoData.LWR,
FractionsGround.fveg,FractionsGround.fbare,FractionsGround.fimp,
PropOpticalWall.emissivity,PropOpticalGround.eveg,
PropOpticalGround.ebare,PropOpticalGround.eimp,Tgrimp,Tgbare,
Tgveg,Twsun,Twshade,ViewFactor)
# Longwave radiation absorbed with trees
LWRin_T, LWRout_T, LWRabs_T, LWREB_T = self.LWRabsorbedWithTrees(geometry.hcanyon, geometry.wcanyon, geometry.radius_tree,
MeteoData.LWR, FractionsGround.fveg, FractionsGround.fbare,
FractionsGround.fimp, PropOpticalWall.emissivity,
PropOpticalTree.emissivity, PropOpticalGround.eveg,
PropOpticalGround.ebare, PropOpticalGround.eimp, Tgrimp,
Tgbare, Tgveg, Twsun, Twshade, Ttree, ViewFactor)
class LWRin_t_Def():
pass
LWRin_t = LWRin_t_Def()
LWRin_t.LWRinGroundImp = ParTree.ftree*LWRin_T.LWRinGroundImp + (1-ParTree.ftree)*LWRin_nT.LWRinGroundImp
LWRin_t.LWRinGroundBare = ParTree.ftree*LWRin_T.LWRinGroundBare + (1-ParTree.ftree)*LWRin_nT.LWRinGroundBare
LWRin_t.LWRinGroundVeg = ParTree.ftree*LWRin_T.LWRinGroundVeg + (1-ParTree.ftree)*LWRin_nT.LWRinGroundVeg
LWRin_t.LWRinTree = ParTree.ftree*LWRin_T.LWRinTree + (1-ParTree.ftree)*LWRin_nT.LWRinTree
LWRin_t.LWRinWallSun = ParTree.ftree*LWRin_T.LWRinWallSun + (1-ParTree.ftree)*LWRin_nT.LWRinWallSun
LWRin_t.LWRinWallShade = ParTree.ftree*LWRin_T.LWRinWallShade + (1-ParTree.ftree)*LWRin_nT.LWRinWallShade
LWRin_t.LWRinTotalGround = ParTree.ftree*LWRin_T.LWRinTotalGround + (1-ParTree.ftree)*LWRin_nT.LWRinTotalGround
LWRin_t.LWRinTotalCanyon = ParTree.ftree*LWRin_T.LWRinTotalCanyon + (1-ParTree.ftree)*LWRin_nT.LWRinTotalCanyon
class LWRout_t_Def():
pass
LWRout_t = LWRout_t_Def()
LWRout_t.LWRoutGroundImp = ParTree.ftree*LWRout_T.LWRoutGroundImp + (1-ParTree.ftree)*LWRout_nT.LWRoutGroundImp
LWRout_t.LWRoutGroundBare = ParTree.ftree*LWRout_T.LWRoutGroundBare + (1-ParTree.ftree)*LWRout_nT.LWRoutGroundBare
LWRout_t.LWRoutGroundVeg = ParTree.ftree*LWRout_T.LWRoutGroundVeg + (1-ParTree.ftree)*LWRout_nT.LWRoutGroundVeg
LWRout_t.LWRoutTree = ParTree.ftree*LWRout_T.LWRoutTree + (1-ParTree.ftree)*LWRout_nT.LWRoutTree
LWRout_t.LWRoutWallSun = ParTree.ftree*LWRout_T.LWRoutWallSun + (1-ParTree.ftree)*LWRout_nT.LWRoutWallSun
LWRout_t.LWRoutWallShade = ParTree.ftree*LWRout_T.LWRoutWallShade + (1-ParTree.ftree)*LWRout_nT.LWRoutWallShade
LWRout_t.LWRoutTotalGround = ParTree.ftree*LWRout_T.LWRoutTotalGround + (1-ParTree.ftree)*LWRout_nT.LWRoutTotalGround
LWRout_t.LWRoutTotalCanyon = ParTree.ftree*LWRout_T.LWRoutTotalCanyon + (1-ParTree.ftree)*LWRout_nT.LWRoutTotalCanyon
class LWRabs_t_Def():
pass
LWRabs_t = LWRabs_t_Def()
LWRabs_t.LWRabsGroundImp = ParTree.ftree*LWRabs_T.LWRabsGroundImp + (1-ParTree.ftree)*LWRabs_nT.LWRabsGroundImp
LWRabs_t.LWRabsGroundBare = ParTree.ftree*LWRabs_T.LWRabsGroundBare + (1-ParTree.ftree)*LWRabs_nT.LWRabsGroundBare
LWRabs_t.LWRabsGroundVeg = ParTree.ftree*LWRabs_T.LWRabsGroundVeg + (1-ParTree.ftree)*LWRabs_nT.LWRabsGroundVeg
LWRabs_t.LWRabsTree = ParTree.ftree*LWRabs_T.LWRabsTree + (1-ParTree.ftree)*LWRabs_nT.LWRabsTree
LWRabs_t.LWRabsWallSun = ParTree.ftree*LWRabs_T.LWRabsWallSun + (1-ParTree.ftree)*LWRabs_nT.LWRabsWallSun
LWRabs_t.LWRabsWallShade = ParTree.ftree*LWRabs_T.LWRabsWallShade + (1-ParTree.ftree)*LWRabs_nT.LWRabsWallShade
LWRabs_t.LWRabsTotalGround = ParTree.ftree*LWRabs_T.LWRabsTotalGround + (1-ParTree.ftree)*LWRabs_nT.LWRabsTotalGround
LWRabs_t.LWRabsTotalCanyon = ParTree.ftree*LWRabs_T.LWRabsTotalCanyon + (1-ParTree.ftree)*LWRabs_nT.LWRabsTotalCanyon
class LWREB_t_Def():
pass
LWREB_t = LWREB_t_Def()
LWREB_t.LWREBGroundImp = ParTree.ftree*LWREB_T.LWREBGroundImp + (1-ParTree.ftree)*LWREB_nT.LWREBGroundImp
LWREB_t.LWREBGroundBare = ParTree.ftree*LWREB_T.LWREBGroundBare + (1-ParTree.ftree)*LWREB_nT.LWREBGroundBare
LWREB_t.LWREBGroundVeg = ParTree.ftree*LWREB_T.LWREBGroundVeg + (1-ParTree.ftree)*LWREB_nT.LWREBGroundVeg
LWREB_t.LWREBTree = ParTree.ftree*LWREB_T.LWREBTree + (1-ParTree.ftree)*LWREB_nT.LWREBTree
LWREB_t.LWREBWallSun = ParTree.ftree*LWREB_T.LWREBWallSun + (1-ParTree.ftree)*LWREB_nT.LWREBWallSun
LWREB_t.LWREBWallShade = ParTree.ftree*LWREB_T.LWREBWallShade + (1-ParTree.ftree)*LWREB_nT.LWREBWallShade
LWREB_t.LWREBTotalGround = ParTree.ftree*LWREB_T.LWREBTotalGround + (1-ParTree.ftree)*LWREB_nT.LWREBTotalGround
LWREB_t.LWREBTotalCanyon = ParTree.ftree*LWREB_T.LWREBTotalCanyon + (1-ParTree.ftree)*LWREB_nT.LWREBTotalCanyon
# The absorbed radiation by the tree is not averaged as it is per tree surface
LWRin_t.LWRinTree = LWRin_T.LWRinTree
LWRout_t.LWRoutTree = LWRout_T.LWRoutTree
LWRabs_t.LWRabsTree = LWRabs_T.LWRabsTree
LWREB_t.LWREBTree = LWREB_T.LWREBTree
elif ParTree.trees == 0:
# Longwave radiation absorbed without trees
LWRin_nT, LWRout_nT, LWRabs_nT, LWREB_nT = self.LWRabsorbedNoTree(geometry.hcanyon,geometry.wcanyon,MeteoData.LWR,
FractionsGround.fveg,FractionsGround.fbare,FractionsGround.fimp,
PropOpticalWall.emissivity,PropOpticalGround.eveg,
PropOpticalGround.ebare,PropOpticalGround.eimp,
Tgrimp,Tgbare,Tgveg,Twsun,Twshade,ViewFactor)
LWRin_t = LWRin_nT
LWRout_t = LWRout_nT
LWRabs_t = LWRabs_nT
LWREB_t = LWREB_nT
return LWRin_t,LWRout_t,LWRabs_t,LWREB_t
def TotalSWRabsorbed(self,geometry,FractionsGround,ParTree,PropOpticalGround,PropOpticalWall,PropOpticalTree,
ParVegTree,MeteoData,SunPosition,ViewFactor):
"""
------
INPUT:
geometry: Normalized geometric parameters
FractionsGround:
ParTree: Trees parameter
PropOpticalGround: Optical properties of the ground (albedo and emissivity)
PropOpticalWall: Optical properties of the wall (albedo and emissivity)
PropOpticalTree: Optical properties of the trees (albedo and emissivity)
ParVegTree: Trees parameter
MeteoData: Forcing variables
SunPosition: Sun angles
ViewFactor:
-------
OUTPUT:
SWRin_t: Incoming shortwave radiations [W m^-2]
SWRout_t: Outgoing shortwave radiations [W m^-2]
SWRabs_t: Absorbed shortwave radiations [W m^-2]
SWRabsDir_t: Absorbed direct shortwave radiations [W m^-2]
SWRabsDiff_t: Absorbed diffuse shortwave radiations [W m^-2]
SWREB_t: Energy conservation of shortwave radiations [W m^-2]
"""
if ParTree.trees == 1:
# Shortwave radiation absorbed without trees
SWRin_nT, SWRout_nT, SWRabs_nT, SWRabsDir_nT, SWRabsDiff_nT, SWREB_nT = \
self.SWRabsorbedNoTrees(geometry.hcanyon,geometry.wcanyon,FractionsGround.fveg,FractionsGround.fbare,FractionsGround.fimp,
PropOpticalWall.albedo,PropOpticalGround.aveg,PropOpticalGround.abare,PropOpticalGround.aimp,
MeteoData.SW_dir,MeteoData.SW_diff,SunPosition.theta_Z,SunPosition.theta_n,ViewFactor,ParVegTree)
# Shortwave radiation absorbed with trees
SWRin_T, SWRout_T, SWRabs_T, SWRabsDir_T, SWRabsDiff_T, SWREB_T = \
self.SWRabsorbedWithTrees(geometry.hcanyon,geometry.wcanyon,geometry.htree,geometry.radius_tree,geometry.distance_tree,
FractionsGround.fveg,FractionsGround.fbare,FractionsGround.fimp,PropOpticalWall.albedo,
PropOpticalGround.aveg,PropOpticalGround.abare,PropOpticalGround.aimp,PropOpticalTree.albedo,
ParVegTree.LAI,MeteoData.SW_dir,MeteoData.SW_diff,SunPosition.theta_Z,SunPosition.theta_n,
ViewFactor,ParVegTree)
class SWRin_t_Def():
pass
SWRin_t = SWRin_t_Def()
SWRin_t.SWRinGroundImp = ParTree.ftree * SWRin_T.SWRinGroundImp + (1 - ParTree.ftree) * SWRin_nT.SWRinGroundImp
SWRin_t.SWRinGroundBare = ParTree.ftree * SWRin_T.SWRinGroundBare + (1 - ParTree.ftree) * SWRin_nT.SWRinGroundBare
SWRin_t.SWRinGroundVeg = ParTree.ftree * SWRin_T.SWRinGroundVeg + (1 - ParTree.ftree) * SWRin_nT.SWRinGroundVeg
SWRin_t.SWRinTree = ParTree.ftree * SWRin_T.SWRinTree + (1 - ParTree.ftree) * SWRin_nT.SWRinTree
SWRin_t.SWRinWallSun = ParTree.ftree * SWRin_T.SWRinWallSun + (1 - ParTree.ftree) * SWRin_nT.SWRinWallSun
SWRin_t.SWRinWallShade = ParTree.ftree * SWRin_T.SWRinWallShade + (1 - ParTree.ftree) * SWRin_nT.SWRinWallShade
SWRin_t.SWRinTotalGround = ParTree.ftree * SWRin_T.SWRinTotalGround + (1 - ParTree.ftree) * SWRin_nT.SWRinTotalGround
SWRin_t.SWRinTotalCanyon = ParTree.ftree * SWRin_T.SWRinTotalCanyon + (1 - ParTree.ftree) * SWRin_nT.SWRinTotalCanyon
class SWRout_t_Def():
pass
SWRout_t = SWRout_t_Def()
SWRout_t.SWRoutGroundImp = ParTree.ftree * SWRout_T.SWRoutGroundImp + (1 - ParTree.ftree) * SWRout_nT.SWRoutGroundImp
SWRout_t.SWRoutGroundBare = ParTree.ftree * SWRout_T.SWRoutGroundBare + (1 - ParTree.ftree) * SWRout_nT.SWRoutGroundBare
SWRout_t.SWRoutGroundVeg = ParTree.ftree * SWRout_T.SWRoutGroundVeg + (1 - ParTree.ftree) * SWRout_nT.SWRoutGroundVeg
SWRout_t.SWRoutTree = ParTree.ftree * SWRout_T.SWRoutTree + (1 - ParTree.ftree) * SWRout_nT.SWRoutTree
SWRout_t.SWRoutWallSun = ParTree.ftree * SWRout_T.SWRoutWallSun + (1 - ParTree.ftree) * SWRout_nT.SWRoutWallSun
SWRout_t.SWRoutWallShade = ParTree.ftree * SWRout_T.SWRoutWallShade + (1 - ParTree.ftree) * SWRout_nT.SWRoutWallShade
SWRout_t.SWRoutTotalGround = ParTree.ftree * SWRout_T.SWRoutTotalGround + (1 - ParTree.ftree) * SWRout_nT.SWRoutTotalGround
SWRout_t.SWRoutTotalCanyon = ParTree.ftree * SWRout_T.SWRoutTotalCanyon + (1 - ParTree.ftree) * SWRout_nT.SWRoutTotalCanyon
class SWRabs_t_Def():
pass
SWRabs_t = SWRabs_t_Def()
SWRabs_t.SWRabsGroundImp = ParTree.ftree * SWRabs_T.SWRabsGroundImp + (1 - ParTree.ftree) * SWRabs_nT.SWRabsGroundImp
SWRabs_t.SWRabsGroundBare = ParTree.ftree * SWRabs_T.SWRabsGroundBare + (1 - ParTree.ftree) * SWRabs_nT.SWRabsGroundBare
SWRabs_t.SWRabsGroundVeg = ParTree.ftree * SWRabs_T.SWRabsGroundVeg + (1 - ParTree.ftree) * SWRabs_nT.SWRabsGroundVeg
SWRabs_t.SWRabsTree = ParTree.ftree * SWRabs_T.SWRabsTree + (1 - ParTree.ftree) * SWRabs_nT.SWRabsTree
SWRabs_t.SWRabsWallSun = ParTree.ftree * SWRabs_T.SWRabsWallSun + (1 - ParTree.ftree) * SWRabs_nT.SWRabsWallSun
SWRabs_t.SWRabsWallShade = ParTree.ftree * SWRabs_T.SWRabsWallShade + (1 - ParTree.ftree) * SWRabs_nT.SWRabsWallShade
SWRabs_t.SWRabsTotalGround = ParTree.ftree * SWRabs_T.SWRabsTotalGround + (1 - ParTree.ftree) * SWRabs_nT.SWRabsTotalGround
SWRabs_t.SWRabsTotalCanyon = ParTree.ftree * SWRabs_T.SWRabsTotalCanyon + (1 - ParTree.ftree) * SWRabs_nT.SWRabsTotalCanyon
class SWRabsDir_t_Def():
pass
SWRabsDir_t = SWRabsDir_t_Def()
SWRabsDir_t.SWRabsGroundImp = ParTree.ftree * SWRabsDir_T.SWRabsGroundImp + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsGroundImp
SWRabsDir_t.SWRabsGroundBare = ParTree.ftree * SWRabsDir_T.SWRabsGroundBare + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsGroundBare
SWRabsDir_t.SWRabsGroundVeg = ParTree.ftree * SWRabsDir_T.SWRabsGroundVeg + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsGroundVeg
SWRabsDir_t.SWRabsTree = ParTree.ftree * SWRabsDir_T.SWRabsTree + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsTree
SWRabsDir_t.SWRabsWallSun = ParTree.ftree * SWRabsDir_T.SWRabsWallSun + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsWallSun
SWRabsDir_t.SWRabsWallShade = ParTree.ftree * SWRabsDir_T.SWRabsWallShade + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsWallShade
SWRabsDir_t.SWRabsTotalGround = ParTree.ftree * SWRabsDir_T.SWRabsTotalGround + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsTotalGround
SWRabsDir_t.SWRabsTotalCanyon = ParTree.ftree * SWRabs_T.SWRabsTotalCanyon + (1 - ParTree.ftree) * SWRabsDir_nT.SWRabsTotalCanyon
class SWRabsDiff_t_Def():
pass
SWRabsDiff_t = SWRabsDiff_t_Def()
SWRabsDiff_t.SWRabsGroundImp = ParTree.ftree * SWRabsDiff_T.SWRabsGroundImp + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsGroundImp
SWRabsDiff_t.SWRabsGroundBare = ParTree.ftree * SWRabsDiff_T.SWRabsGroundBare + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsGroundBare
SWRabsDiff_t.SWRabsGroundVeg = ParTree.ftree * SWRabsDiff_T.SWRabsGroundVeg + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsGroundVeg
SWRabsDiff_t.SWRabsTree = ParTree.ftree * SWRabsDiff_T.SWRabsTree + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsTree
SWRabsDiff_t.SWRabsWallSun = ParTree.ftree * SWRabsDiff_T.SWRabsWallSun + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsWallSun
SWRabsDiff_t.SWRabsWallShade = ParTree.ftree * SWRabsDiff_T.SWRabsWallShade + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsWallShade
SWRabsDiff_t.SWRabsTotalGround = ParTree.ftree * SWRabsDiff_T.SWRabsTotalGround + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsTotalGround
SWRabsDiff_t.SWRabsTotalCanyon = ParTree.ftree * SWRabs_T.SWRabsTotalCanyon + (1 - ParTree.ftree) * SWRabsDiff_nT.SWRabsTotalCanyon
class SWREB_t_Def():
pass
SWREB_t = SWREB_t_Def()
SWREB_t.SWREBGroundImp = ParTree.ftree * SWREB_T.SWREBGroundImp + (1 - ParTree.ftree) * SWREB_nT.SWREBGroundImp
SWREB_t.SWREBGroundBare = ParTree.ftree * SWREB_T.SWREBGroundBare + (1 - ParTree.ftree) * SWREB_nT.SWREBGroundBare
SWREB_t.SWREBGroundVeg = ParTree.ftree * SWREB_T.SWREBGroundVeg + (1 - ParTree.ftree) * SWREB_nT.SWREBGroundVeg
SWREB_t.SWREBTree = ParTree.ftree * SWREB_T.SWREBTree + (1 - ParTree.ftree) * SWREB_nT.SWREBTree
SWREB_t.SWREBWallSun = ParTree.ftree * SWREB_T.SWREBWallSun + (1 - ParTree.ftree) * SWREB_nT.SWREBWallSun
SWREB_t.SWREBWallShade = ParTree.ftree * SWREB_T.SWREBWallShade + (1 - ParTree.ftree) * SWREB_nT.SWREBWallShade
SWREB_t.SWREBTotalGround = ParTree.ftree * SWREB_T.SWREBTotalGround + (1 - ParTree.ftree) * SWREB_nT.SWREBTotalGround
SWREB_t.SWREBTotalCanyon = ParTree.ftree * SWREB_T.SWREBTotalCanyon + (1 - ParTree.ftree) * SWREB_nT.SWREBTotalCanyon
# The absorbed radiation by the tree is not averaged as it is per tree surface
SWRin_t.SWRinTree = SWRin_T.SWRinTree
SWRout_t.SWRoutTree = SWRout_T.SWRoutTree
SWRabs_t.SWRabsTree = SWRabs_T.SWRabsTree
SWRabsDir_t.SWRabsTree = SWRabsDir_T.SWRabsTree
SWRabsDiff_t.SWRabsTree = SWRabsDiff_T.SWRabsTree
SWREB_t.SWREBTree = SWREB_T.SWREBTree
elif ParTree.trees == 0:
SWRin_nT, SWRout_nT, SWRabs_nT, SWRabsDir_nT, SWRabsDiff_nT, SWREB_nT = \
self.SWRabsorbedNoTrees(geometry.hcanyon,geometry.wcanyon,FractionsGround.fveg,FractionsGround.fbare,FractionsGround.fimp,
PropOpticalWall.albedo,PropOpticalGround.aveg,PropOpticalGround.abare,PropOpticalGround.aimp,
MeteoData.SW_dir,MeteoData.SW_diff,SunPosition.theta_Z,SunPosition.theta_n,ViewFactor,ParVegTree)
SWRin_t = SWRin_nT
SWRout_t = SWRout_nT
SWRabs_t = SWRabs_nT
SWRabsDir_t = SWRabsDir_nT
SWRabsDiff_t = SWRabsDiff_nT
SWREB_t = SWREB_nT
return SWRin_t,SWRout_t,SWRabs_t,SWRabsDir_t,SWRabsDiff_t,SWREB_t
def LWRabsorbedNoTree(self,h_can,w_can,LWR,fgveg,fgbare,fgimp,ew,egveg,egbare,egimp,Tgimp,Tgbare,Tgveg,Twsun,Twshade,ViewFactor):
F_gs_nT = ViewFactor.F_gs_nT
F_gw_nT = ViewFactor.F_gw_nT
F_ww_nT = ViewFactor.F_ww_nT
F_wg_nT = ViewFactor.F_wg_nT
F_ws_nT = ViewFactor.F_ws_nT
F_sg_nT = ViewFactor.F_sg_nT
F_sw_nT = ViewFactor.F_sw_nT
# normalized surface areas
A_s = w_can
A_g = w_can
A_w = h_can
bolzm = 5.67 * 10 ** (-8)
SVF = numpy.zeros(3)
SVF[0] = F_gs_nT + 2 * F_gw_nT
SVF[1] = F_ww_nT + F_wg_nT + F_ws_nT
SVF[2] = F_sg_nT + 2 * F_sw_nT
SVF2 = numpy.zeros(3)
SVF2[0] = F_gs_nT + 2 * F_ws_nT * h_can
SVF2[1] = F_sg_nT + 2 * F_wg_nT * h_can
SVF2[2] = F_ww_nT + F_sw_nT / h_can + F_gw_nT / h_can
for i in range(0,len(SVF)):
if SVF[i] < 0.999 or SVF[i] > 1.001:
print('The view factors do not add up to 1 for a canyon with trees')
# Solve for infinite reflections equation A*X=C
if fgimp > 0:
Cimp = 1
else:
Cimp = 0
if fgbare > 0:
Cbare = 1
else:
Cbare = 0
if fgveg > 0:
Cveg = 1
else:
Cveg = 0
# View factor matrix to solve for infinite reflections equation
Tij = numpy.array([[1,0,0, -(1-egveg)*F_gw_nT*Cveg, -(1-egveg)*F_gw_nT*Cveg, -(1-egveg)*F_gs_nT*Cveg],
[0,1,0, -(1-egbare)*F_gw_nT*Cbare, -(1-egbare)*F_gw_nT*Cbare, -(1-egbare)*F_gs_nT*Cbare],
[0,0,1, -(1-egimp)*F_gw_nT*Cimp, -(1-egimp)*F_gw_nT*Cimp, -(1-egimp)*F_gs_nT*Cimp],
[-(1-ew)*F_wg_nT*fgveg*Cveg, -(1-ew)*F_wg_nT*fgbare*Cbare, -(1-ew)*F_wg_nT*fgimp*Cimp, 1, -(1-ew)*F_ww_nT, -(1-ew)*F_ws_nT],
[-(1-ew)*F_wg_nT*fgveg*Cveg, -(1-ew)*F_wg_nT*fgbare*Cbare, -(1-ew)*F_wg_nT*fgimp*Cimp, -(1-ew)*F_ww_nT, 1, -(1-ew)*F_ws_nT],
[0, 0, 0, 0, 0, 1]])
# Emitted radiation per surface
Omega_i = numpy.array([(egveg*bolzm*(Tgveg)**4*Cveg),(egbare * bolzm * (Tgbare) ** 4 * Cbare),(egimp * bolzm * (Tgimp) ** 4 * Cimp),(ew * bolzm * (Twsun) ** 4),(ew * bolzm * (Twshade) ** 4),LWR])
# Outgoing radiation per surface
# Outgoing radiation [W/m^2] per m^2 surface area
B_i = numpy.linalg.solve(Tij,Omega_i)
if B_i[5] != LWR:
print('Incoming lonwave radiation and emitted longwave radiation from the sky after the matrix inversion are not equal')
# Incoming longwave radiation at each surface A_i
Tij2 = numpy.array([[0, 0, 0, F_gw_nT*Cveg, F_gw_nT*Cveg, F_gs_nT*Cveg],
[0, 0, 0, F_gw_nT*Cbare, F_gw_nT*Cbare, F_gs_nT*Cbare],
[0, 0, 0, F_gw_nT*Cimp, F_gw_nT*Cimp, F_gs_nT*Cimp],
[F_wg_nT*fgveg*Cveg, F_wg_nT*fgbare*Cbare, F_wg_nT*fgimp*Cimp, 0, F_ww_nT, F_ws_nT],
[F_wg_nT*fgveg*Cveg, F_wg_nT*fgbare*Cbare, F_wg_nT*fgimp*Cimp, F_ww_nT, 0, F_ws_nT],
[0, 0, 0, 0, 0, 0]])
A_i = numpy.dot(Tij2,B_i)
e_i = [egveg,egbare,egimp,ew,ew,0]
A_i2 = [(B_i[i]-Omega_i[i])/(1-e_i[i]) for i in range(0,len(e_i))]
Qnet_i2 = A_i - B_i
# Absorbed longwave radiation
e_i = [egveg, egbare, egimp, ew, ew, 0]
Qnet_i = [(e_i[i] * B_i[i] - Omega_i[i]) / (1 - e_i[i]) for i in range(0,len(e_i))]
for i in range(0,len(e_i)):
if e_i[i] == 1:
Qnet_i[i] = A_i[i] - Omega_i[i]
# Assumption: The sky has a fixed emission of LWR. Hence, Qnet is 0
Qnet_i[5] = 0
# Assignment
LWRout_i = B_i # Outgoing radiation [W/m^2] per m^2 surface area
LWRemit_i = Omega_i # Emitted radiation [W/m^2] per m^2 surface area
LWRin_i = A_i # Incoming radiation [W/m^2] per m^2 surface area
LWRnet_i = Qnet_i # Net absorbed radiation [W/m^2] per m^2 surface area
# Energy Balance
LWRin_atm = LWR
TotalLWRSurface_in = LWRin_i[0]*fgveg*A_g/A_g + LWRin_i[1]*fgbare*A_g/A_g + LWRin_i[2]*fgimp*A_g/A_g +\
LWRin_i[3]*A_w/A_g + LWRin_i[4]*A_w/A_g
TotalLWRSurface_abs = LWRnet_i[0]*fgveg*A_g/A_g + LWRnet_i[1]*fgbare*A_g/A_g + LWRnet_i[2]*fgimp*A_g/A_g +\
LWRnet_i[3]*A_w/A_g + LWRnet_i[4]*A_w/A_g
TotalLWRSurface_out = LWRout_i[0]*fgveg*A_g/A_s+LWRout_i[1]*fgbare*A_g/A_s+LWRout_i[2]*fgimp*A_g/A_s +\
LWRout_i[3]*A_w/A_s+LWRout_i[4]*A_w/A_s
TotalLWRref_to_atm = LWRout_i[0]*F_sg_nT*fgveg + LWRout_i[1]*F_sg_nT*fgbare + LWRout_i[2]*F_sg_nT*fgimp + \
LWRout_i[3]*F_sw_nT + LWRout_i[4]*F_sw_nT
EBSurface = TotalLWRSurface_in - TotalLWRSurface_abs - TotalLWRSurface_out
EBCanyon = LWRin_atm - TotalLWRSurface_abs - TotalLWRref_to_atm
# Energy balance
if abs(EBSurface)>=10**(-1):
print('EBSurface is not 0',EBSurface)
if abs(EBCanyon)>=10**(-1):
print('EBCanyon is not 0',EBCanyon)
class LWRin_nT_Def():
pass
LWRin_nT = LWRin_nT_Def()
LWRin_nT.LWRinGroundImp = LWRin_i[2] * Cimp
LWRin_nT.LWRinGroundBare = LWRin_i[1] * Cbare
LWRin_nT.LWRinGroundVeg = LWRin_i[0] * Cveg
LWRin_nT.LWRinTree = 0
LWRin_nT.LWRinWallSun = LWRin_i[3]
LWRin_nT.LWRinWallShade = LWRin_i[4]
LWRin_nT.LWRinTotalGround = fgveg * LWRin_i[0] + fgbare * LWRin_i[1] + fgimp * LWRin_i[2]
LWRin_nT.LWRinTotalCanyon = LWRin_i[0] * fgveg * A_g / A_g + LWRin_i[1] * fgbare * A_g / A_g + \
LWRin_i[2] * fgimp * A_g / A_g + LWRin_i[3] * A_w / A_g + LWRin_i[4] * A_w / A_g
# Outgoing longwave radiation
class LWRout_nT_Def():
pass
LWRout_nT = LWRout_nT_Def()
LWRout_nT.LWRoutGroundImp = LWRout_i[2] * Cimp
LWRout_nT.LWRoutGroundBare = LWRout_i[1] * Cbare
LWRout_nT.LWRoutGroundVeg = LWRout_i[0] * Cveg
LWRout_nT.LWRoutTree = 0
LWRout_nT.LWRoutWallSun = LWRout_i[3]
LWRout_nT.LWRoutWallShade = LWRout_i[4]
LWRout_nT.LWRoutTotalGround = fgveg * LWRout_i[0] + fgbare * LWRout_i[1] + fgimp * LWRout_i[2]
LWRout_nT.LWRoutTotalCanyon = LWRout_i[0] * fgveg * A_g / A_g + LWRout_i[1] * fgbare * A_g / A_g + \
LWRout_i[2] * fgimp * A_g / A_g + LWRout_i[3] * A_w / A_g + LWRout_i[4] * A_w / A_g
# Absorbed longwave radiation
class LWRabs_nT_Def():
pass
LWRabs_nT = LWRabs_nT_Def()
LWRabs_nT.LWRabsGroundImp = LWRnet_i[2] * Cimp
LWRabs_nT.LWRabsGroundBare = LWRnet_i[1] * Cbare
LWRabs_nT.LWRabsGroundVeg = LWRnet_i[0] * Cveg
LWRabs_nT.LWRabsTree = 0
LWRabs_nT.LWRabsWallSun = LWRnet_i[3]
LWRabs_nT.LWRabsWallShade = LWRnet_i[4]
LWRabs_nT.LWRabsTotalGround = fgveg * LWRnet_i[0] + fgbare * LWRnet_i[1] + fgimp * LWRnet_i[2]
LWRabs_nT.LWRabsTotalCanyon = LWRnet_i[0] * fgveg * A_g / A_g + LWRnet_i[1] * fgbare * A_g / A_g + \
LWRnet_i[2] * fgimp * A_g / A_g + LWRnet_i[3] * A_w / A_g + LWRnet_i[4] * A_w / A_g
# Energy Balance of longwave radiation
class LWREB_nT_Def():
pass
LWREB_nT = LWREB_nT_Def()
LWREB_nT.LWREBGroundImp = LWRin_nT.LWRinGroundImp - LWRout_nT.LWRoutGroundImp - LWRabs_nT.LWRabsGroundImp
LWREB_nT.LWREBGroundBare = LWRin_nT.LWRinGroundBare - LWRout_nT.LWRoutGroundBare - LWRabs_nT.LWRabsGroundBare
LWREB_nT.LWREBGroundVeg = LWRin_nT.LWRinGroundVeg - LWRout_nT.LWRoutGroundVeg - LWRabs_nT.LWRabsGroundVeg
LWREB_nT.LWREBTree = 0
LWREB_nT.LWREBWallSun = LWRin_nT.LWRinWallSun - LWRout_nT.LWRoutWallSun - LWRabs_nT.LWRabsWallSun
LWREB_nT.LWREBWallShade = LWRin_nT.LWRinWallShade - LWRout_nT.LWRoutWallShade - LWRabs_nT.LWRabsWallShade
LWREB_nT.LWREBTotalGround = LWRin_nT.LWRinTotalGround - LWRout_nT.LWRoutTotalGround - LWRabs_nT.LWRabsTotalGround
LWREB_nT.LWREBTotalCanyon = LWRin_nT.LWRinTotalCanyon - LWRout_nT.LWRoutTotalCanyon - LWRabs_nT.LWRabsTotalCanyon
if abs(LWREB_nT.LWREBGroundImp)>=10**(-6):
print('LWREB_nT.LWREBGroundImp is not 0.')
if abs(LWREB_nT.LWREBGroundBare)>=10**(-6):
print('LWREB_nT.LWREBGroundBare is not 0.')
if abs(LWREB_nT.LWREBGroundVeg )>=10**(-6):
print('LWREB_nT.LWREBGroundVeg is not 0.')
if abs(LWREB_nT.LWREBWallSun)>=10**(-6):
print('LWREB_nT.LWREBWallSun is not 0.')
if abs(LWREB_nT.LWREBWallShade)>=10**(-6):
print('LWREB_nT.LWREBWallShade is not 0.')
if abs(LWREB_nT.LWREBTotalGround)>=10**(-6):
print('LWREB_nT.LWREBTotalGround is not 0.')
if abs(LWREB_nT.LWREBTotalCanyon)>=10**(-6):
print('LWREB_nT.LWREBTotalCanyon is not 0.')
return LWRin_nT,LWRout_nT,LWRabs_nT,LWREB_nT
def LWRabsorbedWithTrees(self,h_can,w_can,r_tree,LWR,fgveg,fgbare,fgimp,ew,et,egveg,egbare,egimp,Tgimp,Tgbare,Tgveg,
Twsun,Twshade,Ttree,ViewFactor):
F_gs_T = ViewFactor.F_gs_T
F_gt_T = ViewFactor.F_gt_T
F_gw_T = ViewFactor.F_gw_T
F_ww_T = ViewFactor.F_ww_T
F_wt_T = ViewFactor.F_wt_T
F_wg_T = ViewFactor.F_wg_T
F_ws_T = ViewFactor.F_ws_T
F_sg_T = ViewFactor.F_sg_T
F_sw_T = ViewFactor.F_sw_T
F_st_T = ViewFactor.F_st_T
F_tg_T = ViewFactor.F_tg_T
F_tw_T = ViewFactor.F_tw_T
F_ts_T = ViewFactor.F_ts_T
F_tt_T = ViewFactor.F_tt_T
# normalized surface areas
A_s = w_can
A_g = w_can
A_w = h_can
# There are 2 trees. Hence, the area of tree is twice a circle
A_t = 2 * 2 * numpy.pi * r_tree
bolzm = 5.67 * 10 ** (-8)
# Check if view factors add up to 1
SVF = numpy.zeros(4)
SVF[0] = F_gs_T + F_gt_T + 2 * F_gw_T
SVF[1] = F_ww_T + F_wt_T + F_wg_T + F_ws_T
SVF[2] = F_sg_T + 2 * F_sw_T + F_st_T
SVF[3] = F_ts_T + 2 * F_tw_T + F_tt_T + F_tg_T
SVF2 = numpy.zeros(4)
SVF2[0] = F_gs_T + 2 * F_ws_T + F_ts_T
SVF2[1] = F_sg_T + 2 * F_wg_T + F_tg_T
SVF2[2] = F_ww_T + F_sw_T + F_gw_T + F_tw_T
SVF2[3] = F_gt_T + 2 * F_wt_T + F_tt_T + F_st_T
for i in range(0,len(SVF)):
if SVF[i]<0.999 or SVF[i]>1.001:
print('The view factors do not add up to 1 for a canyon with trees')
# Solve for infinite reflections equation A*X=C
if fgimp > 0:
Cimp = 1
else:
Cimp = 0
if fgbare > 0:
Cbare = 1
else:
Cbare = 0
if fgveg > 0:
Cveg = 1
else:
Cveg = 0
Tij = numpy.array([[1,0,0, -(1-egveg)*F_gw_T*Cveg, -(1-egveg)*F_gw_T*Cveg, -(1-egveg)*F_gt_T*Cveg, -(1-egveg)*F_gs_T*Cveg],
[0,1,0, -(1-egbare)*F_gw_T*Cbare, -(1-egbare)*F_gw_T*Cbare, -(1-egbare)*F_gt_T*Cbare, -(1-egbare)*F_gs_T*Cbare],
[0,0,1, -(1-egimp)*F_gw_T*Cimp, -(1-egimp)*F_gw_T*Cimp, -(1-egimp)*F_gt_T*Cimp, -(1-egimp)*F_gs_T*Cimp],
[-(1-ew)*F_wg_T*fgveg*Cveg, -(1-ew)*F_wg_T*fgbare*Cbare, -(1-ew)*F_wg_T*fgimp*Cimp, 1, -(1-ew)*F_ww_T, -(1-ew)*F_wt_T, -(1-ew)*F_ws_T],
[-(1-ew)*F_wg_T*fgveg*Cveg, -(1-ew)*F_wg_T*fgbare*Cbare, -(1-ew)*F_wg_T*fgimp*Cimp, -(1-ew)*F_ww_T, 1, -(1-ew)*F_wt_T, -(1-ew)*F_ws_T],
[-(1-et)*F_tg_T*fgveg*Cveg, -(1-et)*F_tg_T*fgbare*Cbare, -(1-et)*F_tg_T*fgimp*Cimp, -(1-et)*F_tw_T, -(1-et)*F_tw_T, 1-(1-et)*F_tt_T, -(1-et)*F_ts_T],
[0, 0, 0, 0, 0, 0, 1]])
Omega_i = numpy.array([(egveg*bolzm*(Tgveg)**4*Cveg),
(egbare * bolzm * (Tgbare) ** 4 * Cbare),
(egimp * bolzm * (Tgimp) ** 4 * Cimp),
(ew * bolzm * (Twsun) ** 4),
(ew * bolzm * (Twshade) ** 4),
(et * bolzm * (Ttree) ** 4),
LWR])
# Outgoing radiation per surface
# Outgoing radiation [W/m^2] per m^2 surface area
B_i = numpy.linalg.solve(Tij,Omega_i)
if B_i[6] != LWR:
print('Incoming lonwave radiation and emitted longwave radiation from the sky after the matrix inversion are not equal')
Tij2 = numpy.array([[0, 0, 0, F_gw_T*Cveg, F_gw_T*Cveg, F_gt_T*Cveg, F_gs_T*Cveg],
[0, 0, 0, F_gw_T*Cbare, F_gw_T*Cbare, F_gt_T*Cbare, F_gs_T*Cbare],
[0, 0, 0, F_gw_T*Cimp, F_gw_T*Cimp, F_gt_T*Cimp, F_gs_T*Cimp],
[F_wg_T*fgveg*Cveg, F_wg_T*fgbare*Cbare, F_wg_T*fgimp*Cimp, 0, F_ww_T, F_wt_T, F_ws_T],
[F_wg_T*fgveg*Cveg, F_wg_T*fgbare*Cbare, F_wg_T*fgimp*Cimp, F_ww_T, 0, F_wt_T, F_ws_T],
[F_tg_T*fgveg*Cveg, F_tg_T*fgbare*Cbare, F_tg_T*fgimp*Cimp, F_tw_T, F_tw_T, F_tt_T, F_ts_T],
[0, 0, 0, 0, 0, 0, 0]])
A_i = numpy.dot(Tij2, B_i)
e_i = [egveg, egbare, egimp, ew, ew, et, 0]
A_i2 = [(B_i[i] - Omega_i[i]) / (1 - e_i[i]) for i in range(0,len(e_i))]
Qnet_i2 = A_i - B_i
# Absorbed longwave radiation
e_i = [egveg, egbare, egimp, ew, ew, et, 0]
Qnet_i = [(e_i[i] * B_i[i] - Omega_i[i]) / (1 - e_i[i]) for i in range(0,len(e_i))]
Qnet_i = [(e_i[i] * B_i[i] - Omega_i[i]) / (1 - e_i[i]) for i in range(0, len(e_i))]
for i in range(0, len(e_i)):
if e_i[i] == 1:
Qnet_i[i] = A_i[i] - Omega_i[i]
# Assumption: The sky has a fixed emission of LWR. Hence, Qnet is 0
Qnet_i[6] = 0
# Assignment
LWRout_i = B_i # Outgoing radiation [W/m^2] per m^2 surface area
LWRemit_i = Omega_i # Emitted radiation [W/m^2] per m^2 surface area
LWRin_i = A_i # Incoming radiation [W/m^2] per m^2 surface area
LWRnet_i = Qnet_i # Net absorbed radiation [W/m^2] per m^2 surface area
# Energy balance
LWRin_atm = LWR
TotalLWRSurface_in = LWRin_i[0]*fgveg*A_g/A_g + LWRin_i[1]*fgbare*A_g/A_g + LWRin_i[2]*fgimp*A_g/A_g +\
LWRin_i[3]*A_w/A_g + LWRin_i[4]*A_w/A_g + LWRin_i[5]*A_t/A_g
TotalLWRSurface_abs = LWRnet_i[0]*fgveg*A_g/A_g + LWRnet_i[1]*fgbare*A_g/A_g + LWRnet_i[2]*fgimp*A_g/A_g +\
LWRnet_i[3]*A_w/A_g + LWRnet_i[4]*A_w/A_g + LWRnet_i[5]*A_t/A_g
TotalLWRSurface_out = LWRout_i[0]*fgveg*A_g/A_s+LWRout_i[1]*fgbare*A_g/A_s+LWRout_i[2]*fgimp*A_g/A_s +\
LWRout_i[3]*A_w/A_s+LWRout_i[4]*A_w/A_s+LWRout_i[5]*A_t/A_s
TotalLWRref_to_atm = LWRout_i[0]*F_sg_T*fgveg + LWRout_i[1]*F_sg_T*fgbare + LWRout_i[2]*F_sg_T*fgimp + \
LWRout_i[3]*F_sw_T + LWRout_i[4]*F_sw_T + LWRout_i[5]*F_st_T
EBSurface = TotalLWRSurface_in - TotalLWRSurface_abs - TotalLWRSurface_out
EBCanyon = LWRin_atm - TotalLWRSurface_abs - TotalLWRref_to_atm
# Energy balance
if abs(EBSurface)>=10**(-1):
print('EBSurface is not 0',EBSurface)
if abs(EBCanyon)>=10**(-1):
print('EBCanyon is not 0',EBCanyon)
class LWRin_T_Def():
pass
LWRin_T = LWRin_T_Def()
LWRin_T.LWRinGroundImp = LWRin_i[2] * Cimp
LWRin_T.LWRinGroundBare = LWRin_i[1] * Cbare
LWRin_T.LWRinGroundVeg = LWRin_i[0] * Cveg
LWRin_T.LWRinTree = LWRin_i[5]
LWRin_T.LWRinWallSun = LWRin_i[3]
LWRin_T.LWRinWallShade = LWRin_i[4]
LWRin_T.LWRinTotalGround = fgveg * LWRin_i[0] + fgbare * LWRin_i[1] + fgimp * LWRin_i[2]
LWRin_T.LWRinTotalCanyon = LWRin_i[0]*fgveg*A_g/A_g + LWRin_i[1]*fgbare*A_g/A_g + \
LWRin_i[2]*fgimp*A_g/A_g + LWRin_i[3]*A_w/A_g + LWRin_i[4]*A_w/A_g + LWRin_i[5]*A_t/A_g
# Outgoing longwave radiation
class LWRout_T_Def():
pass
LWRout_T = LWRout_T_Def()
LWRout_T.LWRoutGroundImp = LWRout_i[2] * Cimp
LWRout_T.LWRoutGroundBare = LWRout_i[1] * Cbare
LWRout_T.LWRoutGroundVeg = LWRout_i[0] * Cveg
LWRout_T.LWRoutTree = LWRout_i[5]
LWRout_T.LWRoutWallSun = LWRout_i[3]
LWRout_T.LWRoutWallShade = LWRout_i[4]
LWRout_T.LWRoutTotalGround = fgveg * LWRout_i[0] + fgbare * LWRout_i[1] + fgimp * LWRout_i[2]
LWRout_T.LWRoutTotalCanyon = LWRout_i[0] * fgveg * A_g / A_g + LWRout_i[1] * fgbare * A_g / A_g + \
LWRout_i[2] * fgimp * A_g / A_g + LWRout_i[3] * A_w / A_g + LWRout_i[4] * A_w / A_g + LWRout_i[5]*A_t/A_g
# Absorbed longwave radiation
class LWRabs_T_Def():
pass
LWRabs_T = LWRabs_T_Def()
LWRabs_T.LWRabsGroundImp = LWRnet_i[2] * Cimp
LWRabs_T.LWRabsGroundBare = LWRnet_i[1] * Cbare
LWRabs_T.LWRabsGroundVeg = LWRnet_i[0] * Cveg
LWRabs_T.LWRabsTree = LWRnet_i[5]
LWRabs_T.LWRabsWallSun = LWRnet_i[3]
LWRabs_T.LWRabsWallShade = LWRnet_i[4]
LWRabs_T.LWRabsTotalGround = fgveg * LWRnet_i[0] + fgbare * LWRnet_i[1] + fgimp * LWRnet_i[2]
LWRabs_T.LWRabsTotalCanyon = LWRnet_i[0] * fgveg * A_g / A_g + LWRnet_i[1] * fgbare * A_g / A_g + \
LWRnet_i[2] * fgimp * A_g / A_g + LWRnet_i[3] * A_w / A_g + LWRnet_i[4] * A_w / A_g + LWRnet_i[5]*A_t/A_g
# Energy Balance of longwave radiation
class LWREB_T_Def():
pass
LWREB_T = LWREB_T_Def()
LWREB_T.LWREBGroundImp = LWRin_T.LWRinGroundImp - LWRout_T.LWRoutGroundImp - LWRabs_T.LWRabsGroundImp
LWREB_T.LWREBGroundBare = LWRin_T.LWRinGroundBare - LWRout_T.LWRoutGroundBare - LWRabs_T.LWRabsGroundBare
LWREB_T.LWREBGroundVeg = LWRin_T.LWRinGroundVeg - LWRout_T.LWRoutGroundVeg - LWRabs_T.LWRabsGroundVeg
LWREB_T.LWREBTree = LWRin_T.LWRinTree - LWRout_T.LWRoutTree - LWRabs_T.LWRabsTree
LWREB_T.LWREBWallSun = LWRin_T.LWRinWallSun - LWRout_T.LWRoutWallSun - LWRabs_T.LWRabsWallSun
LWREB_T.LWREBWallShade = LWRin_T.LWRinWallShade - LWRout_T.LWRoutWallShade - LWRabs_T.LWRabsWallShade
LWREB_T.LWREBTotalGround = LWRin_T.LWRinTotalGround - LWRout_T.LWRoutTotalGround - LWRabs_T.LWRabsTotalGround
LWREB_T.LWREBTotalCanyon = LWRin_T.LWRinTotalCanyon - LWRout_T.LWRoutTotalCanyon - LWRabs_T.LWRabsTotalCanyon
if abs(LWREB_T.LWREBGroundImp)>=10**(-6):
print('LWREB_T.LWREBGroundImp is not 0.')
if abs(LWREB_T.LWREBGroundBare)>=10**(-6):
print('LWREB_T.LWREBGroundBare is not 0.')
if abs(LWREB_T.LWREBGroundVeg )>=10**(-6):
print('LWREB_T.LWREBGroundVeg is not 0.')
if abs(LWREB_T.LWREBWallSun)>=10**(-6):
print('LWREB_T.LWREBWallSun is not 0.')
if abs(LWREB_T.LWREBWallShade)>=10**(-6):
print('LWREB_T.LWREBWallShade is not 0.')
if abs(LWREB_T.LWREBTotalGround)>=10**(-6):
print('LWREB_T.LWREBTotalGround is not 0.')
if abs(LWREB_T.LWREBTotalCanyon)>=10**(-6):
print('LWREB_T.LWREBTotalCanyon is not 0.')
if abs(LWREB_T.LWREBTree)>=10**(-6):
print('LWREB_T.LWREBTree is not 0.')
return LWRin_T,LWRout_T,LWRabs_T,LWREB_T
def SWRabsorbedNoTrees(self,h_can,w_can,fgveg,fgbare,fgimp,aw,agveg,agbare,agimp,SWR_dir,SWR_diff,theta_Z,theta_n,
ViewFactor,ParVegTree):
F_gs_nT = ViewFactor.F_gs_nT
F_gw_nT = ViewFactor.F_gw_nT
F_ww_nT = ViewFactor.F_ww_nT
F_wg_nT = ViewFactor.F_wg_nT
F_ws_nT = ViewFactor.F_ws_nT
F_sg_nT = ViewFactor.F_sg_nT
F_sw_nT = ViewFactor.F_sw_nT
# normalized surface areas
A_s = w_can
A_g = w_can
A_w = h_can
SVF = numpy.zeros(3)
SVF[0] = F_gs_nT + 2 * F_gw_nT
SVF[1] = F_ww_nT + F_wg_nT + F_ws_nT
SVF[2] = F_sg_nT + 2 * F_sw_nT
SVF2 = numpy.zeros(3)
SVF2[0] = F_gs_nT + 2 * F_ws_nT * h_can
SVF2[1] = F_sg_nT + 2 * F_wg_nT * h_can
SVF2[2] = F_ww_nT + F_sw_nT / h_can + F_gw_nT / h_can
for i in range(0, len(SVF)):
if SVF[i] < 0.999 or SVF[i] > 1.001:
print('The view factors do not add up to 1 for a canyon with trees')
SWRdir_ground, SWRdir_wallsun, SWRdir_wallshade, NOTUSED = self.DirectSWRSurfaces(h_can, w_can, numpy.nan, numpy.nan, numpy.nan,
theta_Z,theta_n, SWR_dir, numpy.nan,
0,ParVegTree)
# Balance direct shortwave radiation in
EB_SWRdir = SWR_dir - (SWRdir_ground * A_g / A_g + SWRdir_wallsun * A_w / A_g + SWRdir_wallshade * A_w / A_g)
EB_SWRdiff = SWR_diff - (F_sg_nT * SWR_diff + F_sw_nT * SWR_diff + F_sw_nT * SWR_diff)
if abs(EB_SWRdir)>=10**(-1):
print('EB_SWRdir is not 0',EB_SWRdir)
if abs(EB_SWRdiff)>=10**(-1):
print('EB_SWRdiff is not 0',EB_SWRdiff)
if fgimp > 0:
Cimp = 1
else:
Cimp = 0
if fgbare > 0:
Cbare = 1
else:
Cbare = 0
if fgveg > 0:
Cveg = 1
else:
Cveg = 0
ai = [agveg,agbare,agimp,aw,aw,0]
# View factor matrix to solve for infinite reflections equation
Tij = numpy.array([[1,0,0, -agveg*F_gw_nT*Cveg, -agveg*F_gw_nT*Cveg, -agveg*F_gs_nT*Cveg],
[0,1,0, -agbare*F_gw_nT*Cbare, -agbare*F_gw_nT*Cbare, -agbare*F_gs_nT*Cbare],
[0,0,1, -agimp*F_gw_nT*Cimp, -agimp*F_gw_nT*Cimp, -agimp*F_gs_nT*Cimp],
[-aw*F_wg_nT*fgveg*Cveg,-aw*F_wg_nT*fgbare*Cbare,-aw*F_wg_nT*fgimp*Cimp, 1, -aw*F_ww_nT, -aw*F_ws_nT],
[-aw*F_wg_nT*fgveg*Cveg,-aw*F_wg_nT*fgbare*Cbare,-aw*F_wg_nT*fgimp*Cimp, -aw*F_ww_nT, 1, -aw*F_ws_nT],
[0, 0, 0, 0, 0, 1]])
# Incoming shortwave radiation from sky
Omega_i = numpy.array([agveg*SWRdir_ground*Cveg,
agbare * SWRdir_ground * Cbare,
agimp * SWRdir_ground * Cimp,
aw * SWRdir_wallsun,
aw * 0,
SWR_diff])
# Outgoing radiation per surface
B_i = numpy.linalg.solve(Tij,Omega_i)
if B_i[5] != SWR_diff:
print('Incoming lonwave radiation and emitted longwave radiation from the sky after the matrix inversion are not equal')
# Incoming shortwave radiation at each surface A_i
Tij2 = numpy.array([[0, 0, 0, F_gw_nT*Cveg, F_gw_nT*Cveg, F_gs_nT*Cveg],
[0, 0, 0, F_gw_nT*Cbare, F_gw_nT*Cbare, F_gs_nT*Cbare],
[0, 0, 0, F_gw_nT*Cimp, F_gw_nT*Cimp, F_gs_nT*Cimp],
[F_wg_nT*fgveg*Cveg, F_wg_nT*fgbare*Cbare, F_wg_nT*fgimp*Cimp, 0, F_ww_nT, F_ws_nT],
[F_wg_nT*fgveg*Cveg, F_wg_nT*fgbare*Cbare, F_wg_nT*fgimp*Cimp, F_ww_nT, 0, F_ws_nT],
[0, 0, 0, 0, 0, 0]])
SWRdir_i = numpy.array([SWRdir_ground*Cveg,
SWRdir_ground * Cbare,
SWRdir_ground * Cimp,
SWRdir_wallsun,
0,
0])
A_i1 = numpy.dot(Tij2,B_i)+SWRdir_i # Incoming radiation [W/m^2] per m^2 surface area
A_i = B_i/ai # Incoming radiation [W/m^2] per m^2 surface area
for i in range(0,len(ai)):
if ai[i] == 0:
A_i[i] = A_i1[i]
A_i[5] = 0 # Assumption: The sky has a fixed emission of LWR. Hence, Qnet is 0.
# Absorbed shortwave radiation at ech surface Qnet_i
Qnet_i = A_i-B_i
# Assignment
SWRout_i = B_i # Outgoing radiation [W/m^2] per m^2 surface area
SWRin_i = A_i # Incoming radiation [W/m^2] per m^2 surface area
SWRnet_i = Qnet_i # Net absorbed radiation [W/m^2] per m^2 surface area
# Energy balance
SWRin_atm = SWR_dir + SWR_diff
TotalSWRSurface_in = SWRin_i[0] * fgveg * A_g / A_g + SWRin_i[1] * fgbare * A_g / A_g + SWRin_i[2] * fgimp * A_g / A_g +\
SWRin_i[3] * A_w / A_g + SWRin_i[4] * A_w / A_g
TotalSWRSurface_abs = SWRnet_i[0] * fgveg * A_g / A_g + SWRnet_i[1] * fgbare * A_g / A_g + SWRnet_i[2] * fgimp * A_g / A_g +\
SWRnet_i[3] * A_w / A_g + SWRnet_i[4] * A_w / A_g
TotalSWRSurface_out = SWRout_i[0] * fgveg * A_g / A_s + SWRout_i[1] * fgbare * A_g / A_s + SWRout_i[2] * fgimp * A_g / A_s + \
SWRout_i[3] * A_w / A_s + SWRout_i[4] * A_w / A_s
TotalSWRref_to_atm = SWRout_i[0] * F_sg_nT * fgveg + SWRout_i[1] * F_sg_nT * fgbare + SWRout_i[2] * F_sg_nT * fgimp +\
SWRout_i[3] * F_sw_nT + SWRout_i[4] * F_sw_nT
EBSurface = TotalSWRSurface_in - TotalSWRSurface_abs - TotalSWRSurface_out
EBCanyon = SWRin_atm - TotalSWRSurface_abs - TotalSWRref_to_atm
if abs(EBSurface)>=10**(-1):
print('EBSurface is not 0',EBSurface)
if abs(EBCanyon)>=10**(-1):
print('EBCanyon is not 0',EBCanyon)
# Incoming shortwave radiation
class SWRin_nT_Def():
pass
SWRin_nT = SWRin_nT_Def()
SWRin_nT.SWRinGroundImp = SWRin_i[2] * Cimp
SWRin_nT.SWRinGroundBare = SWRin_i[1] * Cbare
SWRin_nT.SWRinGroundVeg = SWRin_i[0] * Cveg
SWRin_nT.SWRinTree = 0
SWRin_nT.SWRinWallSun = SWRin_i[3]
SWRin_nT.SWRinWallShade = SWRin_i[4]
SWRin_nT.SWRinTotalGround = fgveg * SWRin_i[0] + fgbare * SWRin_i[1] + fgimp * SWRin_i[2]
SWRin_nT.SWRinTotalCanyon = SWRin_i[0] * fgveg * A_g / A_g + SWRin_i[1] * fgbare * A_g / A_g + SWRin_i[2] * fgimp * A_g / A_g + \
SWRin_i[3] * A_w / A_g + SWRin_i[4] * A_w / A_g
# Outgoing shortwave radiation
class SWRout_nT_Def():
pass
SWRout_nT = SWRout_nT_Def()
SWRout_nT.SWRoutGroundImp = SWRout_i[2] * Cimp
SWRout_nT.SWRoutGroundBare = SWRout_i[1] * Cbare
SWRout_nT.SWRoutGroundVeg = SWRout_i[0] * Cveg
SWRout_nT.SWRoutTree = 0
SWRout_nT.SWRoutWallSun = SWRout_i[3]
SWRout_nT.SWRoutWallShade = SWRout_i[4]
SWRout_nT.SWRoutTotalGround = fgveg * SWRout_i[0] + fgbare * SWRout_i[1] + fgimp * SWRout_i[2]
SWRout_nT.SWRoutTotalCanyon = SWRout_i[0] * fgveg * A_g / A_g + SWRout_i[1] * fgbare * A_g / A_g + SWRout_i[2] * fgimp * A_g / A_g + \
SWRout_i[3] * A_w / A_g + SWRout_i[4] * A_w / A_g
# Absorbed shortwave radiation
class SWRabs_nT_Def():
pass
SWRabs_nT = SWRabs_nT_Def()
SWRabs_nT.SWRabsGroundImp = SWRnet_i[2] * Cimp
SWRabs_nT.SWRabsGroundBare = SWRnet_i[1] * Cbare
SWRabs_nT.SWRabsGroundVeg = SWRnet_i[0] * Cveg
SWRabs_nT.SWRabsTree = 0
SWRabs_nT.SWRabsWallSun = SWRnet_i[3]
SWRabs_nT.SWRabsWallShade = SWRnet_i[4]
SWRabs_nT.SWRabsTotalGround = fgveg * SWRnet_i[0] + fgbare * SWRnet_i[1] + fgimp * SWRnet_i[2]
SWRabs_nT.SWRabsTotalCanyon = SWRnet_i[0] * fgveg * A_g / A_g + SWRnet_i[1] * fgbare * A_g / A_g + SWRnet_i[2] * fgimp * A_g / A_g + \
SWRnet_i[3] * A_w / A_g + SWRnet_i[4] * A_w / A_g
# Direct absorbed shortwave radiation
class SWRabsDir_nT_Def():
pass
SWRabsDir_nT = SWRabsDir_nT_Def()
SWRabsDir_nT.SWRabsGroundImp = (1-agimp)*SWRdir_ground*Cimp
SWRabsDir_nT.SWRabsGroundBare = (1-agbare)*SWRdir_ground*Cbare
SWRabsDir_nT.SWRabsGroundVeg = (1-agveg)*SWRdir_ground*Cveg
SWRabsDir_nT.SWRabsTree = 0
SWRabsDir_nT.SWRabsWallSun = (1-aw)*SWRdir_wallsun
SWRabsDir_nT.SWRabsWallShade = (1-aw)*SWRdir_wallshade
SWRabsDir_nT.SWRabsTotalGround = fgveg*(1-agveg)*SWRdir_ground+fgbare*(1-agbare)*SWRdir_ground+fgimp*(1-agimp)*SWRdir_ground
SWRabsDir_nT.SWRabsTotalCanyon = fgveg*(1-agveg)*SWRdir_ground*A_g/A_g+fgbare*(1-agbare)*SWRdir_ground*A_g/A_g+\
fgimp*(1-agimp)*SWRdir_ground*A_g/A_g + (1-aw)*SWRdir_wallsun*A_w/A_g + (1-aw)*SWRdir_wallshade*A_w/A_g
# Diffuse absorbed shortwave radiation
class SWRabsDiff_nT_Def():
pass
SWRabsDiff_nT = SWRabsDiff_nT_Def()
SWRabsDiff_nT.SWRabsGroundImp = (SWRabs_nT.SWRabsGroundImp-SWRabsDir_nT.SWRabsGroundImp)*Cimp
SWRabsDiff_nT.SWRabsGroundBare = (SWRabs_nT.SWRabsGroundBare-SWRabsDir_nT.SWRabsGroundBare)*Cbare
SWRabsDiff_nT.SWRabsGroundVeg = (SWRabs_nT.SWRabsGroundVeg-SWRabsDir_nT.SWRabsGroundVeg)*Cveg
SWRabsDiff_nT.SWRabsTree = 0
SWRabsDiff_nT.SWRabsWallSun = SWRabs_nT.SWRabsWallSun-SWRabsDir_nT.SWRabsWallSun
SWRabsDiff_nT.SWRabsWallShade = SWRabs_nT.SWRabsWallShade-SWRabsDir_nT.SWRabsWallShade
SWRabsDiff_nT.SWRabsTotalGround = SWRabs_nT.SWRabsTotalGround-SWRabsDir_nT.SWRabsTotalGround
SWRabsDiff_nT.SWRabsTotalCanyon = SWRabs_nT.SWRabsTotalCanyon-SWRabsDir_nT.SWRabsTotalCanyon
# Energy Balance of shortwave radiation
class SWREB_nT_Def():
pass
SWREB_nT = SWREB_nT_Def()
SWREB_nT.SWREBGroundImp = SWRin_nT.SWRinGroundImp - SWRout_nT.SWRoutGroundImp - SWRabs_nT.SWRabsGroundImp
SWREB_nT.SWREBGroundBare = SWRin_nT.SWRinGroundBare - SWRout_nT.SWRoutGroundBare - SWRabs_nT.SWRabsGroundBare
SWREB_nT.SWREBGroundVeg = SWRin_nT.SWRinGroundVeg - SWRout_nT.SWRoutGroundVeg - SWRabs_nT.SWRabsGroundVeg
SWREB_nT.SWREBTree = 0
SWREB_nT.SWREBWallSun = SWRin_nT.SWRinWallSun-SWRout_nT.SWRoutWallSun - SWRabs_nT.SWRabsWallSun
SWREB_nT.SWREBWallShade = SWRin_nT.SWRinWallShade-SWRout_nT.SWRoutWallShade - SWRabs_nT.SWRabsWallShade
SWREB_nT.SWREBTotalGround = SWRin_nT.SWRinTotalGround-SWRout_nT.SWRoutTotalGround - SWRabs_nT.SWRabsTotalGround
SWREB_nT.SWREBTotalCanyon = SWRin_nT.SWRinTotalCanyon-SWRout_nT.SWRoutTotalCanyon - SWRabs_nT.SWRabsTotalCanyon
if abs(SWREB_nT.SWREBGroundImp) >= 10**(-6):
print('SWREB_nT.SWREBGroundImp is not 0')
elif abs(SWREB_nT.SWREBGroundBare) >= 10**(-6):
print('SWREB_nT.SWREBGroundBare is not 0')
elif abs(SWREB_nT.SWREBGroundVeg) >= 10**(-6):
print('SWREB_nT.SWREBGroundVeg is not 0')
elif abs(SWREB_nT.SWREBWallSun) >= 10**(-6):
print('SWREB_nT.SWREBWallSun is not 0')
elif abs(SWREB_nT.SWREBWallShade) >= 10**(-6):
print('SWREB_nT.SWREBWallShade is not 0')
elif abs(SWREB_nT.SWREBTotalGround) >= 10**(-6):
print('SWREB_nT.SWREBTotalGround is not ')
elif abs(SWREB_nT.SWREBTotalCanyon) >= 10**(-6):
print('SWREB_nT.SWREBTotalCanyon is not 0')
return SWRin_nT,SWRout_nT,SWRabs_nT,SWRabsDir_nT,SWRabsDiff_nT,SWREB_nT
def SWRabsorbedWithTrees(self,h_can,w_can,h_tree,r_tree,d_tree,fgveg,fgbare,fgimp,aw,agveg,agbare,agimp,at,LAIt,
SWR_dir,SWR_diff,theta_Z,theta_n,ViewFactor,ParVegTree):
F_gs_T = ViewFactor.F_gs_T
F_gt_T = ViewFactor.F_gt_T
F_gw_T = ViewFactor.F_gw_T
F_ww_T = ViewFactor.F_ww_T
F_wt_T = ViewFactor.F_wt_T
F_wg_T = ViewFactor.F_wg_T
F_ws_T = ViewFactor.F_ws_T
F_sg_T = ViewFactor.F_sg_T
F_sw_T = ViewFactor.F_sw_T
F_st_T = ViewFactor.F_st_T
F_tg_T = ViewFactor.F_tg_T
F_tw_T = ViewFactor.F_tw_T
F_ts_T = ViewFactor.F_ts_T
F_tt_T = ViewFactor.F_tt_T
# normalized surface areas
A_s = w_can
A_g = w_can
A_w = h_can
# There are 2 trees. Hence, the area of tree is twice a circle
A_t = 2 * 2 * numpy.pi * r_tree
# load shortwave radiation
SWRdir_ground, SWRdir_wallsun, SWRdir_wallshade, SWRdir_tree = self.DirectSWRSurfaces(h_can, w_can, d_tree, h_tree,
r_tree, theta_Z, theta_n,
SWR_dir, LAIt, 1,ParVegTree)
# Balance direct shortwave radiation in
EB_SWRdir = SWR_dir-(SWRdir_ground*A_g/A_g+SWRdir_wallsun*A_w/A_g+SWRdir_wallshade*A_w/A_g+SWRdir_tree*A_t/A_g)
EB_SWRdiff = SWR_diff-(F_sg_T*SWR_diff+F_sw_T*SWR_diff+F_sw_T*SWR_diff+F_st_T*SWR_diff)
if abs(EB_SWRdir) >= 10 ** (-1):
print('EB_SWRdir is not 0',EB_SWRdir)
if abs(EB_SWRdiff) >= 10 ** (-1):
print('EB_SWRdiff is not 0',EB_SWRdiff)
SVF = numpy.zeros(4)
SVF[0] = F_gs_T+F_gt_T+2*F_gw_T
SVF[1] = F_ww_T+F_wt_T+F_wg_T+F_ws_T
SVF[2] = F_sg_T+2*F_sw_T+F_st_T
SVF[3] = F_ts_T+2*F_tw_T+F_tt_T+F_tg_T
SVF2 = numpy.zeros(4)
SVF2[0] = F_gs_T+2*F_ws_T*A_w+F_ts_T*A_t
SVF2[1] = F_sg_T+2*F_wg_T*A_w+F_tg_T*A_t
SVF2[2] = F_ww_T+F_sw_T*A_g/A_w+F_gw_T*A_g/A_w+F_tw_T*A_t/A_w
SVF2[3] = F_gt_T*A_g/A_t+2*F_wt_T*A_w/A_t+F_tt_T
for i in range(0, len(SVF)):
if SVF[i] < 0.999 or SVF[i] > 1.001:
print('The view factors do not add up to 1 for a canyon with trees')
if fgimp > 0:
Cimp = 1
else:
Cimp = 0
if fgbare > 0:
Cbare = 1
else:
Cbare = 0
if fgveg > 0:
Cveg = 1
else:
Cveg = 0
ai = [agveg, agbare, agimp, aw, aw, at, 0]
# View factor matrix to solve for infinite reflections equation
Tij = numpy.array([[1,0,0, -agveg*F_gw_T*Cveg, -agveg*F_gw_T*Cveg, -agveg*F_gt_T*Cveg, -agveg*F_gs_T*Cveg],
[0,1,0, -agbare*F_gw_T*Cbare, -agbare*F_gw_T*Cbare, -agbare*F_gt_T*Cbare, -agbare*F_gs_T*Cbare],
[0,0,1, -agimp*F_gw_T*Cimp, -agimp*F_gw_T*Cimp, -agimp*F_gt_T*Cimp, -agimp*F_gs_T*Cimp],
[-aw*F_wg_T*fgveg*Cveg,-aw*F_wg_T*fgbare*Cbare,-aw*F_wg_T*fgimp*Cimp, 1, -aw*F_ww_T, -aw*F_wt_T, -aw*F_ws_T],
[-aw*F_wg_T*fgveg*Cveg,-aw*F_wg_T*fgbare*Cbare,-aw*F_wg_T*fgimp*Cimp, -aw*F_ww_T, 1, -aw*F_wt_T, -aw*F_ws_T],
[-at*F_tg_T*fgveg*Cveg,-at*F_tg_T*fgbare*Cbare,-at*F_tg_T*fgimp*Cimp, -at*F_tw_T, -at*F_tw_T, 1-at*F_tt_T, -at*F_ts_T],
[0, 0, 0, 0, 0, 0, 1]])
# Incoming shortwave radiation from sky
Omega_i = numpy.array([agveg * SWRdir_ground * Cveg,
agbare * SWRdir_ground * Cbare,
agimp * SWRdir_ground * Cimp,
aw * SWRdir_wallsun,
aw * 0,
at * SWRdir_tree,
SWR_diff])
# Outgoing radiation per surface
B_i = numpy.linalg.solve(Tij,Omega_i)
if B_i[6] != SWR_diff:
print('Incoming lonwave radiation and emitted longwave radiation from the sky after the matrix inversion are not equal')
# Incoming shortwave radiation at each surface A_i
Tij2 = numpy.array([[0, 0, 0, F_gw_T*Cveg, F_gw_T*Cveg, F_gt_T*Cveg, F_gs_T*Cveg],
[0, 0, 0, F_gw_T*Cbare, F_gw_T*Cbare, F_gt_T*Cbare, F_gs_T*Cbare],
[0, 0, 0, F_gw_T*Cimp, F_gw_T*Cimp, F_gt_T*Cimp, F_gs_T*Cimp],
[F_wg_T*fgveg*Cveg, F_wg_T*fgbare*Cbare, F_wg_T*fgimp*Cimp, 0, F_ww_T, F_wt_T, F_ws_T],
[F_wg_T*fgveg*Cveg, F_wg_T*fgbare*Cbare, F_wg_T*fgimp*Cimp, F_ww_T, 0, F_wt_T, F_ws_T],
[F_tg_T*fgveg*Cveg, F_tg_T*fgbare*Cbare, F_tg_T*fgimp*Cimp, F_tw_T, F_tw_T, F_tt_T, F_ts_T],
[0, 0, 0, 0, 0, 0, 0]])
SWRdir_i = numpy.array([SWRdir_ground * Cveg,
SWRdir_ground * Cbare,
SWRdir_ground * Cimp,