.....

7 files changed, 521 insertions, 49 deletions

diff --git a/DisplayClass.cpp b/DisplayClass.cpp
index 8d79e74..174e4ac 100644
--- a/DisplayClass.cpp
+++ b/DisplayClass.cpp
@@ -61,18 +61,19 @@ DisplayClass::~DisplayClass() {
 }
 
 void DisplayClass::renderScene(unsigned long delta) {
-  angle += delta*0.025;
+  angle += delta*0.055;
   if(angle >= 360)
     angle -= 360;
 
   glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
   
   glLoadIdentity(); // Clean up matrix
-  glTranslatef(0.0, 0.0, -10.0); // Then set up transformation
+  glTranslatef(0.0, -2.0, -10.0); // Then set up transformation
+  glRotatef(15, 1.0, 0.0, 0.0);
   glRotatef(angle, 0.0, 1.0, 0.0);
-  glRotatef(angle*2, 1.0, 0.0, 0.0);
-  glRotatef(angle*3, 0.0, 0.0, 1.0);
-  glRotatef(-angle*5, 1.0, 1.0, 1.0);
+//  glRotatef(angle*2, 1.0, 0.0, 0.0);
+//  glRotatef(angle*3, 0.0, 0.0, 1.0);
+//  glRotatef(-angle*5, 1.0, 1.0, 1.0);
   
   
   static Temparray temp(20, x, y, z);
diff --git a/House.cpp b/House.cpp
index 4b2d83c..82e114a 100644
--- a/House.cpp
+++ b/House.cpp
@@ -1,8 +1,391 @@
 #include "House.h"
 
-House::House (float inittemp0, int earthx0, int earthy0, int earthz0, float collectortemp0,
-              float latitude0):
-              collectortemp(collectortemp0), latitude(latitude0) {
+House::House(float inittemp0, int earthx0, int earthy0, int earthz0, float collectortemp0,
+             float radiation0, float latitude0, float eta00, float k10, float k20, float collectorarea0,
+             float watthoursperhuman0, int humans0):
+             collectortemp(collectortemp0), latitude(latitude0), radiation(radiation0), eta0(eta00),
+             k1(k10), k2(k20), watthoursperhuman(watthoursperhuman0), humans(humans0),
+             collectorarea(collectorarea0){
   Temparray temp0(inittemp0, earthx0, earthy0, earthz0);
   temp = temp0;
+  calcefficiency();
+  calctemp();
+}
+
+void House::calcefficiency() {
+  int ioffset;
+  float x[11], y[11], a, b, j;
+
+  for(int i=0; i<11; ++i) x[i] = i*10;
+
+  y[0] = 0.755;
+  y[1] = 0.695;
+  y[2] = 0.635;
+  y[3] = 0.575;
+  y[4] = 0.523;
+  y[5] = 0.468;
+  y[6] = 0.404;
+  y[7] = 0.34;
+  y[8] = 0.282;
+  y[9] = 0.222;
+  y[10] = 0.155;
+  
+  for(int i=0; i<100; ++i) {
+    if(i >= x[0] && i < x[1]) {
+      a = (y[1]-y[0])/10;
+      b = y[0];
+      ioffset = 0;
+    }
+    if(i >= x[1] && i < x[2]) {
+      a = (y[2]-y[1])/10;
+      b = y[1];
+      ioffset = 10;
+    }
+    if(i >= x[2] && i < x[3]) {
+      a = (y[3]-y[2])/10;
+      b = y[2];
+      ioffset = 20;
+    }
+    if(i >= x[3] && i < x[4]) {
+      a = (y[4]-y[3])/10;
+      b = y[3];
+      ioffset = 30;
+    }
+    if(i >= x[4] && i < x[5]) {
+      a = (y[5]-y[4])/10;
+      b = y[4];
+      ioffset = 40;
+    }
+    if(i >= x[5] && i < x[6]) {
+      a = (y[6]-y[5])/10;
+      b = y[5];
+      ioffset = 50;
+    }
+    if(i >= x[6] && i < x[7]) {
+      a = (y[7]-y[6])/10;
+      b = y[6];
+      ioffset = 60;
+    }
+    if(i >= x[7] && i < x[8]) {
+      a = (y[8]-y[7])/10;
+      b = y[7];
+      ioffset = 70;
+    }
+    if(i >= x[8] && i < x[9]) {
+      a = (y[9]-y[8])/10;
+      b = y[8];
+      ioffset = 80;
+    }
+    if(i >= x[9] && i < x[10]) {
+      a = (y[10]-y[9])/10;
+      b = y[9];
+      ioffset = 90;
+    }
+    efficiencyflat[i] = a*(i-ioffset)+b;
+  }
+
+  for(int i=0; i<11; ++i){y[i] = eta0 - (k1/radiation)*i*10 - (k2/radiation)*i*i*100;}
+
+  for(int i=0; i<100; ++i) {
+    if(i >= x[0] && i < x[1]) {
+      a = (y[1]-y[0])/10;
+      b = y[0];
+      ioffset = 0;
+    }
+    if(i >= x[1] && i < x[2]) {
+      a = (y[2]-y[1])/10;
+      b = y[1];
+      ioffset = 10;
+    }
+    if(i >= x[2] && i < x[3]) {
+      a = (y[3]-y[2])/10;
+      b = y[2];
+      ioffset = 20;
+    }
+    if(i >= x[3] && i < x[4]) {
+      a = (y[4]-y[3])/10;
+      b = y[3];
+      ioffset = 30;
+    }
+    if(i >= x[4] && i < x[5]) {
+      a = (y[5]-y[4])/10;
+      b = y[4];
+      ioffset = 40;
+    }
+    if(i >= x[5] && i < x[6]) {
+      a = (y[6]-y[5])/10;
+      b = y[5];
+      ioffset = 50;
+    }
+    if(i >= x[6] && i < x[7]) {
+      a = (y[7]-y[6])/10;
+      b = y[6];
+      ioffset = 60;
+    }
+    if(i >= x[7] && i < x[8]) {
+      a = (y[8]-y[7])/10;
+      b = y[7];
+      ioffset = 70;
+    }
+    if(i >= x[8] && i < x[9]) {
+      a = (y[9]-y[8])/10;
+      b = y[8];
+      ioffset = 80;
+    }
+    if(i >= x[9] && i < x[10]) {
+      a = (y[10]-y[9])/10;
+      b = y[9];
+      ioffset = 90;
+    }
+    efficiencyvacuum[i] = a*(i-ioffset)+b;
+  }
+  anglecorrectionfactortransversal[0] = 1;
+  anglecorrectionfactortransversal[1] = 0.998;
+  anglecorrectionfactortransversal[2] = 0.996;
+  anglecorrectionfactortransversal[3] = 0.994;
+  anglecorrectionfactortransversal[4] = 0.992;
+  anglecorrectionfactortransversal[5] = 0.99;
+  anglecorrectionfactortransversal[6] = 0.988;
+  anglecorrectionfactortransversal[7] = 0.986;
+  anglecorrectionfactortransversal[8] = 0.984;
+  anglecorrectionfactortransversal[9] = 0.982;
+  anglecorrectionfactortransversal[10] = 0.98;
+  anglecorrectionfactortransversal[11] = 0.98;
+  anglecorrectionfactortransversal[12] = 0.98;
+  anglecorrectionfactortransversal[13] = 0.98;
+  anglecorrectionfactortransversal[14] = 0.98;
+  anglecorrectionfactortransversal[15] = 0.98;
+  anglecorrectionfactortransversal[16] = 0.982;
+  anglecorrectionfactortransversal[17] = 0.984;
+  anglecorrectionfactortransversal[18] = 0.986;
+  anglecorrectionfactortransversal[19] = 0.988;
+  anglecorrectionfactortransversal[20] = 0.99;
+  anglecorrectionfactortransversal[21] = 0.992;
+  anglecorrectionfactortransversal[22] = 0.994;
+  anglecorrectionfactortransversal[23] = 0.996;
+  anglecorrectionfactortransversal[24] = 0.998;
+  for(int i=25; i<61; i++) anglecorrectionfactortransversal[i] = 1;
+  anglecorrectionfactortransversal[61] = 1.003;
+  anglecorrectionfactortransversal[62] = 1.006;
+  anglecorrectionfactortransversal[63] = 1.009;
+  anglecorrectionfactortransversal[64] = 1.012;
+  anglecorrectionfactortransversal[65] = 1.015;
+  anglecorrectionfactortransversal[66] = 1.022;
+  anglecorrectionfactortransversal[67] = 1.029;
+  anglecorrectionfactortransversal[68] = 1.036;
+  anglecorrectionfactortransversal[69] = 1.043;
+  anglecorrectionfactortransversal[70] = 1.05;
+  anglecorrectionfactortransversal[71] = 1.07;
+  anglecorrectionfactortransversal[72] = 1.09;
+  anglecorrectionfactortransversal[73] = 1.11;
+  anglecorrectionfactortransversal[74] = 1.13;
+  anglecorrectionfactortransversal[75] = 1.15;
+  anglecorrectionfactortransversal[76] = 1.175;
+  anglecorrectionfactortransversal[77] = 1.205;
+  anglecorrectionfactortransversal[78] = 1.240;
+  anglecorrectionfactortransversal[79] = 1.280;
+  anglecorrectionfactortransversal[80] = 1.285;
+  anglecorrectionfactortransversal[81] = 1.265;
+  anglecorrectionfactortransversal[82] = 1.20;
+  anglecorrectionfactortransversal[83] = 1.07;
+  anglecorrectionfactortransversal[84] = 0.95;
+  anglecorrectionfactortransversal[85] = 0.8;
+  anglecorrectionfactortransversal[86] = 0.64;
+  anglecorrectionfactortransversal[87] = 0.48;
+  anglecorrectionfactortransversal[88] = 0.32;
+  anglecorrectionfactortransversal[89] = 0.16;
+  anglecorrectionfactortransversal[90] = 0;
+  
+  for(int i=0; i<25; ++i) anglecorrectionfactorlongitudinal[i] = 1;
+  anglecorrectionfactorlongitudinal[26] = 0.999;
+  anglecorrectionfactorlongitudinal[27] = 0.998;
+  anglecorrectionfactorlongitudinal[28] = 0.997;
+  anglecorrectionfactorlongitudinal[29] = 0.996;
+  anglecorrectionfactorlongitudinal[30] = 0.995;
+  anglecorrectionfactorlongitudinal[31] = 0.992;
+  anglecorrectionfactorlongitudinal[32] = 0.989;
+  anglecorrectionfactorlongitudinal[33] = 0.986;
+  anglecorrectionfactorlongitudinal[34] = 0.983;
+  anglecorrectionfactorlongitudinal[35] = 0.98;
+  anglecorrectionfactorlongitudinal[36] = 0.978;
+  anglecorrectionfactorlongitudinal[37] = 0.976;
+  anglecorrectionfactorlongitudinal[38] = 0.974;
+  anglecorrectionfactorlongitudinal[39] = 0.972;
+  anglecorrectionfactorlongitudinal[40] = 0.97;
+  anglecorrectionfactorlongitudinal[41] = 0.968;
+  anglecorrectionfactorlongitudinal[42] = 0.966;
+  anglecorrectionfactorlongitudinal[43] = 0.964;
+  anglecorrectionfactorlongitudinal[44] = 0.962;
+  anglecorrectionfactorlongitudinal[45] = 0.96;
+  anglecorrectionfactorlongitudinal[46] = 0.958;
+  anglecorrectionfactorlongitudinal[47] = 0.956;
+  anglecorrectionfactorlongitudinal[48] = 0.954;
+  anglecorrectionfactorlongitudinal[49] = 0.952;
+  anglecorrectionfactorlongitudinal[50] = 0.95;
+  anglecorrectionfactorlongitudinal[51] = 0.948;
+  anglecorrectionfactorlongitudinal[52] = 0.946;
+  anglecorrectionfactorlongitudinal[53] = 0.944;
+  anglecorrectionfactorlongitudinal[54] = 0.942;
+  anglecorrectionfactorlongitudinal[55] = 0.94;
+  anglecorrectionfactorlongitudinal[56] = 0.938;
+  anglecorrectionfactorlongitudinal[57] = 0.936;
+  anglecorrectionfactorlongitudinal[58] = 0.934;
+  anglecorrectionfactorlongitudinal[59] = 0.932;
+  anglecorrectionfactorlongitudinal[60] = 0.93;
+  anglecorrectionfactorlongitudinal[61] = 0.922;
+  anglecorrectionfactorlongitudinal[62] = 0.914;
+  anglecorrectionfactorlongitudinal[63] = 0.906;
+  anglecorrectionfactorlongitudinal[64] = 0.898;
+  anglecorrectionfactorlongitudinal[65] = 0.89;
+  anglecorrectionfactorlongitudinal[66] = 0.88;
+  anglecorrectionfactorlongitudinal[67] = 0.87;
+  anglecorrectionfactorlongitudinal[68] = 0.86;
+  anglecorrectionfactorlongitudinal[69] = 0.85;
+  anglecorrectionfactorlongitudinal[70] = 0.84;
+  anglecorrectionfactorlongitudinal[71] = 0.824;
+  anglecorrectionfactorlongitudinal[72] = 0.808;
+  anglecorrectionfactorlongitudinal[73] = 0.792;
+  anglecorrectionfactorlongitudinal[74] = 0.776;
+  anglecorrectionfactorlongitudinal[75] = 0.76;
+  anglecorrectionfactorlongitudinal[76] = 0.735;
+  anglecorrectionfactorlongitudinal[77] = 0.71;
+  anglecorrectionfactorlongitudinal[78] = 0.685;
+  anglecorrectionfactorlongitudinal[79] = 0.66;
+  anglecorrectionfactorlongitudinal[80] = 0.635;
+  anglecorrectionfactorlongitudinal[81] = 0.583;
+  anglecorrectionfactorlongitudinal[82] = 0.531;
+  anglecorrectionfactorlongitudinal[83] = 0.479;
+  anglecorrectionfactorlongitudinal[84] = 0.427;
+  anglecorrectionfactorlongitudinal[85] = 0.375;
+  anglecorrectionfactorlongitudinal[86] = 0.3;
+  anglecorrectionfactorlongitudinal[87] = 0.225;
+  anglecorrectionfactorlongitudinal[88] = 0.15;
+  anglecorrectionfactorlongitudinal[89] = 0.075;
+  anglecorrectionfactorlongitudinal[90] = 0;
+  
+  for(int i=0; i<25; ++i) anglecorrectionfactorflat[i] = 1;
+  anglecorrectionfactorflat[26] = 0.999;
+  anglecorrectionfactorflat[27] = 0.998;
+  anglecorrectionfactorflat[28] = 0.997;
+  anglecorrectionfactorflat[29] = 0.996;
+  anglecorrectionfactorflat[30] = 0.995;
+  anglecorrectionfactorflat[31] = 0.992;
+  anglecorrectionfactorflat[32] = 0.989;
+  anglecorrectionfactorflat[33] = 0.986;
+  anglecorrectionfactorflat[34] = 0.983;
+  anglecorrectionfactorflat[35] = 0.98;
+  anglecorrectionfactorflat[36] = 0.982;
+  anglecorrectionfactorflat[37] = 0.974;
+  anglecorrectionfactorflat[38] = 0.966;
+  anglecorrectionfactorflat[39] = 0.958;
+  anglecorrectionfactorflat[40] = 0.95;
+  anglecorrectionfactorflat[41] = 0.946;
+  anglecorrectionfactorflat[42] = 0.942;
+  anglecorrectionfactorflat[43] = 0.938;
+  anglecorrectionfactorflat[44] = 0.934;
+  anglecorrectionfactorflat[45] = 0.93;
+  anglecorrectionfactorflat[46] = 0.924;
+  anglecorrectionfactorflat[47] = 0.918;
+  anglecorrectionfactorflat[48] = 0.912;
+  anglecorrectionfactorflat[49] = 0.906;
+  anglecorrectionfactorflat[50] = 0.90;
+  anglecorrectionfactorflat[51] = 0.89;
+  anglecorrectionfactorflat[52] = 0.88;
+  anglecorrectionfactorflat[53] = 0.87;
+  anglecorrectionfactorflat[54] = 0.86;
+  anglecorrectionfactorflat[55] = 0.85;
+  anglecorrectionfactorflat[56] = 0.837;
+  anglecorrectionfactorflat[57] = 0.824;
+  anglecorrectionfactorflat[58] = 0.811;
+  anglecorrectionfactorflat[59] = 0.798;
+  anglecorrectionfactorflat[60] = 0.785;
+  anglecorrectionfactorflat[61] = 0.77;
+  anglecorrectionfactorflat[62] = 0.755;
+  anglecorrectionfactorflat[63] = 0.74;
+  anglecorrectionfactorflat[64] = 0.725;
+  anglecorrectionfactorflat[65] = 0.71;
+  anglecorrectionfactorflat[66] = 0.693;
+  anglecorrectionfactorflat[67] = 0.676;
+  anglecorrectionfactorflat[68] = 0.659;
+  anglecorrectionfactorflat[69] = 0.642;
+  anglecorrectionfactorflat[70] = 0.625;
+  anglecorrectionfactorflat[71] = 0.605;
+  anglecorrectionfactorflat[72] = 0.585;
+  anglecorrectionfactorflat[73] = 0.565;
+  anglecorrectionfactorflat[74] = 0.545;
+  anglecorrectionfactorflat[75] = 0.525;
+  anglecorrectionfactorflat[76] = 0.501;
+  anglecorrectionfactorflat[77] = 0.477;
+  anglecorrectionfactorflat[78] = 0.453;
+  anglecorrectionfactorflat[79] = 0.429;
+  anglecorrectionfactorflat[80] = 0.405;
+  anglecorrectionfactorflat[81] = 0.371;
+  anglecorrectionfactorflat[82] = 0.337;
+  anglecorrectionfactorflat[83] = 0.303;
+  anglecorrectionfactorflat[84] = 0.269;
+  anglecorrectionfactorflat[85] = 0.235;
+  anglecorrectionfactorflat[86] = 0.188;
+  anglecorrectionfactorflat[87] = 0.141;
+  anglecorrectionfactorflat[88] = 0.094;
+  anglecorrectionfactorflat[89] = 0.047;
+  anglecorrectionfactorflat[90] = 0;
+}
+
+void House::calctemp() {
+  float temperatureday[12], temperaturenight[12];
+
+  temperatureday[0] = 3;
+  temperatureday[1] = 4;
+  temperatureday[2] = 8;
+  temperatureday[3] = 13;
+  temperatureday[4] = 18;
+  temperatureday[5] = 20;
+  temperatureday[6] = 22;
+  temperatureday[7] = 22;
+  temperatureday[8] = 18;
+  temperatureday[9] = 13;
+  temperatureday[10] = 7;
+  temperatureday[11] = 4;
+  temperaturenight[0] = -1;
+  temperaturenight[1] = -1;
+  temperaturenight[2] = 2;
+  temperaturenight[3] = 4;
+  temperaturenight[4] = 8;
+  temperaturenight[5] = 11;
+  temperaturenight[6] = 13;
+  temperaturenight[7] = 13;
+  temperaturenight[8] = 10;
+  temperaturenight[9] = 7;
+  temperaturenight[10] = 3;
+  temperaturenight[11] = 0;
+  for(int i=0; i<12; ++i) {
+    for(int j=0; j<24; ++j) {
+      tempoutside[i][j] = ((temperatureday[i]-temperaturenight[i])/2)*
+                           (cos(((j-12)/24)*2*M_PI)+1)+temperaturenight[i];
+    }
+  }
+  tempinside[0] = 18;
+  tempinside[1] = 18;
+  tempinside[2] = 18;
+  tempinside[3] = 18;
+  tempinside[4] = 18;
+  tempinside[5] = 18;
+  tempinside[6] = 18;
+  tempinside[7] = 20;
+  tempinside[8] = 20;
+  tempinside[9] = 20;
+  tempinside[10] = 20;
+  tempinside[11] = 20;
+  tempinside[12] = 20;
+  tempinside[13] = 20;
+  tempinside[14] = 20;
+  tempinside[15] = 20;
+  tempinside[16] = 20;
+  tempinside[17] = 20;
+  tempinside[18] = 20;
+  tempinside[19] = 20;
+  tempinside[20] = 20;
+  tempinside[21] = 20;
+  tempinside[22] = 20;
+  tempinside[23] = 18;
 }
diff --git a/House.h b/House.h
index 67d90f2..0278726 100644
--- a/House.h
+++ b/House.h
@@ -11,15 +11,23 @@ class House
 {
   public:
     House(float inittemp0, int earthx0, int earthy0, int earthz0, float collectortemp0,
-          float latitude0);
+          float radiation0, float latitude0, float eta00, float k10, float k20, float collectorarea0,
+          float watthoursperhuman, int humans);
     std::list<Triangle> getTriangles(){
       std::list<Triangle> triangles = temp.getTriangles();
       return triangles;
     }
+
+    void calcefficiency();
+    void calctemp();
     
   private:
     Temparray temp;
-    float collectortemp, latitude;
+    float collectortemp, latitude, radiation, eta0, k1, k2, watthoursperhuman, collectorarea;
+    float tempinside[24], tempoutside[12][24], efficiencyflat[100], efficiencyvacuum[100],
+          anglecorrectionfactorflat[91], anglecorrectionfactortransversal[91],
+          anglecorrectionfactorlongitudinal[91];
+    int humans;
 };
 
 #endif /* _HOUSE_H_ */
diff --git a/Sun.h b/Sun.h
index 218f7fa..e2191c9 100644
--- a/Sun.h
+++ b/Sun.h
@@ -10,12 +10,26 @@ class Sun
     
     float declination(int day) {return 23.4*sin(2*M_PI*(284+day)/365);}
     float hourangle(float time) {return 15*(time-12);}
-    float evelationangle(int day, float time) {return asin(sin(declination(day)*M_PI/180)*
-                                                           sin(latitude*M_PI/180)*
-                                                           cos(hourangle(time)*M_PI/180)*
-                                                           cos(declination(day)*M_PI/180)*
-                                                           cos(latitude*M_PI/180))
-                                                           *180/M_PI;}
+    float evelationangle(int day, float time) {
+      return asin(sin(declination(day)*M_PI/180)*
+             sin(latitude*M_PI/180)*
+             cos(hourangle(time)*M_PI/180)*
+             cos(declination(day)*M_PI/180)*
+             cos(latitude*M_PI/180))
+             *180/M_PI;
+    }
+    float azimutangle(int day, float time){
+      float tau, sin_tau, cos_tau, sin_phi, cos_phi, tan_dec;
+
+      tau = hourangle(time)*M_PI/180;
+      sin_tau = sin(tau);
+      cos_tau = cos(tau);
+      sin_phi = sin(latitude*M_PI/180);
+      cos_phi = cos(latitude*M_PI/180);
+      tan_dec = tan(declination(day)*M_PI/180);
+      
+      return atan2(sin_tau, cos_tau*sin_phi-tan_dec*cos_phi)/180*M_PI;      
+    }
   private:
     float latitude;
 };
diff --git a/Temparray.cpp b/Temparray.cpp
index 6385434..83008cc 100644
--- a/Temparray.cpp
+++ b/Temparray.cpp
@@ -35,7 +35,7 @@ Temparray::Temparray(float initialtemp, int x0, int y0, int z0){
       }
     }
   }
-  cubehole(2, 2, 2, 1).setColor(vmml::vec4f(0.0, 0.0, 0.0, 1.0),
+/*  cubehole(2, 2, 2, 1).setColor(vmml::vec4f(0.0, 0.0, 0.0, 1.0),
                                 vmml::vec4f(0.0, 0.0, 0.0, 1.0),
                                 vmml::vec4f(0.0, 0.0, 0.0, 1.0),
                                 vmml::vec4f(0.0, 0.0, 0.0, 1.0));
@@ -46,7 +46,7 @@ Temparray::Temparray(float initialtemp, int x0, int y0, int z0){
   cubehole(2, 2, 2, 3).setColor(vmml::vec4f(0.5, 0.5, 0.5, 1.0),
                                 vmml::vec4f(0.5, 0.5, 0.5, 1.0),
                                 vmml::vec4f(0.5, 0.5, 0.5, 1.0),
-                                vmml::vec4f(0.5, 0.5, 0.5, 1.0));
+                                vmml::vec4f(0.5, 0.5, 0.5, 1.0));*/
 }
 
 void Temparray::calcTemp(){
@@ -256,6 +256,27 @@ void Temparray::calcTemp(){
                 temperaturenew(i, j, k, l, 0) + ((temperatureold(i, j, k, l, 3) -
                   temperatureold(i, j, k, l, 0))*(1-exp(-(1/tau21))));
             }
+            if(j == sy-1) {
+              if(m % 2 == 0)area = (((depth - innerdepth)/2)*((width + innerwidth)/2))
+                                    / 0.9 * meterperunit;
+              if(m % 2 == 1)area = (((width - innerwidth)/2)*((depth + innerdepth)/2))
+                                    / 0.9 * meterperunit;
+
+              thermalresistance = 3 + (height/2 / 0.9 * meterperunit /(conductivity * area));
+    
+              if(m % 2 == 0)capacity = specificcapacity *((((depth - innerdepth)/2)*
+                                    ((width + innerwidth)/2)*height)
+                                    / 0.9 * meterperunit);
+              if(m % 2 == 1)capacity = specificcapacity *((((width - innerwidth)/2)*
+                                    ((depth + innerdepth)/2)*height)
+                                    / 0.9 * meterperunit);
+                                    
+              tau12 = capacity * thermalresistance;
+              
+              temperaturenew(i, j, k, 0, m) =
+                temperaturenew(i, j, k, 0, m) - ((temperatureold(i, j, k, 0, m) - earthtemp)
+                  *(1-exp(-(1/tau12))));
+            }
           }
         }
         if(i > 0) {
@@ -270,7 +291,7 @@ void Temparray::calcTemp(){
           thermalresistance = ((width - innerwidth)/2) / 0.9 * meterperunit /
                                (conductivity * area);
 
-          capacity = specificcapacity *((((depth - innerdepth)/2)*((width + innerwidth)/2)*
+          capacity = specificcapacity *((((depth + innerdepth)/2)*((width - innerwidth)/2)*
                                               height)/0.9*meterperunit);
           tau12 = capacity * thermalresistance;
           
@@ -281,6 +302,29 @@ void Temparray::calcTemp(){
             temperaturenew(i-1, j, k, 0, 1) + ((temperatureold(i, j, k, 0, 3) -
               temperatureold(i-1, j, k, 0, 1))*(1-exp(-(1/tau12))));       
         }
+        if (i == 0 || i == sx-1) {
+          int h=3;
+          if(i == sx-1) h=1;
+          
+          width = cubehole(i, j, k, 0).getWidth();
+          height = cubehole(i, j, k, 0).getHeight();
+          depth = cubehole(i, j, k, 0).getDepth();
+          innerwidth = cubehole(i, j, k, 0).getInnerWidth();
+          innerdepth = cubehole(i, j, k, 0).getInnerDepth();
+          
+          area = (((depth + innerdepth)/2)*height) / 0.9 * meterperunit;
+
+          thermalresistance = 3 + ((width - innerwidth)/4) / 0.9 * meterperunit /
+                               (conductivity * area);
+
+          capacity = specificcapacity *((((depth + innerdepth)/2)*((width - innerwidth)/2)*
+                                              height)/0.9*meterperunit);
+          tau12 = capacity * thermalresistance;
+          
+          temperaturenew(i, j, k, 0, h) =
+            temperaturenew(i, j, k, 0, h) - ((temperatureold(i, j, k, 0, h) - earthtemp)
+              *(1-exp(-(1/tau12))));
+        }
         if(k > 0) {
           width = cubehole(i, j, k, 0).getWidth();
           height = cubehole(i, j, k, 0).getHeight();
@@ -304,10 +348,34 @@ void Temparray::calcTemp(){
             temperaturenew(i, j, k-1, 0, 2) + ((temperatureold(i, j, k, 0, 0) -
               temperatureold(i, j, k-1, 0, 2))*(1-exp(-(1/tau12))));       
         }
+        if(k == 0 || k == sz-1) {
+          int h=0;
+          if(i == sz-1) h=2;
+          
+          width = cubehole(i, j, k, 0).getWidth();
+          height = cubehole(i, j, k, 0).getHeight();
+          depth = cubehole(i, j, k, 0).getDepth();
+          innerwidth = cubehole(i, j, k, 0).getInnerWidth();
+          innerdepth = cubehole(i, j, k, 0).getInnerDepth();
+          
+          area = (((width + innerwidth)/2)*height) / 0.9 * meterperunit;
+
+          thermalresistance = 3 + ((depth - innerdepth)/4) / 0.9 * meterperunit /
+                               (conductivity * area);
+
+          capacity = specificcapacity *((((width + innerwidth)/2)*((depth - innerdepth)/2)*
+                                              height)/0.9*meterperunit);
+          tau12 = capacity * thermalresistance;
+          
+          temperaturenew(i, j, k, 0, h) =
+            temperaturenew(i, j, k, 0, h) - ((temperatureold(i, j, k, 0, h) - earthtemp)
+              *(1-exp(-(1/tau12))));
+        }
       }
     }
   }
   mergetemperature();
+  coloring();
 //  std::cerr << temperaturenew(2, 2, 2, 1, 2) << std::endl;
 //  std::cerr << temperaturenew(2, 2, 2, 2, 2) << std::endl;
 //  std::cerr << temperaturenew(2, 2, 2, 3, 2) << "\n" << std::endl;
@@ -329,6 +397,28 @@ std::list<Triangle> Temparray::getTriangles(){
       }
     }
   }
-
   return triangles;
 }
+
+void Temparray::coloring() {
+  float r[4], g[4], b[4];
+  for(int i=0; i<sx; ++i){
+    for(int j=0; j<sy; ++j){
+      for(int k=0; k<sz; ++k){
+        for(int l=0; l<6; ++l){
+          for(int m=0; m<4; ++m){
+            if(temperatureold(i, j, k, l, m) >= 25 ) b[m]=0;
+            else b[m]=(25-temperatureold(i, j, k, l, m))/25;
+            g[m]=fabs((25-temperatureold(i, j, k, l, m))/25);
+            if(temperatureold(i, j, k, l, m) <= 25 ) r[m]=0;
+            else r[m]=(100-temperatureold(i, j, k, l, m))/25;
+          }
+          cubehole(i, j, k, l).setColor(vmml::vec4f(r[0], g[0], b[0], 1.0),
+                                        vmml::vec4f(r[1], g[1], b[1], 1.0),
+                                        vmml::vec4f(r[2], g[2], b[2], 1.0),
+                                        vmml::vec4f(r[3], g[3], b[3], 1.0));
+        }
+      }
+    }
+  }
+}
diff --git a/Temparray.h b/Temparray.h
index 24cb95d..ba4e8e8 100644
--- a/Temparray.h
+++ b/Temparray.h
@@ -20,20 +20,14 @@ class Temparray
     float& temperaturenew(int x, int y, int z, int line, int pos) {
       return temparraynew[x*sy*sz*6*4 + y*sz*6*4 + z*6*4 + line*4 + pos];
     }
-/*    void deactivateCube(int x, int y, int z){
-      for(int i=0; i<6; i++){
-        for(int j=0; j<4; j++){
-          setTemp(-100, x, y, z, i, j);
-        }
-      }
-    }*/
     std::list<Triangle> getTriangles();
+    void coloring();
 
   private:
     float *temparrayold, *temparraynew;
     Cubehole *cubearray;
     int sx, sy, sz;
-    float averagetemp;
+    float averagetemp, earthtemp;
 
     Cubehole& cubehole(int x, int y, int z, int line) {
       return cubearray[x*sy*sz*6 + y*sz*6 + z*6 + line];
@@ -58,28 +52,10 @@ class Temparray
           }
         }
       }
+      earthtemp=8;
       temperatureold(2, 2, 2, 2, 2) = 90;
       temperaturenew(2, 2, 2, 2, 2) = 90;
     }
-/*    void calcAverage(){
-      float tempcache = 0;
-      int times = 0;
-      for(int i=0; i<sx; i++){
-        for(int j=0; j<sy; j++){
-          for(int k=0; k<sz; k++){
-            for(int l=0; l<6; l++){
-              for(int m=0; m<4; m++){
-                if(temperature(i, j, k, l, m) != -100){
-                  tempcache += temperature(i, j, k, l, m);
-                  ++times;
-                }
-              }
-            }
-          }
-        }
-      }
-      averagetemp = tempcache/times;
-    }*/
 };
 
 #endif /* _TEMPARRAY_H_ */
diff --git a/main.cpp b/main.cpp
index b27c596..45a9c88 100644
--- a/main.cpp
+++ b/main.cpp
@@ -12,7 +12,7 @@
 #endif
 
 //#define MIN_FRAME_DELTA 16
-#define MIN_FRAME_DELTA 250
+#define MIN_FRAME_DELTA 41
 #define DEFAULT_WIDTH 640
 #define DEFAULT_HEIGHT 640
 
@@ -203,7 +203,7 @@ int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine
 
   bool running = true;
   MSG msg;
-
+  
   while(running)
   {
     unsigned long delta = 0;