diff --git a/.gitignore b/.gitignore
index f93affb683f92fce850d1bb29ec7cf273b915141..37f53ca0a07b155c41630ad29365790285e2be0a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,13 +4,17 @@
*.png
+# PyCharm project folder
.idea/
*.csv
-
+# STICS related files
*.sti
*param_files/*
*default/*
*logs/*
-*STICS/*
\ No newline at end of file
+*STICS/*
+
+# Profiling outputs
+*.prof
\ No newline at end of file
diff --git a/INPUTS/CIE_standard_skies.csv b/INPUTS/CIE_standard_skies.csv
new file mode 100644
index 0000000000000000000000000000000000000000..463d374b6483b463cf1d29d2263babd34d19bd03
--- /dev/null
+++ b/INPUTS/CIE_standard_skies.csv
@@ -0,0 +1,16 @@
+Type,Gradation,Indikatrix,a,b,c,d,e,Description of luminance distribution
+1,1,1,4.0 ,-0.70 ,0,-1.0 ,0.00 ,"CIE Standard Overcast Sky, alternative form. Steep luminance gradation towards zenith, azimuthal uniformity"
+2,1,2,4.0 ,-0.70 ,2,-1.5 ,0.15 ,"Overcast, with steep luminance gradation and slight brightening towards the sun"
+3,2,1,1.1 ,-0.8 ,0,-1.0 ,0.00 ,"Overcast, moderately graded with azimuthal uniformity"
+4,2,2,1.1 ,-0.8 ,2,-1.5 ,0.15 ,"Overcast, moderately graded and slight brightening towards the sun"
+5,3,1,0.0 ,-1.0 ,0,-1.0 ,0.00 ,"Sky of uniform luminance"
+6 ,3,2,0.0 ,-1.0 ,2,-1.5 ,0.15 ,"Partly cloudy sky, no gradation towards zenith, slight brightening towards the sun"
+7,3,3,0.0 ,-1.0 ,5,-2.5 ,0.30 ,"Partly cloudy sky, no gradation towards zenith, brighter circumsolar region"
+8,3,4,0.0 ,-1.0 ,10,-3.0 ,0.45 ,"Partly cloudy sky, no gradation towards zenith, distinct solar corona"
+9,4,2,-1.0 ,-0.55 ,2,-1.5 ,0.15 ,"Partly cloudy, with the obscured sun"
+10,4,3,-1.0 ,-0.55 ,5,-2.5 ,0.30 ,"Partly cloudy, with brighter circumsolar region"
+11,4,4,-1.0 ,-0.55 ,10,-3.0 ,0.45 ,"White-blue sky with distinct solar corona"
+12,5,4,-1.0 ,-0.32 ,10,-3.0 ,0.45 ,"CIE Standard Clear Sky, low illuminance turbidity"
+13 ,5,5,-1.0 ,-0.32 ,16,-3.0 ,0.30 ,"CIE Standard Clear Sky, polluted atmosphere"
+14 ,6,5,-1.0 ,-0.15 ,16,-3.0 ,0.30 ,"Cloudless turbid sky with broad solar corona"
+15,6,6,-1.0 ,-0.15 ,24,-2.8 ,0.15 ,"White-blue turbid sky with broad solar corona"
diff --git a/INPUTS/Igawa-5_sky_types_lut.csv b/INPUTS/Igawa-5_sky_types_lut.csv
new file mode 100644
index 0000000000000000000000000000000000000000..7a0f0c925a8d6fed9546f35e22c8dbc72e1afd71
--- /dev/null
+++ b/INPUTS/Igawa-5_sky_types_lut.csv
@@ -0,0 +1,599 @@
+;Kc;Cle;Sky category;CIE Sky Type
+0;0.01;0.0;O;1
+1;0.05;0.0;O;1
+2;0.1;0.0;O;1
+3;0.15;0.0;O;1
+4;0.2;0.0;O;1
+5;0.25;0.0;O;1
+6;0.3;0.0;O;1
+7;0.35;0.0;O;1
+8;0.4;0.0;O;1
+9;0.45;0.0;IO;4
+10;0.5;0.0;IO;4
+11;0.55;0.0;IO;4
+12;0.6;0.0;IO;4
+13;0.65;0.0;IO;4
+14;0.7;0.0;IO;4
+15;0.75;0.0;;5
+16;0.8;0.0;;5
+17;0.85;0.0;;5
+18;0.9;0.0;;5
+19;0.95;0.0;;5
+20;1.0;0.0;;5
+21;1.05;0.0;;5
+22;1.1;0.0;;5
+23;1.15;0.0;;5
+24;1.2;0.0;;5
+25;1.25;0.0;;5
+26;0.01;0.05;O;1
+27;0.05;0.05;O;1
+28;0.1;0.05;O;1
+29;0.15;0.05;O;1
+30;0.2;0.05;IO;4
+31;0.25;0.05;IO;4
+32;0.3;0.05;IO;4
+33;0.35;0.05;IO;4
+34;0.4;0.05;IO;4
+35;0.45;0.05;IO;4
+36;0.5;0.05;IO;4
+37;0.55;0.05;I;7
+38;0.6;0.05;I;7
+39;0.65;0.05;I;7
+40;0.7;0.05;I;7
+41;0.75;0.05;I;7
+42;0.8;0.05;I;7
+43;0.85;0.05;I;7
+44;0.9;0.05;;5
+45;0.95;0.05;;5
+46;1.0;0.05;;5
+47;1.05;0.05;;5
+48;1.1;0.05;;5
+49;1.15;0.05;;5
+50;1.2;0.05;;5
+51;1.25;0.05;;5
+52;0.01;0.1;O;1
+53;0.05;0.1;O;1
+54;0.1;0.1;IO;4
+55;0.15;0.1;IO;4
+56;0.2;0.1;IO;4
+57;0.25;0.1;IO;4
+58;0.3;0.1;IO;4
+59;0.35;0.1;IO;4
+60;0.4;0.1;IO;4
+61;0.45;0.1;I;7
+62;0.5;0.1;I;7
+63;0.55;0.1;I;7
+64;0.6;0.1;I;7
+65;0.65;0.1;I;7
+66;0.7;0.1;I;7
+67;0.75;0.1;I;7
+68;0.8;0.1;I;7
+69;0.85;0.1;I;7
+70;0.9;0.1;I;7
+71;0.95;0.1;;5
+72;1.0;0.1;;5
+73;1.05;0.1;;5
+74;1.1;0.1;;5
+75;1.15;0.1;;5
+76;1.2;0.1;;5
+77;1.25;0.1;;5
+78;0.01;0.15;O;1
+79;0.05;0.15;IO;4
+80;0.1;0.15;IO;4
+81;0.15;0.15;IO;4
+82;0.2;0.15;IO;4
+83;0.25;0.15;IO;4
+84;0.3;0.15;IO;4
+85;0.35;0.15;I;7
+86;0.4;0.15;I;7
+87;0.45;0.15;I;7
+88;0.5;0.15;I;7
+89;0.55;0.15;I;7
+90;0.6;0.15;I;7
+91;0.65;0.15;I;7
+92;0.7;0.15;I;7
+93;0.75;0.15;I;7
+94;0.8;0.15;I;7
+95;0.85;0.15;I;7
+96;0.9;0.15;I;7
+97;0.95;0.15;;5
+98;1.0;0.15;;5
+99;1.05;0.15;;5
+100;1.1;0.15;;5
+101;1.15;0.15;;5
+102;1.2;0.15;;5
+103;1.25;0.15;;5
+104;0.01;0.2;IO;4
+105;0.05;0.2;IO;4
+106;0.1;0.2;IO;4
+107;0.15;0.2;IO;4
+108;0.2;0.2;IO;4
+109;0.25;0.2;IO;4
+110;0.3;0.2;I;7
+111;0.35;0.2;I;7
+112;0.4;0.2;I;7
+113;0.45;0.2;I;7
+114;0.5;0.2;I;7
+115;0.55;0.2;I;7
+116;0.6;0.2;I;7
+117;0.65;0.2;I;7
+118;0.7;0.2;I;7
+119;0.75;0.2;I;7
+120;0.8;0.2;I;7
+121;0.85;0.2;I;7
+122;0.9;0.2;I;7
+123;0.95;0.2;I;7
+124;1.0;0.2;;5
+125;1.05;0.2;;5
+126;1.1;0.2;;5
+127;1.15;0.2;;5
+128;1.2;0.2;;5
+129;1.25;0.2;;5
+130;0.01;0.25;;5
+131;0.05;0.25;;5
+132;0.1;0.25;;5
+133;0.15;0.25;;5
+134;0.2;0.25;;5
+135;0.25;0.25;IO;4
+136;0.3;0.25;I;7
+137;0.35;0.25;I;7
+138;0.4;0.25;I;7
+139;0.45;0.25;I;7
+140;0.5;0.25;I;7
+141;0.55;0.25;I;7
+142;0.6;0.25;I;7
+143;0.65;0.25;I;7
+144;0.7;0.25;I;7
+145;0.75;0.25;I;7
+146;0.8;0.25;I;7
+147;0.85;0.25;I;7
+148;0.9;0.25;I;7
+149;0.95;0.25;I;7
+150;1.0;0.25;I;7
+151;1.05;0.25;I;7
+152;1.1;0.25;I;7
+153;1.15;0.25;;5
+154;1.2;0.25;;5
+155;1.25;0.25;;5
+156;0.01;0.3;;5
+157;0.05;0.3;;5
+158;0.1;0.3;;5
+159;0.15;0.3;;5
+160;0.2;0.3;;5
+161;0.25;0.3;I;7
+162;0.3;0.3;I;7
+163;0.35;0.3;I;7
+164;0.4;0.3;I;7
+165;0.45;0.3;I;7
+166;0.5;0.3;I;7
+167;0.55;0.3;I;7
+168;0.6;0.3;I;7
+169;0.65;0.3;I;7
+170;0.7;0.3;I;7
+171;0.75;0.3;I;7
+172;0.8;0.3;I;7
+173;0.85;0.3;I;7
+174;0.9;0.3;I;7
+175;0.95;0.3;I;7
+176;1.0;0.3;I;7
+177;1.05;0.3;I;7
+178;1.1;0.3;I;7
+179;1.15;0.3;I;7
+180;1.2;0.3;;5
+181;1.25;0.3;;5
+182;0.01;0.35;;5
+183;0.05;0.35;;5
+184;0.1;0.35;;5
+185;0.15;0.35;;5
+186;0.2;0.35;;5
+187;0.25;0.35;;5
+188;0.3;0.35;I;7
+189;0.35;0.35;I;7
+190;0.4;0.35;I;7
+191;0.45;0.35;I;7
+192;0.5;0.35;I;7
+193;0.55;0.35;I;7
+194;0.6;0.35;I;7
+195;0.65;0.35;I;7
+196;0.7;0.35;I;7
+197;0.75;0.35;I;7
+198;0.8;0.35;I;7
+199;0.85;0.35;I;7
+200;0.9;0.35;I;7
+201;0.95;0.35;I;7
+202;1.0;0.35;I;7
+203;1.05;0.35;I;7
+204;1.1;0.35;I;7
+205;1.15;0.35;I;7
+206;1.2;0.35;;5
+207;1.25;0.35;;5
+208;0.01;0.4;;5
+209;0.05;0.4;;5
+210;0.1;0.4;;5
+211;0.15;0.4;;5
+212;0.2;0.4;;5
+213;0.25;0.4;;5
+214;0.3;0.4;I;7
+215;0.35;0.4;I;7
+216;0.4;0.4;I;7
+217;0.45;0.4;I;7
+218;0.5;0.4;I;7
+219;0.55;0.4;I;7
+220;0.6;0.4;I;7
+221;0.65;0.4;I;7
+222;0.7;0.4;I;7
+223;0.75;0.4;I;7
+224;0.8;0.4;I;7
+225;0.85;0.4;I;7
+226;0.9;0.4;I;7
+227;0.95;0.4;I;7
+228;1.0;0.4;I;7
+229;1.05;0.4;I;7
+230;1.1;0.4;I;7
+231;1.15;0.4;I;7
+232;1.2;0.4;I;7
+233;1.25;0.4;;5
+234;0.01;0.45;;5
+235;0.05;0.45;;5
+236;0.1;0.45;;5
+237;0.15;0.45;;5
+238;0.2;0.45;;5
+239;0.25;0.45;;5
+240;0.3;0.45;;5
+241;0.35;0.45;I;7
+242;0.4;0.45;I;7
+243;0.45;0.45;I;7
+244;0.5;0.45;I;7
+245;0.55;0.45;I;7
+246;0.6;0.45;I;7
+247;0.65;0.45;I;7
+248;0.7;0.45;I;7
+249;0.75;0.45;I;7
+250;0.8;0.45;I;7
+251;0.85;0.45;I;7
+252;0.9;0.45;I;7
+253;0.95;0.45;I;7
+254;1.0;0.45;I;7
+255;1.05;0.45;I;7
+256;1.1;0.45;I;7
+257;1.15;0.45;I;7
+258;1.2;0.45;I;7
+259;1.25;0.45;I;7
+260;0.01;0.5;;5
+261;0.05;0.5;;5
+262;0.1;0.5;;5
+263;0.15;0.5;;5
+264;0.2;0.5;;5
+265;0.25;0.5;;5
+266;0.3;0.5;;5
+267;0.35;0.5;;5
+268;0.4;0.5;;5
+269;0.45;0.5;I;7
+270;0.5;0.5;I;7
+271;0.55;0.5;I;7
+272;0.6;0.5;I;7
+273;0.65;0.5;I;7
+274;0.7;0.5;I;7
+275;0.75;0.5;I;7
+276;0.8;0.5;I;7
+277;0.85;0.5;I;7
+278;0.9;0.5;I;7
+279;0.95;0.5;IC;11
+280;1.0;0.5;IC;11
+281;1.05;0.5;IC;11
+282;1.1;0.5;I;7
+283;1.15;0.5;I;7
+284;1.2;0.5;I;7
+285;1.25;0.5;I;7
+286;0.01;0.55;;5
+287;0.05;0.55;;5
+288;0.1;0.55;;5
+289;0.15;0.55;;5
+290;0.2;0.55;;5
+291;0.25;0.55;;5
+292;0.3;0.55;;5
+293;0.35;0.55;;5
+294;0.4;0.55;;5
+295;0.45;0.55;;5
+296;0.5;0.55;I;7
+297;0.55;0.55;I;7
+298;0.6;0.55;I;7
+299;0.65;0.55;I;7
+300;0.7;0.55;I;7
+301;0.75;0.55;I;7
+302;0.8;0.55;I;7
+303;0.85;0.55;IC;11
+304;0.9;0.55;IC;11
+305;0.95;0.55;IC;11
+306;1.0;0.55;IC;11
+307;1.05;0.55;IC;11
+308;1.1;0.55;IC;11
+309;1.15;0.55;IC;11
+310;1.2;0.55;I;7
+311;1.25;0.55;I;7
+312;0.01;0.6;;5
+313;0.05;0.6;;5
+314;0.1;0.6;;5
+315;0.15;0.6;;5
+316;0.2;0.6;;5
+317;0.25;0.6;;5
+318;0.3;0.6;;5
+319;0.35;0.6;;5
+320;0.4;0.6;;5
+321;0.45;0.6;;5
+322;0.5;0.6;I;7
+323;0.55;0.6;I;7
+324;0.6;0.6;I;7
+325;0.65;0.6;I;7
+326;0.7;0.6;I;7
+327;0.75;0.6;I;7
+328;0.8;0.6;I;7
+329;0.85;0.6;IC;11
+330;0.9;0.6;IC;11
+331;0.95;0.6;IC;11
+332;1.0;0.6;IC;11
+333;1.05;0.6;IC;11
+334;1.1;0.6;IC;11
+335;1.15;0.6;IC;11
+336;1.2;0.6;I;7
+337;1.25;0.6;I;7
+338;0.01;0.65;;5
+339;0.05;0.65;;5
+340;0.1;0.65;;5
+341;0.15;0.65;;5
+342;0.2;0.65;;5
+343;0.25;0.65;;5
+344;0.3;0.65;;5
+345;0.35;0.65;;5
+346;0.4;0.65;;5
+347;0.45;0.65;;5
+348;0.5;0.65;;5
+349;0.55;0.65;I;7
+350;0.6;0.65;I;7
+351;0.65;0.65;I;7
+352;0.7;0.65;I;7
+353;0.75;0.65;I;7
+354;0.8;0.65;IC;11
+355;0.85;0.65;IC;11
+356;0.9;0.65;IC;11
+357;0.95;0.65;IC;11
+358;1.0;0.65;IC;11
+359;1.05;0.65;IC;11
+360;1.1;0.65;IC;11
+361;1.15;0.65;IC;11
+362;1.2;0.65;IC;11
+363;1.25;0.65;I;7
+364;0.01;0.7;;5
+365;0.05;0.7;;5
+366;0.1;0.7;;5
+367;0.15;0.7;;5
+368;0.2;0.7;;5
+369;0.25;0.7;;5
+370;0.3;0.7;;5
+371;0.35;0.7;;5
+372;0.4;0.7;;5
+373;0.45;0.7;;5
+374;0.5;0.7;;5
+375;0.55;0.7;I;7
+376;0.6;0.7;I;7
+377;0.65;0.7;I;7
+378;0.7;0.7;I;7
+379;0.75;0.7;IC;11
+380;0.8;0.7;IC;11
+381;0.85;0.7;IC;11
+382;0.9;0.7;IC;11
+383;0.95;0.7;IC;11
+384;1.0;0.7;IC;11
+385;1.05;0.7;IC;11
+386;1.1;0.7;IC;11
+387;1.15;0.7;IC;11
+388;1.2;0.7;IC;11
+389;1.25;0.7;IC;11
+390;0.01;0.75;;5
+391;0.05;0.75;;5
+392;0.1;0.75;;5
+393;0.15;0.75;;5
+394;0.2;0.75;;5
+395;0.25;0.75;;5
+396;0.3;0.75;;5
+397;0.35;0.75;;5
+398;0.4;0.75;;5
+399;0.45;0.75;;5
+400;0.5;0.75;;5
+401;0.55;0.75;;5
+402;0.6;0.75;I;7
+403;0.65;0.75;I;7
+404;0.7;0.75;I;7
+405;0.75;0.75;IC;11
+406;0.8;0.75;IC;11
+407;0.85;0.75;IC;11
+408;0.9;0.75;IC;11
+409;0.95;0.75;C;13
+410;1.0;0.75;C;13
+411;1.05;0.75;C;13
+412;1.1;0.75;IC;11
+413;1.15;0.75;IC;11
+414;1.2;0.75;IC;11
+415;1.25;0.75;IC;11
+416;0.01;0.8;;5
+417;0.05;0.8;;5
+418;0.1;0.8;;5
+419;0.15;0.8;;5
+420;0.2;0.8;;5
+421;0.25;0.8;;5
+422;0.3;0.8;;5
+423;0.35;0.8;;5
+424;0.4;0.8;;5
+425;0.45;0.8;;5
+426;0.5;0.8;;5
+427;0.55;0.8;;5
+428;0.6;0.8;I;7
+429;0.65;0.8;I;7
+430;0.7;0.8;I;7
+431;0.75;0.8;IC;11
+432;0.8;0.8;IC;11
+433;0.85;0.8;IC;11
+434;0.9;0.8;C;13
+435;0.95;0.8;C;13
+436;1.0;0.8;C;13
+437;1.05;0.8;C;13
+438;1.1;0.8;C;13
+439;1.15;0.8;IC;11
+440;1.2;0.8;IC;11
+441;1.25;0.8;IC;11
+442;0.01;0.85;;5
+443;0.05;0.85;;5
+444;0.1;0.85;;5
+445;0.15;0.85;;5
+446;0.2;0.85;;5
+447;0.25;0.85;;5
+448;0.3;0.85;;5
+449;0.35;0.85;;5
+450;0.4;0.85;;5
+451;0.45;0.85;;5
+452;0.5;0.85;;5
+453;0.55;0.85;;5
+454;0.6;0.85;;5
+455;0.65;0.85;I;7
+456;0.7;0.85;I;7
+457;0.75;0.85;IC;11
+458;0.8;0.85;IC;11
+459;0.85;0.85;IC;11
+460;0.9;0.85;C;13
+461;0.95;0.85;C;13
+462;1.0;0.85;C;13
+463;1.05;0.85;C;13
+464;1.1;0.85;C;13
+465;1.15;0.85;IC;11
+466;1.2;0.85;IC;11
+467;1.25;0.85;IC;11
+468;0.01;0.9;;5
+469;0.05;0.9;;5
+470;0.1;0.9;;5
+471;0.15;0.9;;5
+472;0.2;0.9;;5
+473;0.25;0.9;;5
+474;0.3;0.9;;5
+475;0.35;0.9;;5
+476;0.4;0.9;;5
+477;0.45;0.9;;5
+478;0.5;0.9;;5
+479;0.55;0.9;;5
+480;0.6;0.9;;5
+481;0.65;0.9;;5
+482;0.7;0.9;I;7
+483;0.75;0.9;IC;11
+484;0.8;0.9;IC;11
+485;0.85;0.9;IC;11
+486;0.9;0.9;C;13
+487;0.95;0.9;C;13
+488;1.0;0.9;C;13
+489;1.05;0.9;C;13
+490;1.1;0.9;C;13
+491;1.15;0.9;IC;11
+492;1.2;0.9;IC;11
+493;1.25;0.9;IC;11
+494;0.01;0.95;;5
+495;0.05;0.95;;5
+496;0.1;0.95;;5
+497;0.15;0.95;;5
+498;0.2;0.95;;5
+499;0.25;0.95;;5
+500;0.3;0.95;;5
+501;0.35;0.95;;5
+502;0.4;0.95;;5
+503;0.45;0.95;;5
+504;0.5;0.95;;5
+505;0.55;0.95;;5
+506;0.6;0.95;;5
+507;0.65;0.95;;5
+508;0.7;0.95;I;7
+509;0.75;0.95;IC;11
+510;0.8;0.95;IC;11
+511;0.85;0.95;IC;11
+512;0.9;0.95;C;13
+513;0.95;0.95;C;13
+514;1.0;0.95;C;13
+515;1.05;0.95;C;13
+516;1.1;0.95;C;13
+517;1.15;0.95;IC;11
+518;1.2;0.95;IC;11
+519;1.25;0.95;;5
+520;0.01;1.0;;5
+521;0.05;1.0;;5
+522;0.1;1.0;;5
+523;0.15;1.0;;5
+524;0.2;1.0;;5
+525;0.25;1.0;;5
+526;0.3;1.0;;5
+527;0.35;1.0;;5
+528;0.4;1.0;;5
+529;0.45;1.0;;5
+530;0.5;1.0;;5
+531;0.55;1.0;;5
+532;0.6;1.0;;5
+533;0.65;1.0;;5
+534;0.7;1.0;;5
+535;0.75;1.0;IC;11
+536;0.8;1.0;IC;11
+537;0.85;1.0;IC;11
+538;0.9;1.0;C;13
+539;0.95;1.0;C;13
+540;1.0;1.0;C;13
+541;1.05;1.0;C;13
+542;1.1;1.0;C;13
+543;1.15;1.0;IC;11
+544;1.2;1.0;;5
+545;1.25;1.0;;5
+546;0.01;1.05;;5
+547;0.05;1.05;;5
+548;0.1;1.05;;5
+549;0.15;1.05;;5
+550;0.2;1.05;;5
+551;0.25;1.05;;5
+552;0.3;1.05;;5
+553;0.35;1.05;;5
+554;0.4;1.05;;5
+555;0.45;1.05;;5
+556;0.5;1.05;;5
+557;0.55;1.05;;5
+558;0.6;1.05;;5
+559;0.65;1.05;;5
+560;0.7;1.05;;5
+561;0.75;1.05;IC;11
+562;0.8;1.05;IC;11
+563;0.85;1.05;IC;11
+564;0.9;1.05;C;13
+565;0.95;1.05;C;13
+566;1.0;1.05;C;13
+567;1.05;1.05;C;13
+568;1.1;1.05;C;13
+569;1.15;1.05;IC;11
+570;1.2;1.05;;5
+571;1.25;1.05;;5
+572;0.01;1.1;;5
+573;0.05;1.1;;5
+574;0.1;1.1;;5
+575;0.15;1.1;;5
+576;0.2;1.1;;5
+577;0.25;1.1;;5
+578;0.3;1.1;;5
+579;0.35;1.1;;5
+580;0.4;1.1;;5
+581;0.45;1.1;;5
+582;0.5;1.1;;5
+583;0.55;1.1;;5
+584;0.6;1.1;;5
+585;0.65;1.1;;5
+586;0.7;1.1;;5
+587;0.75;1.1;;5
+588;0.8;1.1;IC;11
+589;0.85;1.1;IC;11
+590;0.9;1.1;C;13
+591;0.95;1.1;C;13
+592;1.0;1.1;C;13
+593;1.05;1.1;C;13
+594;1.1;1.1;C;13
+595;1.15;1.1;IC;11
+596;1.2;1.1;;5
+597;1.25;1.1;;5
diff --git a/INPUTS/SCENARIOS/Example1_loc.yaml b/INPUTS/SCENARIOS/Example1_loc.yaml
index 78d08b575d1005336638862985a2d5b4a325b26b..dee079f35d43a2de97c5f7a663fb4bf617849d03 100644
--- a/INPUTS/SCENARIOS/Example1_loc.yaml
+++ b/INPUTS/SCENARIOS/Example1_loc.yaml
@@ -71,6 +71,14 @@ FibonacciSamples:
Type: integer
Limit: [1,100000]
Definition: Samples in the Fibonacci sky subdivision scheme.
+DiffuseSkyType:
+ Value: From weather data
+ Type: string
+ Possibilities: [From weather data, Uniform]
+ Definition: Determine how to compute sky diffuse radiance distribution.
+ "From weather data" uses a simplified implementation of Igawa (2014)
+ to derive the sky type from weather data.
+ "Uniform" assumes uniform radiance distribution (faster but less accurate).
Xmin_InterestZone:
Value: -2
Type: float
diff --git a/MODULES/ENVIRONMENT/light.py b/MODULES/ENVIRONMENT/light.py
index c09c7a213c1799f72662ced061200ba47867f0d1..207fdde09ab1dd10acb1d8314c6fa838b116fa29 100644
--- a/MODULES/ENVIRONMENT/light.py
+++ b/MODULES/ENVIRONMENT/light.py
@@ -5,38 +5,23 @@
#Authors : Roxane Bruhwyler (roxane.bruhwyler@uliege.be or roxane.bruhwyler@hotmail.com) and Nicolas De Cock (nicolas.decock1@gmail.com)
#This file is part of the PASE software, and is distributed under the MIT license.
+from calendar import isleap
import logging
import numpy as np
import pandas as pd
import pyvista as pyV
import pvlib.solarposition as pvlibSP
+from pvlib.atmosphere import get_relative_airmass
+from pvlib.irradiance import get_extra_radiation
import matplotlib.pyplot as plt
import os
import MODULES.conversion_functions as cf
-from MODULES.ENVIRONMENT.sky_model import ReinhartSky
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky, fibonacci_half_sphere
+from MODULES.ENVIRONMENT.sky_model import CIEStandardSky
logger = logging.getLogger(__name__)
-def fibonacci_half_sphere(samples=18):
- """
- Function computing a number of direction sampling a virtual upper half sphere with a center at 0,0,0
-
- Parameters:
- samples (int): Number of directions
-
- Returns:
- Directions (np.array): Matrix (samples x 3) providing the direction
- """
- phi = np.pi * (3. - np.sqrt(5.))
- i = np.linspace(0,samples-1,num=samples)
- yp = (1 - i/float(samples-1))
- radius = np.sqrt(1-yp**2)
- theta = phi * i
- xp = np.cos(theta) * radius
- zp = np.sin(theta) * radius
- return np.column_stack([xp,zp,yp])
-
class Sun_positions:
@@ -103,7 +88,12 @@ class Sun_positions:
return solar_vector
def get_top_of_atm_radiation(self, index, n):
-
+ """
+
+ :param index: datetime index
+ :param n: seems unused ?
+ :return: irradiance at top of atmosphere on a horizontal surface [W/m²]
+ """
day_of_year = np.array(index.dayofyear)
n_days_in_year = day_of_year[len(index)-1]
solar_declination = self.get_solar_declination(day_of_year)
@@ -137,21 +127,37 @@ class Light:
def __init__(self, WD, SP, ghi_multiplier=1):
self.data = {}
+ sky_type_lut_path = os.path.join('INPUTS', 'Igawa-5_sky_types_lut.csv')
+ self.sky_type_lut = pd.read_csv(sky_type_lut_path, sep=';')
for year in WD.keys():
GHI = WD[year]['G(h)'].to_numpy() * ghi_multiplier
- if int(year)%4 == 0:
+ if isleap(int(year)):
rad_top_atm = SP.sp_leapY['Top_atm_radiation'].to_numpy()
+ apparent_sun_zenith = SP.sp_leapY['apparent_zenith'].to_numpy()
+ sun_elevation = SP.sp_leapY['elevation'].to_numpy()
else:
rad_top_atm = SP.sp_nonleapY['Top_atm_radiation'].to_numpy()
-
+ apparent_sun_zenith = SP.sp_nonleapY['apparent_zenith'].to_numpy()
+ sun_elevation = SP.sp_nonleapY['elevation'].to_numpy()
+
kt = self.get_clearness_sky_index(rad_top_atm, GHI)
DHI = self.get_diffuse_horizontal_radiation(kt, GHI)
BHI = self.get_beam_horizontal_radiation(GHI, DHI)
Ai = self.get_anisotropy_index(rad_top_atm, BHI)
f = self.get_modulating_factor(GHI, BHI)
+
+ # Extraterrestrial Normal Irradiance (used for the Kc and Cle below)
+ ENI = get_extra_radiation(WD[year]['DateTime'].dt.dayofyear.values)
+
+ # Meteorological indices Clear Sky Index (Kc) and Cloudless index
+ # (Cle) from Igawa (2014)
+ Kc = self.get_clear_sky_index(ENI, GHI, apparent_sun_zenith)
+ Cle = self.get_cloudless_index(DHI, GHI, sun_elevation)
+
+ cie_sky_type = self.get_sky_type(Kc, Cle)
df = pd.DataFrame({'GHI': GHI.tolist(),
'BHI': BHI.tolist(),
@@ -159,9 +165,12 @@ class Light:
'rad_top_atm': rad_top_atm.tolist(),
'kt': kt.tolist(),
'Ai': Ai.tolist(),
- 'f': f.tolist()}, index=WD[year].index)
+ 'f': f.tolist(),
+ 'Kc': Kc.tolist(),
+ 'Cle': Cle.tolist(),
+ 'CIE Sky Type': cie_sky_type}, index=WD[year].index)
- self.data[year] = df
+ self.data[year] = df
def get_clearness_sky_index(self, rad_top_atm, GHI):
@@ -214,8 +223,63 @@ class Light:
f[GHI>0] = np.sqrt(BHI[GHI>0]/GHI[GHI>0])
return f
-
-
+
+ def get_clear_sky_index(self, ENI, GHI, apparent_sun_zenith):
+
+ # Get relative airmass from pvlib method
+ # (takes apparent zenith angle in [deg] for Kasten and Young 1989 model)
+ m = get_relative_airmass(apparent_sun_zenith, model='kastenyoung1989')
+
+ GHI_non_zero = GHI != 0
+
+ Kc = np.empty(len(GHI))
+ Kc[:] = np.nan
+
+ Kc = np.divide(GHI, 0.84 * ENI / m * np.exp(-0.054 * m), out=Kc, where=GHI_non_zero)
+
+ return Kc
+
+ def get_cloudless_index(self, DHI, GHI, sun_elevation):
+
+ sun_elevation_rad = np.deg2rad(sun_elevation)
+ GHI_non_zero = GHI != 0
+
+ # Cloud ratio [-]
+ Ce = self.get_cloud_ratio(DHI, GHI, GHI_non_zero)
+
+ # Standard cloud ratio [-]
+ Ces = (0.08302
+ + 0.5358 * np.exp(-17.3 * sun_elevation_rad)
+ + 0.3818 * np.exp(-3.2899 * sun_elevation_rad))
+
+ Cle = np.zeros(len(GHI))
+ Cle = np.divide((1 - Ce), (1 - Ces), out=Cle, where=GHI_non_zero)
+
+ return Cle
+
+ def get_cloud_ratio(self, DHI, GHI, GHI_non_zero):
+
+ Ce = np.zeros((len(GHI),))
+
+ Ce = np.divide(DHI, GHI, out=Ce, where=GHI_non_zero)
+
+ return Ce
+
+ def get_sky_type(self, Kc_target, Cle_target):
+
+ mydf = pd.DataFrame({'Kc': Kc_target, 'Cle': Cle_target})
+ temp_list = []
+ for k, v in mydf.iterrows():
+ if np.isnan(v['Kc']):
+ temp_list.append(np.nan)
+ else:
+ i = ((self.sky_type_lut['Kc']-v['Kc']) *
+ (self.sky_type_lut['Cle']-v['Cle'])).abs().idxmin()
+ temp_list.append(int(self.sky_type_lut['CIE Sky Type'].iloc[i]))
+
+ return temp_list
+
+
class Sun_positions_sampled:
def __init__(self, lat, long, precision_lvl, loc_name, freq_deter, TZ):
@@ -371,17 +435,19 @@ class Ray_casting_scene:
The class is initiate with a mesh containing the source points and a geometry
'''
#The class light shade scene init with a geometry (pyvista.polydata) and a mesh instance
- def __init__(self,mesh,geometry):
+ def __init__(self, mesh, geometry, discrete_sky):
self.mesh = mesh
self.sourcepoints = self.mesh.get_sourcepoints()
+
self.geometry = geometry
self.n_sourcepoints = self.sourcepoints.shape[0]
self.sources_flag_dict = mesh.get_sources_flag_dict()
-
-
- def get_light_maps(self, sun_P, scheme='Reinhart', MF=1, n_small_suns=180):
+
+ self.discrete_sky = discrete_sky
+
+ def get_light_maps(self, sun_P):
- if type(self.geometry) == list:
+ if type(self.geometry) == list: # sun tracking
sv = np.zeros((1,3))
self.dir_map = np.zeros((len(self.sourcepoints),len(sun_P[:,0])))
self.diff_map = np.zeros((len(self.sourcepoints),len(sun_P[:,0])))
@@ -389,19 +455,13 @@ class Ray_casting_scene:
for time in range(len(sun_P[:,0])):
print(time)
geometry = self.geometry[time]
- difff_map = self.diffuse_map(geometry,
- scheme=scheme,
- MF=MF,
- n_small_suns=n_small_suns)
+ difff_map = self.diffuse_map(geometry)
sv[0,:] = sun_P[time,:]
dirrr_map = self.direct_map(sv, geometry)
self.dir_map[:,time] = dirrr_map
self.diff_map[:,time] = difff_map
- else:
- self.diff_map = self.diffuse_map(self.geometry,
- scheme=scheme,
- MF=MF,
- n_small_suns=n_small_suns)
+ else: # no sun tracking
+ self.diff_map = self.diffuse_map(self.geometry)
self.dir_map = self.direct_map(sun_P, self.geometry)
def self_intercept(self,SourcePoints,intercept_points,id_rays_stopped,tol = 0.01):
@@ -424,45 +484,33 @@ class Ray_casting_scene:
return np.unique(id_rays_stopped[delta>tol])
- def diffuse_map(self, geometry, scheme='Reinhart', MF=1, n_small_suns=180):
+ def diffuse_map(self, geometry):
"""
- Public method, compute the diffuse light at the point sources defined in the input mesh using
- the approximation of a isotropic half sphere sky.
- The method uses the mesh and the geometry set at the initialization of the instance
-
- Parameters:
- n_small_suns (int): number of sources consider in the sky for the diffuse light computation
- higher number will provide a better accuracy but heavier computation
+ Public method, compute the diffuse light at the point sources defined
+ in the input mesh (see this class' constructor) under a
+ discrete_sky sky model.
+
+ :param geometry: geometry set at the initialization of the instance
+ :return: diffuse_mask: diffuse map binary mask for each source point and discrete_sky element. Shape: (n_sky_patches, n_source_points)
- Returns:
- Diffu (np.array 1 x n): Providing a vector with the fraction ([0-1]) of diffuse light
- for each of the "n" source points defined in the mesh
"""
-
-
+
#Get direction of ray to reach the small suns and compute the sky view of each point
- if scheme.lower() == 'reinhart':
- sky = ReinhartSky(MF=MF)
- pTarget = np.column_stack([sky.reinhart_patches.x,
- sky.reinhart_patches.y,
- sky.reinhart_patches.z])
- n_small_suns = len(sky.reinhart_patches)
- elif scheme.lower() == 'fibonacci':
- pTarget = fibonacci_half_sphere(n_small_suns)
- else:
- raise NotImplementedError('Unrecognized sky discretization scheme')
+ pTarget = np.column_stack([self.discrete_sky.x,
+ self.discrete_sky.y,
+ self.discrete_sky.z])
+ n_sky_elements = len(self.discrete_sky)
-
- #Creation of the source points array (Nx3) with N = len(Source) * len(n_small_suns)
+ #Creation of the source points array (Nx3) with N = len(Source) * len(n_sky_elements)
SourcePoints = np.repeat(np.column_stack((
self.sourcepoints[:,0],
self.sourcepoints[:,1],
self.sourcepoints[:,2]
)),
- n_small_suns,
+ n_sky_elements,
axis=0)
- #Creation of the target points array (Nx3) with N = len(Source) * len(n_small_suns)
+ #Creation of the target points array (Nx3) with N = len(Source) * len(n_sky_elements)
TargetPoints = np.tile(pTarget,[self.n_sourcepoints,1])
#Computation of the ray interception of the N rays
@@ -480,7 +528,7 @@ class Ray_casting_scene:
SourceID = np.repeat(np.linspace(0,
self.n_sourcepoints-1,
self.n_sourcepoints),
- n_small_suns,
+ n_sky_elements,
axis=0)
#Touched provide a list with the sourceID of the intercept ray
@@ -491,24 +539,17 @@ class Ray_casting_scene:
# Compute the cos(zenith angle) of all the small suns for the normalization
_, _zenith_angle_all = cf.get_zenith_angle_from_cart(TargetPoints)
- _cos_zenith_angle_all = np.cos(_zenith_angle_all).reshape(self.n_sourcepoints, n_small_suns)
+ _cos_zenith_angle_all = np.cos(_zenith_angle_all).reshape(self.n_sourcepoints, n_sky_elements)
# Compute the cos(zenith angle) of the small suns that do NOT contribute to the diffuse map
# (i.e. rays that were intercepted)
- _mask = np.ones(_cos_zenith_angle_all.size, bool)
- _mask[id_rays_stopped_filtred] = 0
- _cos_zenith_angle_blocked = _cos_zenith_angle_all.copy()
- _mask = _mask.reshape(self.n_sourcepoints, n_small_suns)
- _cos_zenith_angle_blocked[_mask] = 0
-
- #Creation of the empty matrix of sky view
- Diffu = np.ones(self.n_sourcepoints, dtype=np.float16)
+ diffuse_mask = np.ones(_cos_zenith_angle_all.size, bool)
+ diffuse_mask[id_rays_stopped_filtred] = 0
- #Computation of the sky view by removing the fraction of intercepted ray at each location
- Diffu[unique.astype("int")] = 1 - np.sum(_cos_zenith_angle_blocked[unique.astype("int"), :], axis=1)/np.sum(_cos_zenith_angle_all[unique.astype("int"), :], axis=1)
+ diffuse_mask = diffuse_mask.reshape(self.n_sourcepoints, n_sky_elements)
- return Diffu
+ return diffuse_mask
def get_direct_map_by_flag(self,Flags):
"""
@@ -641,9 +682,9 @@ class Ray_casting_scene:
dfShade = SP_sampled.tz_localize(None).reset_index()
dfShade['doy'] = dfShade['index'].dt.dayofyear
dfShade['RefDate'] = pd.to_datetime(dfShade['index'].dt.date)
- # dfShade['hour'] = dfShade['index'].dt.hour
+ # dfShade['hour'] = dfShade['index'].dt.hour
dfShade['second'] = pd.to_timedelta(dfShade['index'].dt.time.astype(str)).dt.total_seconds()
- # dfShade['RefHour'] = dfShade['hour']
+ # dfShade['RefHour'] = dfShade['hour']
dfShade['RefSecond'] = dfShade['second']
for year in light_data:
@@ -657,10 +698,8 @@ class Ray_casting_scene:
n = 6
- #dfWeather = light_data[year].tz_localize('Etc/GMT+0').reset_index()
dfWeather = light_data[year].tz_localize(None).reset_index()
dfWeather['doy'] = dfWeather['index'].dt.dayofyear
- # dfWeather['hour'] = dfWeather['index'].dt.hour
dfWeather['second'] = pd.to_timedelta(dfWeather['index'].dt.time.astype(str)).dt.total_seconds()
# Jointure entre les données météos et les données d'ombrage en vue de déterminée la date/index liée aux données d'ombrage
@@ -669,34 +708,58 @@ class Ray_casting_scene:
irradianceMap_direct = {}
irradianceMap_diffus = {}
- # Boucle sur les n jours de l'annee afin de calculer l'irradiation journaliere
+ # Loop over days of year (doy) to compute daily irradiance.
doy = dfWeatherMerged['index'].dt.dayofyear.unique()
- doy.sort()
+ doy.sort() # doy = [1, 2, 3, ..., 365]
for day in doy:
- #Creation de sous dataframe comprenant les donnees du jour numero "doy"
+ #Creation de sous dataframe comprenant les donnees du jour numero "day"
sublight_df = dfWeatherMerged.loc[dfWeatherMerged['index'].dt.dayofyear==day,:]
df_subShade = dfShade.loc[sublight_df['RefDate'].unique()[0]==dfShade['RefDate']]
df_subShade_merged = df_subShade.join(sublight_df[['second','BHI','GHI','DHI']].set_index('second'),on='second',how='left',rsuffix='_',lsuffix="__")
#df_subShade_merged = df_subShade.join(sublight_df[['hour','BHI','GHI','DHI']].set_index('hour'),on='hour',how='left',rsuffix='_',lsuffix="__")
#Jointure sur l'instant de la journee la plus proche sur base des secondes écoulées depuis le debut de la journee
- df_subShade_merged = pd.merge_asof(df_subShade,sublight_df[['second','BHI','GHI','DHI','doy']].set_index('second'),on=['second'],direction='nearest',suffixes=('_x','_y')).sort_values('second')
+ df_subShade_merged = pd.merge_asof(df_subShade,sublight_df[['second','BHI','GHI','DHI','doy', 'CIE Sky Type']].set_index('second'),on=['second'],direction='nearest',suffixes=('_x','_y')).sort_values('second')
#Calcul de l'irradiation directe et diffuse dans ce jour et injection de la donnee dans le dictionnaire lie
irradianceMap_direct[day] = np.sum(self.dir_map[:,list(df_subShade.index)] * df_subShade_merged['BHI'].to_numpy(),axis=1)*10**-6*60*60/n
- if type(self.geometry) != list:
- irradianceMap_diffus[day] = np.sum(df_subShade_merged['DHI'].to_numpy())*self.diff_map*10**-6*60*60/n
- else:
- print('day')
+ if type(self.geometry) != list: # no sun tracking
+ irradianceMap_diffus[day] = self.compute_daily_diff_irradiation(df_subShade_merged.dropna())
+ else: # sun tracking -> in this case, the mask changes at each time step
+ print(f'{day=}')
irradianceMap_diffus[day] = np.sum(df_subShade_merged['DHI'].to_numpy()*self.diff_map[:,list(df_subShade.index)], axis=1)*10**-6*60*60/n
-
+
#Conversion des dictionnaires en matrice numpy et ajout dans l attribut ad-hoc
- self.daily_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_diffus).to_numpy() + pd.DataFrame.from_dict(irradianceMap_direct).to_numpy()
- self.daily_dir_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_direct).to_numpy()
+ self.daily_irr_spat[year] = (pd.DataFrame.from_dict(irradianceMap_diffus).to_numpy()
+ + pd.DataFrame.from_dict(irradianceMap_direct).to_numpy())
+ self.daily_dir_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_direct).to_numpy()
self.daily_diff_irr_spat[year] = pd.DataFrame.from_dict(irradianceMap_diffus).to_numpy()
-
+ def compute_daily_diff_irradiation(self, df):
+ # loop over each timestep in the day's dataframes
+ for i in range(len(df)):
+ unweighted_radiance = CIEStandardSky(self.discrete_sky,
+ df['azimuth'].iloc[i],
+ df['elevation'].iloc[i],
+ df['CIE Sky Type'].iloc[i]).rel_radiance_distribution
+
+ # sky patch area and cos(z) weighting
+ weighted_radiance = unweighted_radiance * self.discrete_sky['Normalized surf area'] * self.discrete_sky['cos(z)']
+
+ # normalize to relative radiance contributions
+ rel_radiance_contrib = weighted_radiance/np.sum(weighted_radiance)
+
+ # multiply by the mask
+ if type(self.geometry) == list: # sun tracking
+ raise NotImplementedError('Sun tracking is temporarily broken')
+ # shaded_radiance_contrib = rel_radiance_contrib * self.diff_map[:,time]
+ else:
+ shaded_radiance_contrib = rel_radiance_contrib.values * self.diff_map
+
+ diff_irradiance_map = (shaded_radiance_contrib * df['DHI'].iloc[i]).sum(axis=1)
+
+ return diff_irradiance_map
def show_light_map(light_matrix, msh_grid, PV_central, lim_min, lim_max, lgd_title):
diff --git a/MODULES/ENVIRONMENT/sky_model.py b/MODULES/ENVIRONMENT/sky_model.py
index 04da300cc5f3e20fafcc8e316462dcc9d0ec535b..ac9d75f3ba1a0587ae958a601ae82dca12dc05eb 100644
--- a/MODULES/ENVIRONMENT/sky_model.py
+++ b/MODULES/ENVIRONMENT/sky_model.py
@@ -5,10 +5,13 @@ Paper openly available on Researchgate: https://www.google.com/url?sa=t&source=w
"""
+import cProfile
import logging
import math
import numpy as np
+import os
import pandas as pd
+from pathlib import Path
from matplotlib import pyplot as plt
import matplotlib
matplotlib.use('TkAgg')
@@ -17,6 +20,27 @@ from MODULES.conversion_functions import sph_to_cart
logger = logging.getLogger(__name__)
+
+def fibonacci_half_sphere(samples=18):
+ """
+ Function computing a number of direction sampling a virtual upper half sphere with a center at 0,0,0
+
+ Parameters:
+ samples (int): Number of directions
+
+ Returns:
+ Directions (np.array): Matrix (samples x 3) providing the direction
+ """
+ phi = np.pi * (3. - np.sqrt(5.))
+ i = np.linspace(0, samples - 1, num=samples)
+ yp = (1 - i / float(samples - 1))
+ radius = np.sqrt(1 - yp ** 2)
+ theta = phi * i
+ xp = np.cos(theta) * radius
+ zp = np.sin(theta) * radius
+ return np.column_stack([xp, zp, yp])
+
+
class ReinhartSky:
def __init__(self, MF=1):
@@ -38,6 +62,9 @@ class ReinhartSky:
# Compute x, y, z coordinates of unit vectors pointing towards the sky patches
self.get_patches_xyz()
+ # Compute cos(zenith angle), used in irradiance computations
+ self.compute_cos_zenith()
+
def get_original_tregenza(self):
original_tregenza_segments = pd.DataFrame()
@@ -112,12 +139,14 @@ class ReinhartSky:
return reinhart_sky, reinhart_num_total
def build_patches(self):
- self.reinhart_patches = pd.DataFrame()
self.patch_id = 0
self.cols = ['row', 'patch_id', 'patch_id_in_row', 'd_el', 'd_az', 'el',
'az', 'solid_angle_sr']
+
+ # temporary list that will be appended to (dictionaries) in the loop,
+ # then converted to Dataframe after the loop
+ temp_list = []
for row in self.reinhart_sky['row']:
- # print(f'{row=}')
d_az = 360 / \
self.reinhart_sky[self.reinhart_sky['row'] == row][
'num_segments'].values[
@@ -125,21 +154,22 @@ class ReinhartSky:
el = self.alpha * (row + 1 / 2) # [deg]
for segment in range(self.reinhart_sky['num_segments'].iloc[row]):
- # print(f'{segment=}')
az = d_az * (segment + 1 / 2) # [deg]
solid_angle = self.compute_solid_angle_from_elev_azim(el,
- d_az)
+ d_az) # [sr]
- df = pd.DataFrame(
- [[row, self.patch_id, segment, self.alpha, d_az, el, az,
- solid_angle]],
- columns=self.cols)
+ # Dictionary that will be appended to temp_list
+ dict = {'row': row, 'patch_id': self.patch_id, 'patch_id_in_row': segment,
+ 'd_el': self.alpha, 'd_az': d_az, 'el': el,
+ 'az': az, 'solid_angle_sr': solid_angle}
- self.reinhart_patches = pd.concat([self.reinhart_patches, df],
- ignore_index=True)
+ temp_list.append(dict)
self.patch_id += 1
+ # Convert the list of dictionaries to DataFrame
+ self.reinhart_patches = pd.DataFrame.from_records(temp_list)
+
def add_top_patch(self):
top_row = max(self.reinhart_sky['row']) + 1
top_patch_solid_angle = self.compute_solid_angle_cone(self.alpha / 2)
@@ -191,9 +221,9 @@ class ReinhartSky:
def compute_surface_areas(self):
# Compute surf area of each patch
- self.reinhart_patches['surf area'] = self.reinhart_patches['solid_angle_sr'] / (
+ self.reinhart_patches['Normalized surf area'] = self.reinhart_patches['solid_angle_sr'] / (
2 * np.pi)
- sum_surf_areas = self.reinhart_patches['surf area'].sum()
+ sum_surf_areas = self.reinhart_patches['Normalized surf area'].sum()
logger.debug(f'Sum of surface areas = {sum_surf_areas:.4g}')
@@ -202,9 +232,9 @@ class ReinhartSky:
sum_solid_angles = self.reinhart_patches['solid_angle_sr'].sum()
sum_solid_angles_div_pi = sum_solid_angles / np.pi
logger.debug(f'Sum of solid angles = {sum_solid_angles_div_pi:.4g} * pi')
- if not math.isclose(sum_solid_angles, 2*np.pi):
+ if not math.isclose(sum_solid_angles, 2*np.pi, rel_tol=5e-3):
logger.warning('Sum of solid angles seems off, check the sky discretization scheme')
- logger.warning(f'Sum of solid angles = {sum_solid_angles_div_pi:.4g} * pi')
+ logger.warning(f'Sum of solid angles = {sum_solid_angles_div_pi:.6g} * pi')
def get_patches_xyz(self):
# Compute xyz coords on a unit sphere (actually half sphere since it's the sky dome)
@@ -214,6 +244,9 @@ class ReinhartSky:
azimut=self.reinhart_patches['az'],
elev=self.reinhart_patches['el'])
+ def compute_cos_zenith(self):
+ self.reinhart_patches['cos(z)'] = np.cos(np.deg2rad(90-self.reinhart_patches['el']))
+
def test_MFs(self):
MFs = np.arange(1, 33)
for MF in MFs:
@@ -262,9 +295,215 @@ class ReinhartSky:
plt.show()
+class CIEStandardSky:
+ """
+ Generate a standard radiance distribution among the 15 CIE General Skies.
+ """
+
+ def __init__(self, discrete_sky,
+ sun_az, sun_el,
+ sky_type=5):
+
+ self.sky = discrete_sky
+
+ self.az = discrete_sky['az'].values
+ self.el = discrete_sky['el'].values
+
+ self.az_s = sun_az # [°] Sun azimuth (145° is south-east)
+ self.el_s = sun_el # [°] Sun elevation
+
+ self.sky_type = int(sky_type)
+
+ self.rel_radiance_distribution = self.compute_rel_radiance()
+
+
+ def compute_sph_angular_dist(self, Z, Z_s, az, az_s):
+ """
+
+ :param Z: Sky patch zenith angle [deg]
+ :param Z: float
+ :param Z_s: Solar zenith angle [deg]
+ :param Z_s: float
+ :param az: Sky patch azimuth angle from north [deg]
+ :param az: float
+ :param az_s: Solar azimuth angle from north [deg]
+ :param az_s: float
+ :return: chi = Spherical angular distance [rad]
+ """
+
+ # Convert all angles from degrees to radians
+ Z = np.deg2rad(Z)
+ Z_s = np.deg2rad(Z_s)
+ az = np.deg2rad(az)
+ az_s = np.deg2rad(az_s)
+
+ # Absolute difference in azimuth angles
+ Az = np.abs(az - az_s)
+
+ # Angular distance chi
+ chi = np.arccos(
+ np.cos(Z_s) * np.cos(Z) + np.sin(Z_s) * np.sin(Z) * np.cos(Az))
+
+ return chi
+
+ def gradation_fun(self, a, b, Z):
+ """
+ Gradation function. It is not used as is but vectorized in compute_gradation().
+
+ :param a: Horizon-zenith gradient (gradation parameter) [-]
+ :type a: float
+ :param b: Gradient intensity (gradation parameter) [-]
+ :type b: float
+ :param Z: Zenith angle [deg]
+ :type Z: float
+ :return: gradation function value [-]
+ """
+ if math.isclose(Z, np.pi / 2):
+ return 1
+ else:
+ return 1 + a * np.exp(b / np.cos(Z))
+
+ def compute_gradation(self, a, b, Z):
+ """
+ Vectorizes and computes the gradation function
+
+ :param a: Horizon-zenith gradient (gradation parameter) [-]
+ :type a: float
+ :param b: Gradient intensity (gradation parameter) [-]
+ :type b: float
+ :param Z: Zenith angle [deg]
+ :type Z: Nx1 array of float
+ :return: gradation function [-], same shape as Z
+ """
+
+ v_gradation = np.vectorize(self.gradation_fun)
+
+ Z = np.deg2rad(Z)
+
+ return v_gradation(a, b, Z)
+
+ def scattering_indicatrix_fun(self, c, d, e, x):
+ """
+ Compute the scattering indicatrix value.
+
+ :param c: Circumsolar intensity (scattering parameter) [-]
+ :type c: float
+ :param d: Circumsolar radius (scattering parameter) [-]
+ :type d: float
+ :param e: Backscattering effect (scattering parameter) [-]
+ :type e: float
+ :param x: angle (spherical angular distance or zenith, i.e. 0) [rad]
+ :type x: array of float or float
+ :return: Scattering indicatrix value [-], same shape as x
+ """
+ return 1 + c * (np.exp(d * x) - np.exp(d * np.pi / 2)) + e * (
+ np.cos(x)) ** 2
+
+ def get_sky_params(self, standard_sky):
+ """
+
+ :param standard_sky: Dataframe of length 1 with standard sky information
+ :return a: Horizon-zenith gradient (gradation parameter) [-]
+ :return b: Gradient intensity (gradation parameter) [-]
+ :return c: Circumsolar intensity (scattering parameter) [-]
+ :return d: Circumsolar radius (scattering parameter) [-]
+ :return e: Backscattering effect (scattering parameter) [-]
+ """
+
+ a = standard_sky['a'].values[0]
+ b = standard_sky['b'].values[0]
+ c = standard_sky['c'].values[0]
+ d = standard_sky['d'].values[0]
+ e = standard_sky['e'].values[0]
+
+ return a, b, c, d, e
+
+ def relative_radiance_fun(self, sky_type=None, az=None, el=None):
+ """
+
+ :return: relative quantity (radiance or luminance) with respect to
+ value at the zenith
+ """
+
+ # Read the csv containing the descriptions of the 15 standard skies
+
+ pase_dir_path = Path(__file__).parents[2]
+ standard_skies = pd.read_csv(os.path.join(pase_dir_path,
+ 'INPUTS',
+ 'CIE_standard_skies.csv'))
+
+ # Extract the parameters of the desired sky type
+ if sky_type is None:
+ standard_sky = standard_skies.loc[standard_skies['Type'] == self.sky_type]
+ else:
+ standard_sky = standard_skies.loc[standard_skies['Type'] == sky_type]
+
+ # Sky model parameters
+ a, b, c, d, e = self.get_sky_params(standard_sky)
+
+ # Elevation to zenith (Sun)
+ Z_s = 90 - self.el_s # [deg]
+
+ # Elevation to zenith (sky patch)
+ if el is None:
+ Z = 90 - self.el # [deg]
+ else:
+ Z = 90 - el # [deg]
+
+
+ # Spherical angular dist
+ if az is None:
+ chi = self.compute_sph_angular_dist(Z, Z_s, self.az, self.az_s)
+ else:
+ chi = self.compute_sph_angular_dist(Z, Z_s, az, self.az_s)
+
+ # Gradation
+ phi = self.compute_gradation(a, b, Z)
+ phi_zenith = self.compute_gradation(a, b, Z=0) # TODO do not recompute each time !
+
+ # Scattering indicatrix
+ indicatrix = self.scattering_indicatrix_fun(c, d, e, chi)
+ indicatrix_zenith = self.scattering_indicatrix_fun(c, d, e, np.deg2rad(Z_s)) # TODO do not recompute each time !
+
+ # Relative luminance/radiance (called "rel_quantity" for genericity)
+ rel_quantity = phi * indicatrix / (phi_zenith * indicatrix_zenith)
+
+ return rel_quantity
+
+ def compute_rel_radiance(self, sky_type=None, az=None, el=None):
+ v_radiance_fun = np.vectorize(self.relative_radiance_fun)
+ return v_radiance_fun(sky_type, az, el)
+
+
if __name__ == '__main__':
- sky = ReinhartSky(MF=1)
- sky.scatter_2D()
- sky.scatter_2D(text_id=True)
- sky.scatter_3D()
- sky.scatter_3D(text_id=True)
+
+ profile_flag = False
+
+ if profile_flag:
+ MFs = [1, 2]
+ for MF in MFs:
+ profile_output = os.path.join('..', '..', 'OUTPUTS', 'profiling', f'cprofile_output_sky_model_MF{MF}.prof')
+ cProfile.run(f"ReinhartSky(MF={MF})", profile_output, sort='tottime')
+ else:
+ MF = 1
+
+ sky = ReinhartSky(MF=MF)
+ if MF == 1:
+ text_id = True
+ else:
+ text_id = False
+
+ if MF < 8:
+ sky.scatter_2D(text_id=text_id)
+ sky.scatter_3D(text_id=text_id)
+
+ sky.reinhart_patches['Relative radiance distribution'] = CIEStandardSky(sky.reinhart_patches,
+ sun_az=147.67, sun_el=38.02, sky_type=15).rel_radiance_distribution
+
+ DHI = 100 # [W/m²] arbitrary value for example's sake
+ integral_rel_sky_radiance = (sky.reinhart_patches['Relative radiance distribution']
+ * sky.reinhart_patches['solid_angle_sr']).sum()
+ zenith_radiance = DHI/integral_rel_sky_radiance
+ sky.reinhart_patches['Absolute radiance distribution'] = (sky.reinhart_patches['Relative radiance distribution']
+ * zenith_radiance)
+
diff --git a/TEMPORARY_FILES/parse_cprofile_data.py b/TEMPORARY_FILES/parse_cprofile_data.py
new file mode 100644
index 0000000000000000000000000000000000000000..6ccc77a09d3b377aa4a5421d362b860edabd66a8
--- /dev/null
+++ b/TEMPORARY_FILES/parse_cprofile_data.py
@@ -0,0 +1,18 @@
+import os
+import pstats
+from pstats import SortKey
+
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky
+
+# Profiles of sky_model Reinhart module
+MF = 8
+MFs = [16, 32]
+for MF in MFs:
+ profile_output = os.path.join('..', 'OUTPUTS', 'profiling', f'cprofile_output_sky_model_MF{MF}.prof')
+ p = pstats.Stats(profile_output)
+ p.strip_dirs().sort_stats('tottime').print_stats(10)
+
+# Profile of example.py
+profile_output = os.path.join('..', 'OUTPUTS', 'profiling', f'example.prof')
+p = pstats.Stats(profile_output)
+p.strip_dirs().sort_stats('tottime').print_stats(100)
diff --git a/TEMPORARY_FILES/perez_standard_skies.py b/TEMPORARY_FILES/perez_standard_skies.py
new file mode 100644
index 0000000000000000000000000000000000000000..7937019b714ffdda15c0d696f37f98487da98085
--- /dev/null
+++ b/TEMPORARY_FILES/perez_standard_skies.py
@@ -0,0 +1,348 @@
+import logging
+import math
+import numpy as np
+import os
+import pandas as pd
+import pyvista as pv
+import matplotlib
+matplotlib.use('TkAgg')
+from matplotlib import pyplot as plt
+
+from MODULES.conversion_functions import sph_to_cart
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky
+
+logger = logging.getLogger(__name__)
+logger.setLevel('DEBUG')
+
+def compute_sph_angular_dist(Z, Z_s, az, az_s):
+ """
+
+ Args:
+ Z: Sky patch zenith angle [deg]
+ Z_s: Solar zenith angle [deg]
+ az: Sky patch azimuth angle from north [deg]
+ az_s: Solar azimuth angle from north [deg]
+
+ Returns:
+ chi: Spherical angular distance [rad]
+ """
+
+ # Convert all angles from degrees to radians
+ Z = np.deg2rad(Z)
+ Z_s = np.deg2rad(Z_s)
+ az = np.deg2rad(az)
+ az_s = np.deg2rad(az_s)
+
+ # Absolute difference in azimuth angles
+ Az = np.abs(az-az_s)
+
+ # Angular distance chi
+ chi = np.arccos(np.cos(Z_s) * np.cos(Z) + np.sin(Z_s)*np.sin(Z)*np.cos(Az))
+
+ return chi
+
+def gradation_fun(a, b, Z):
+ """
+ Gradation function defined to be vectorized in compute_gradation().
+ Args:
+ a: Horizon-zenith gradient (gradation parameter) [-]
+ b: Gradient intensity (gradation parameter) [-]
+ Z: zenith angle [rad]
+
+ Returns:
+ Gradation value.
+ """
+ if math.isclose(Z, np.pi/2):
+ return 1
+ else:
+ return 1 + a*np.exp(b/np.cos(Z))
+
+def compute_gradation(a, b, Z):
+ """
+ Computes the gradation function based on vectorized gradation_fun function.
+ Args:
+ a: Horizon-zenith gradient (gradation parameter) [-]
+ b: Gradient intensity (gradation parameter) [-]
+ Z: Zenith angle [deg]. Type : Nx1 array
+
+ Returns:
+ gradation function result
+ """
+ v_gradation = np.vectorize(gradation_fun)
+
+ Z = np.deg2rad(Z)
+
+ return v_gradation(a, b, Z)
+
+def scattering_indicatrix_fun(c, d, e, x):
+ """
+
+ Args:
+ c: Circumsolar intensity (scattering parameter) [-]
+ d: Circumsolar radius (scattering parameter) [-]
+ e: Backscattering effect (scattering parameter) [-]
+ x: angle (sphericla angular distance or zenith, i.e. 0) [rad]
+ Returns:
+ Scattering indicatrix value
+ """
+ return 1 + c * (np.exp(d*x) - np.exp(d*np.pi/2)) + e * (np.cos(x))**2
+
+def get_sky_params(standard_sky):
+ """
+
+ Args:
+ standard_sky: Dataframe of length 1 with standard sky information
+
+ Returns:
+ a: Horizon-zenith gradient (gradation parameter) [-]
+ b: Gradient intensity (gradation parameter) [-]
+ c: Circumsolar intensity (scattering parameter) [-]
+ d: Circumsolar radius (scattering parameter) [-]
+ e: Backscattering effect (scattering parameter) [-]
+ """
+
+ a = standard_sky['a'].values[0]
+ b = standard_sky['b'].values[0]
+ c = standard_sky['c'].values[0]
+ d = standard_sky['d'].values[0]
+ e = standard_sky['e'].values[0]
+
+ return a, b, c, d, e
+
+def luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el):
+ # Standard sky parameters and description
+ standard_sky = standard_skies.loc[standard_skies['Type'] == sky_type]
+
+ # Sky model parameters
+ a, b, c, d, e = get_sky_params(standard_sky)
+
+ # Elevation to zenith (sky patch)
+ Z = 90 - patch_el # [deg]
+
+ # Elevation to zenith (Sun)
+ Z_s = 90 - sun_el # [deg]
+
+ # Spherical angular dist
+ chi = compute_sph_angular_dist(Z, Z_s, patch_az, sun_az)
+
+ # Gradation
+ phi = compute_gradation(a, b, Z)
+ phi_zenith = compute_gradation(a, b, Z=0)
+
+ # Indicatrix
+ indicatrix = scattering_indicatrix_fun(c, d, e, chi)
+ indicatrix_zenith = scattering_indicatrix_fun(c, d, e, np.deg2rad(Z_s))
+
+ # Relative luminance
+ rel_lum = phi * indicatrix / (phi_zenith * indicatrix_zenith)
+
+ return rel_lum
+
+def compute_luminance(sky_type, patch_az, patch_el, sun_az, sun_el):
+ v_luminance_fun = np.vectorize(luminance_fun)
+ return v_luminance_fun(sky_type, patch_az, patch_el, sun_az, sun_el)
+
+def print_n_max_values(vect, n):
+ a = vect.copy()
+ b = np.unique(a)
+ c = np.sort(b)
+ c = c[::-1]
+ print(c[:n])
+
+def plot_luminance_map(az, el, lum, az_s=None, el_s=None, luminance_or_radiance='radiance', polar=True, colormap=None):
+ lum_min = np.min(lum)
+ lum_max = np.max(lum)
+
+ fig = plt.figure()
+ if polar:
+ ax = fig.add_subplot(111, polar=True)
+ c = ax.pcolormesh(np.deg2rad(az), 90-el, lum, vmin=lum_min, vmax=lum_max, cmap=colormap)
+ else:
+ ax = fig.gca()
+
+ c = ax.pcolormesh(az, el, lum, vmin=lum_min, vmax=lum_max, cmap=colormap)
+
+ ax.set_xlabel('Azimuth [deg]')
+ ax.set_ylabel('Zenith angle [deg]')
+
+ cbar = fig.colorbar(c, ax=ax)
+ if luminance_or_radiance == 'luminance':
+ cbar_label = 'Luminance [cd/m²]'
+ elif luminance_or_radiance == 'radiance':
+ cbar_label = 'Radiance [W/(m²⋅sr)]'
+ cbar.set_label(cbar_label)
+
+ if (az_s is not None) and (el_s is not None):
+ if polar:
+ ax.plot(np.deg2rad(az_s), 90-el_s, 'ro', label='Sun position')
+ else:
+ ax.plot(az_s, el_s, 'ro', label='Sun position')
+ ax.legend()
+
+ if polar:
+ # Set N at the top, E to the right
+ ax.set_theta_zero_location('N')
+ ax.set_theta_direction(-1)
+
+ plt.show()
+
+def plot_all_skies(polar=True, luminance_or_radiance='luminance', colormap='jet'):
+ # sky = ReinhartSky().reinhart_patches
+ # az = sky['az']
+ # el = sky['el']
+ az = np.arange(0, 360, step=5) # [°]
+ el = np.arange(0, 90, step=5) # [°]
+
+ az_meshgrid, el_meshgrid = np.meshgrid(az, el)
+
+ el_s = 38.02 # [°]
+ az_s = 147.67 # [°] 145° is south-east
+ zenith_luminance = 4404 # cd/m²
+
+ fig, axs = plt.subplots(3, 5,
+ figsize=(14, 8),
+ subplot_kw=dict(polar=polar))
+
+ sky_types = np.arange(1, 16)
+ for sky_type in sky_types:
+ lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s)
+
+ lum_min = np.min(lum)
+ lum_max = np.max(lum)
+
+ ax_row = int(np.floor((sky_type-1)/5))
+ ax_col = (sky_type-1) % 5
+ print(f'{sky_type=}', f'{ax_row=}', f'{ax_col=}')
+
+ ax = axs[ax_row, ax_col]
+ if polar:
+ c = ax.pcolormesh(np.deg2rad(az), 90-el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
+ else:
+ c = ax.pcolormesh(az, el, lum.reshape(az_meshgrid.shape), vmin=lum_min, vmax=lum_max, cmap=colormap)
+
+ ax.set_xlabel('Azimuth [deg]')
+ ax.set_ylabel('Zenith [deg]')
+
+ ax.title.set_text(f'Sky type {sky_type}')
+
+ # Set N at the top, E to the right
+ ax.set_theta_zero_location('N')
+ ax.set_theta_direction(-1)
+
+ fig.colorbar(c, ax=ax)
+
+ if (az_s is not None) and (el_s is not None):
+ if polar:
+ ax.plot(np.deg2rad(az_s), 90-el_s, 'ro', markeredgecolor='yellow', label='Sun position')
+ else:
+ ax.plot(az_s, 90-el_s, 'ro', markeredgecolor='yellow', label='Sun position')
+
+ if luminance_or_radiance == 'luminance':
+ fig_title = 'Luminance [cd/m²]'
+ elif luminance_or_radiance == 'radiance':
+ fig_title = 'Radiance [W/(m²⋅sr)]'
+ fig.suptitle(fig_title + ' (azimuth-zenith)')
+
+ plt.tight_layout()
+ plt.show()
+
+def plot_3D(sky, rel_lum):
+
+ sky['patch_bound_el_low'] = sky['el'] - sky['d_el'] / 2
+ sky['patch_bound_el_up'] = sky['el'] + sky['d_el'] / 2
+ sky['patch_bound_az_low'] = sky['az'] - sky['d_az'] / 2
+ sky['patch_bound_az_up'] = sky['az'] + sky['d_az'] / 2
+ print(sky['patch_bound_el_up'].iloc[-1])
+ # Fix elevation for top patch
+ sky['patch_bound_el_up'].iloc[-1] = 90
+ print(sky['patch_bound_el_up'].iloc[len(sky) - 1])
+
+ # Vertical levels
+ # in this case a single level slightly above the surface of a sphere
+ RADIUS = 1
+ levels = [RADIUS * 1.01]
+
+ xx_bounds = np.concatenate([sky['patch_bound_az_low'].values,
+ [sky['patch_bound_az_up'].values[
+ -1]]]) # az
+ xx_bounds[-2:-1] += 180
+ yy_bounds = np.concatenate([sky['patch_bound_el_low'].values,
+ [sky['patch_bound_el_up'].values[
+ -1]]]) # el
+
+ grid_scalar = pv.grid_from_sph_coords(xx_bounds, 90 - yy_bounds, levels)
+
+ # And fill its cell arrays with the scalar data
+ grid_scalar.cell_data["Luminance relative to zenith"] = np.array(
+ rel_lum.reshape(az_meshgrid.shape)).swapaxes(-2, -1).ravel("C")
+
+ # Prepare the Sun
+ xyz_sun = sph_to_cart('deg', azimut=az_s, elev=el_s, zenith_angle=None,
+ dist=RADIUS)
+ sun_mesh = pv.Sphere(radius=RADIUS / 20, center=xyz_sun)
+
+ # Make a plot
+ p = pv.Plotter()
+ p.add_mesh(sun_mesh, color='yellow')
+ p.add_mesh(grid_scalar, opacity=0.8, cmap="jet")
+ p.show()
+
+if __name__ == '__main__':
+
+ standard_skies = pd.read_csv(os.path.join('..', 'INPUTS', 'CIE_standard_skies.csv'))
+
+ # Example from the paper
+ sky_type = 12
+
+ sky_mesh = 'arange' # {'arange', 'linspace', 'Reinhart'}
+
+ if sky_mesh.lower() == 'reinhart':
+ sky = ReinhartSky(MF=1).reinhart_patches
+ az = sky['az'].values
+ el = sky['el'].values
+ # meshgrid_flag = False
+ meshgrid_flag = True
+ az_meshgrid, el_meshgrid = np.meshgrid(az, el)
+ else:
+ if sky_mesh.lower() == 'arange':
+ step = 5
+ az = np.arange(0, 360+step, step=step) # [°]
+ el = np.arange(0, 90+step, step=step) # [°]
+ elif sky_mesh.lower() == 'linspace':
+ az = np.linspace(0, 360, num=145) # [°]
+ el = np.linspace(0, 90, num=145) # [°]
+
+ else:
+ raise ValueError('Unrecognized sky discretization')
+
+ meshgrid_flag = True
+ az_meshgrid, el_meshgrid = np.meshgrid(az, el)
+
+ el_s = 38.02 # [°]
+ az_s = 147.67 # [°] 145° is south-east
+ zenith_luminance = 4404 # cd/m²
+
+ if meshgrid_flag:
+ rel_lum = compute_luminance(sky_type, az_meshgrid.flatten(), el_meshgrid.flatten(), az_s, el_s).reshape(az_meshgrid.shape)
+ else:
+ # rel_lum = compute_luminance(sky_type, az, el, az_s, el_s)
+ rel_lum = compute_luminance(sky_type, np.append(az, [360]), np.append(el, [90]), az_s, el_s)
+
+ lum = rel_lum * zenith_luminance
+
+ print_n_max_values(rel_lum, 3)
+
+ # Plot with pcolormesh
+ plot_luminance_map(az, el, rel_lum,
+ az_s, el_s,
+ luminance_or_radiance='radiance',
+ polar=True,
+ colormap='jet')
+
+ # Plot all 15 standard skies in a (5, 3) subplots figure (takes a couple of minutes)
+ # plot_all_skies()
+
+if sky_mesh.lower() == 'reinhart':
+ # Plot with PyVista (only for Reinhart scheme)
+ plot_3D(sky, rel_lum)
+
diff --git a/example.py b/example.py
index 5c384b1725690892176378cd53f1afbe7b8ab3bb..50796fb74d54f29ccfeddc8d2a05ba608838b3be 100644
--- a/example.py
+++ b/example.py
@@ -17,6 +17,7 @@ from MODULES.PHOTOVOLTAICS.configuration import PV_Configuration_3D
from MODULES.ENVIRONMENT.light import Sun_positions_sampled, Sun_positions, Light
from MODULES.ENVIRONMENT.light import show_light_map2, Ray_casting_scene
from MODULES.ENVIRONMENT.mesh import Mesh
+from MODULES.ENVIRONMENT.sky_model import ReinhartSky
from MODULES.PHOTOVOLTAICS.production import PV_Production
from MODULES.CROPS.run_crop_simulations import run_crop_simu
@@ -62,7 +63,7 @@ Sun_positions_complete = Sun_positions(Loc_1['Latitude'],
# Instantiation of the 3D PV central
PV_1_3Dconfig = PV_Configuration_3D(PV_params_dict,
Sun_positions_samp.solar_vector,
- visualization=True) # !!!! Problem with rotation angle that are negative
+ visualization=False) # !!!! Problem with rotation angle that are negative
# Initiation of the object containing points of interest to compute light
M = Mesh()
@@ -76,15 +77,18 @@ M.add_plane_ground_regular_meshes(Loc_1['Xmin_InterestZone'],
Loc_1['dY_InterestZone'],
flag="crop")
+# Discrete sky model
+discrete_sky = ReinhartSky(MF=Loc_1['MF']).reinhart_patches
+
# Computation of sun and light data
Light_instance = Light(WD.nyears_data, Sun_positions_complete)
# Iniation and run of light ray casting model (direct and diffuse) with points of interest and scene
-L = Ray_casting_scene(mesh=M, geometry=PV_1_3Dconfig.PV_central_PD)
-L.get_light_maps(Sun_positions_samp.solar_vector,
- scheme=Loc_1['SkyDiscretizationScheme'],
- MF=Loc_1['MF'],
- n_small_suns=Loc_1['FibonacciSamples'])
+L = Ray_casting_scene(mesh=M,
+ geometry=PV_1_3Dconfig.PV_central_PD,
+ discrete_sky=discrete_sky)
+
+L.get_light_maps(Sun_positions_samp.solar_vector)
# Integration of irradiation along days
L.get_daily_irradiation_map(Sun_positions_samp.SP,
@@ -100,8 +104,11 @@ show_light_map2(L.sourcepoints[:,:-1],
PV_1_3Dconfig.PV_central_PD,
"Direct map [-]")
+# Display the Sky visibility map (= fraction of sky patches viewed by each control point)
+# It is purely derived from the geometry of the central
+# (i.e. computed once, or for each tracking position if tracking is activated)
show_light_map2(L.sourcepoints[:,:-1],
- np.array(L.diff_map,dtype=np.float32),
+ np.array(L.diff_map.sum(axis=1)/len(discrete_sky),dtype=np.float32),
PV_1_3Dconfig.PV_central_PD,
"Sky visibility map [-]")