/* Datei streckung.rgg, gehoert zu pappel.gsz */

/* Aenderungen gegenueber frueherer Version (p9):

Wichtig: Die Basiszeiteinheit dieser Simulation liegt bei einer Stunde,
d.h. pro Ausfuehrung des run()-Regelblocks vergeht eine Stunde Simulationszeit.
Definiert wird dies durch die Konstante NUMBER_OF_TIME_UNITS_DURING_A_DAY = 24.
Alle weiteren Angaben sind relativ zu dieser Angabe.

Die Photosyntheseberechnung wurde auf SI-Einheiten umgestellt.
Kommentare wurden ergaenzt.
Die Statusvariable fuer das Oeffnen des Datenfiles wurde geaendert.

*/
import java.io.*;
// Import der zweiten rgg-Datei:
import static photosynthese.*;

/**************************************************************************/
/************************* Modul-Definitionen des Baumes ******************/
/**************************************************************************/ 


// Der ganze Baum:
module Tree(int age, float pspool) extends Sphere(0.01);

// ein einfaches Internodium mit einer festen, sehr kleinen Startgroesse:
module Internode(int alter, int lnfb) extends F(0.01, 0.0005)
   {
   // Variable zur Speicherung der freien Assimilate:
   float pspool = 0;
   float mass = 0;
   { setShader(GREEN); }
	
   public double getVolume()
      { 
      return Math.PI*(this.diameter)*(this.diameter) / 4 * this.length;
      }
   public double getMass()
      { 
      return this.getVolume()*INTERNODIUM_DENSITY;
      }
}

// ein Meristem als Keimzelle des Wachstums:
module Meristem(int alter, int lnfb) extends Sphere(0.0005)
   {{ setShader(RED); }};
   // lnfb = leaf number from base

// der Blattstiel:
module Petiole(int alter, int lnfb) extends F(0.001, 0.0001)
   {
   // Variable zur Speicherung der freien Assimilate:
   float pspool = 0.0;
   float mass = 0.0;
   { setShader(GREEN); }
   
   public double getVolume()
      { 
      return Math.PI*(this.diameter)*(this.diameter) / 4 * this.length;
      }
   public double getMass()
      { 
      return this.getVolume()*PETIOLE_DENSITY;
      }
   }


// das Blatt: die Groesse wird hier ueber Scale gesetzt.
// Basis der Scale-Berechnung ist 1 Meter Laenge und 0.77 Meter Breite 
// (definiert durch LEAF_LENGTH_WIDTH_RATIO).
// Das Verhaeltnis ist zur Zeit unveraenderlich!
module Blade(int alter, int lnfb) 
   extends Parallelogram(1.0, 1.0/LEAF_LENGTH_WIDTH_RATIO)
   { 
   //Variable zur Speicherung der freien Assimilate
   float pspool = 0;
   float leafarea = 0;
   float mass = 0;
	
   public double getMass()
      { 
      return this.leafarea*LEAF_WEIGHT_PER_SQUAREMETER;
      }
   };

// Module zum Veraendern der Winkel (zwischen Internodium und Blattstiel 
// und zwischen Blattsiel und Blattspreite):
// Die Eigenschaften des Moduls werden durch eine Aktualisierungsregel 
// veraendert, so dass sich die Position der Blaetter waehrend der Simulation 
// veraendert.
module Winkel(super.angle, int max_angle) extends RL(angle); 


/***************  Die Sonne: **************************************/

// Die Sonne, sie umkreist den Baum einmal pro Tag ...
module Sonne(float ppfd) extends Sphere(3.3)
   {{ setShader(RGBAShader.YELLOW);}};

// Hilfssymbol zur Bewegung der Sonne:
module Sonnenwinkel(super.angle) extends RU(angle); 


/**************************************************************/
/************** Konstanten und globale Variablen: *************/
/**************************************************************/

const float NUMBER_OF_TIME_UNITS_DURING_A_DAY = 24; 	
// dient der Umrechnung von Tageswerten (Photosynthese) auf die 
// Bezugsgroesse Basiszeiteinheiten.


const Shader leafmat = shader("Leaf1");  // Textur fuer die Blaetter
const PrintWriter tmpfile;	// Ausgabespeicher: -> Ausgabe in Datei
private Boolean tmpFileOpen = false; 
   // Kontrollvariable, ob Datenfile offen ist.

const Vector3d zz = new Vector3d(0, 0, 1); // Richtung nach oben

const int PLASTOCHRON = 48;
   // Anzahl der Zeiteinheiten zwischen dem Erscheinen zweier Blaetter

// Maximale Wachstumszeiten fuer die drei Metamer-Bestandteile: 
// koennen fuer Stopbedingungen fuer das Wachstum genutzt werden.
const int MAX_GROWTHTIME_INTERNODE = 600;    // in Basiszeiteinheiten
const int MAX_GROWTHTIME_PETIOLE = 30; 	     // in Basiszeiteinheiten
const int MAX_GROWTHTIME_LEAF = 96;          // in Basiszeiteinheiten


// Wachstumsraten fuer die drei Metamer-Bestandteile
// Wachstum pro Zeiteinheit 
const float GROWTH_RATE_INTERNODE_LENGTH = 0.001;       // in Metern
const float GROWTH_RATE_INTERNODE_DIAMETER = 0.000004;  // in Metern
const float GROWTH_RATE_PETIOLE_LENGTH = 0.00117;       // in Metern
const float GROWTH_RATE_PETIOLE_DIAMETER = 0.000031;    // in Metern
const float GROWTH_RATE_LEAF_SCALE = 0.0016; 
   // betrifft beide Achsen: Skalierungsfaktor der Basisgroesse des 
   // Blattes von 1 Meter * 0.77 Meter

// Definition von 3 Wachstumszonen am Baum:
// Die oberste Wachstumszone wird definiert durch den maximalen 
//    Leaf Position Index (LPI): -1 bis TOPSECTION_MAXLPI.
// Die unterste Wachstumszone durch die angegebene maximale 
//    Leaf Number From Base (LNFB): 1 bis BASESECTION_MAXLEAFNUMBER.
// Die mittlere Wachstumszone umfasst alle anderen Blaetter.
// Die Wachstumszone wird errechnet durch die Servicefunktion 
//    getGrowthSection()
const int TOPSECTION = 0;
const int TOPSECTION_MAXLPI = 8;
const int MIDDLESECTION= 1;
const int BASESECTION = 2;
const int BASESECTION_MAXLEAFNUMBER = 5;

// maximale Dimensionen fuer die 3 Metamerbestandteile:
// die Werte werden als Arry getrennt fuer die drei Wachstumszonen 
//    angegeben
const float[] MAX_GWROWTH_INTERNODE = {0.01, 0.03, 0.02};   // in Metern
const float[] MAX_GWROWTH_PETIOLE = {0.01, 0.07, 0.02};     // in Metern
const float[] MAX_GWROWTH_LEAF_SCALE = {0.07, 0.09, 0.03};  // in Metern

const float LEAF_LENGTH_WIDTH_RATIO = 1.3;
   // konstantes Laengen / Breitenverhaeltnis der Blaetter
const float LEAFFORMFACTOR = 0.6;
   // Flaechenanteil der realen Blattflaeche an Laenge * Breite

const float INTERNODIUM_DENSITY = 106.0;        // in kg/m^3
const float PETIOLE_DENSITY = 106.0;            // in kg/m^3
const float LEAF_WEIGHT_PER_SQUAREMETER =100.0; // in kg/m^2


/***************** Winkel: ****************/

const int LEAF_MAX_ANGLE = 55;
const int PETIOLE_MAX_ANGLE = 90;
// Maximaler Divergenzwinkel zwischen zwei Internodien:
const int MAX_ORIENTATION_ANGLE = 5; 
// Azimuthwinkel fuer den Ansatzpunkt der Blaetter am Stamm:
const int MAX_AZIMUTH = 160; 
const int MIN_AZIMUTH = 120; 

// Seitliche Drehung der Blattspreite um die Blattstielachse (nicht gemessen)
const int MAX_ANGLE_LAMINA = 5;    //angle of lamina



/******************** globale Variablen: ***************************/
public float totalPSPool = 0.0; 
public int maxleaf = 0; // aktuelle Zahl aller Blaetter des Baumes
public int time = 0;    // Zaehlvariable fuer aktuelle Zeit



/***************************************************************/
/**************** Regelbloecke: ********************************/
/***************************************************************/

protected void init()
   {
   time = 0;
   maxleaf = 0;
      [
      Axiom ==> Tree(0, 0.03) Meristem(0, 0) 
                   { totalPSPool = 0.03; },
		^ Sonnenwinkel(0) RU(90.0) M(50.0) Sonne(0);
      ]
   }

// Zentrale Methode zur Steuerung des Simulationslaufes:
// Ruft die Regelbloecke fuer die Photosntheseberechnung, den Transport 
// der freien Assimilate und das Wachstum der Metamere und der Bewegung 
// der Sonne auf. Basis der zeitlichen Abfolge ist die Variable time.

public void run()
   {
   if (!tmpFileOpen) newDataFile();
	
   println("Stunde:" + time + "  Gesamt-Produktion: " + totalPSPool); 
   tmpfile.println("Stunde" + time + "Gesamt-Produktion: " + totalPSPool); 
	
   for (apply(1)) moveSun(360 / NUMBER_OF_TIME_UNITS_DURING_A_DAY); 
   for (apply(1)) photosynthese(NUMBER_OF_TIME_UNITS_DURING_A_DAY);
   for (apply(4)) transport();
   for (apply(1)) grow();
   tmpfile.flush(); // Schreiben des Speichers tmpfiles in die Datei.
   time ++;
   }

// Regel zur Kalkulation der Photosynthese in den Blaettern:
// Die Photosyntheseberechnung erfolgt immer fuer einen Tagesgang.
// Daher ist es noetig, den Anteil fuer die Zeitschritte zu berechnen.
// Hierzu wird durch number_of_time_units geteilt.

private void photosynthese(float number_of_time_units)
   {
   tmpfile.println("   Photosynthese"); 
      [	
      s:Blade ::> 
         {
         int lpi = getLPI(maxleaf, s[lnfb]);
         s[leafarea] = s[scale] * s[scale] / LEAF_LENGTH_WIDTH_RATIO * LEAFFORMFACTOR;
         println("blattflaeche" + s[leafarea]);
         float production = 
            calculatePS(s[leafarea], convertMC(lpi)) / number_of_time_units;
         s[pspool] += production;
         totalPSPool += production;
		
         // tmpfile.println("      Blatt: " + s[lnfb] + " LPI " + lpi + 
         //  "  Pool: " + s[pspool] + "  Blattflaeche: " + s[leafarea]);
         }
      ]
   }

// Regelblock fuer das Wachstum der Pflanze:
// es existieren Unter-Regelbloecke fuer das Meristem, Internodiium, 
// Petiole, Blade und Winkel.
private void grow()
   {
   // Regelblock fuer das Meristem:
      [   
      // erste Regel: Das Meristem altert
      m1:Meristem ::> m1[alter] :+= 1;	
      // fuer ein neues Meristem: Bildung des Internodiums:
      m2:Meristem, ((m2[alter] == PLASTOCHRON) || (m2[lnfb] == 0)) ==> 
          Internode(0, m2[lnfb]+1) [ Winkel(0, PETIOLE_MAX_ANGLE) 
             Petiole(0, m2[lnfb]+1) RH(random(MAX_ANGLE_LAMINA, -MAX_ANGLE_LAMINA)) 
             Winkel(1,LEAF_MAX_ANGLE) 
                [ Blade(0, m2[lnfb]+1).(setShader(leafmat), setScale(0.001))
                  { maxleaf++; } ] ]
          RH(random(MAX_AZIMUTH, MIN_AZIMUTH)) 
          RU(random(MAX_ORIENTATION_ANGLE, -MAX_ORIENTATION_ANGLE)) 
          Meristem(0, m2[lnfb]+1);
      ]

   // Regelblock fuer das Veraendern der Winkel von Blattstiel und Blatt:
      [
      x:Winkel, (x[angle] <= x[max_angle]) ::> x[angle] :+= 1;
      ]

   // Regelblock fuer die Streckung des Internodiums
      [
      i:Internode ::> 
         {
         if (i[alter] < MAX_GROWTHTIME_INTERNODE && 
            i[length] <= MAX_GWROWTH_INTERNODE[getGrowthSection(maxleaf, i[lnfb])]) 
            { 
            i[length] += GROWTH_RATE_INTERNODE_LENGTH; 
            }  // Laengenwachstum in den ersten 30 Tagen
         i[diameter] += GROWTH_RATE_INTERNODE_DIAMETER; // Dickenwachstum immer
         i[alter]++;               // Hochsetzen des Alters
         i[mass] = i.getMass();	
         }
      ]

   // Regelblock fuer die Streckung des Blattstiels in den ersten 60 Tagen:
      [
      p:Petiole, (p[alter] < MAX_GROWTHTIME_PETIOLE && 
         p[length] < MAX_GWROWTH_PETIOLE[getGrowthSection(maxleaf, p[lnfb])]) ::> 
            {
            p[length] += GROWTH_RATE_PETIOLE_LENGTH;     // Laengenwachstum 
            p[diameter] += GROWTH_RATE_PETIOLE_DIAMETER; // Dickenwachstum 
            p[alter]++;    // Hochsetzen des Alters
            p[mass] = p.getMass();
            }
      ]

   // Regelblock fuer das Blattwachstum:
      [
      s:Blade, (s[alter] < MAX_GROWTHTIME_LEAF && 
         s[scale] < MAX_GWROWTH_LEAF_SCALE[getGrowthSection(maxleaf, s[lnfb])]) ::> 
            {
            s.setScale(s[scale] + GROWTH_RATE_LEAF_SCALE);
            s[alter]++;
            s[mass] = s.getMass();
            }	
      ]
   }

// Regelblock fuer den Transport von Assimilaten innerhalb der Pflanze:
// 	es existieren mehrere Transportwege
//	Blade (Quelle) -> zum tragenden Internodium
//	Internodium -> zum benachbarten Internodium
//	Tree (Keimling) -> zum ersten Infernodium


private void transport()
   {
   tmpfile.println("  Transport " ); 
   // Transport Blade -> zum tragenden Internodiium
      [
      b:Blade -ancestor-> i:Internode ::> 
         {
         if (getLPI(maxleaf, b[lnfb]) >= 4)
            {
            float r = 0.01 * (b[pspool] - i[pspool]);
//          println("basipetaler PS-transport: " + r + " (rule 2)");
            b[pspool] :-= r;
            i[pspool] :+= r; 
            tmpfile.println("       Blatt: " + b[lnfb] + " zu Internode " + 
               i[lnfb] + "  Menge: " + r);
            }
         }
      ]
	
   // Transport Internodium -> zum benachbarten Internodium
      [
      i_oben:Internode -ancestor-> i_unten:Internode ::> 
         {
         float r = 0.01 * (i_unten[pspool] - i_oben[pspool]);
         tmpfile.println("        Transport von Internodium: " + i_unten[lnfb] + 
            "zu : " + i_oben[lnfb] + " Menge: " + r);
         i_unten[pspool] :-= r;
         i_oben[pspool] :+= r;		
         }
      ]
	
   // Transport Tree -> zum ersten Internodium
      [
      t:Tree -successor-> i:Internode, (t[pspool] > 0.0001) ::>
         {
         float r = 0.001 * t[pspool];
         tmpfile.println("        Transport von Treepool zu Internodium: " + 
            i[lnfb] + " Menge: " + r);
         t[pspool] :-= r;
         i[pspool] :+= r;		
         }
      ]
	
/* Hier moegliche Ergaenzung Transport zum Meristem: Sinnvoll bei Kopplung der 
   Metamerbildung an vorhandene Assimilate  */
   }

private void moveSun(float sonnenbewegung)
   [
   // Bewegung der Sonne: Hat nur Bedeutung für die Animation. 
   // Es existiert kein Bezug zur numerischen Simulation.
   x:Sonnenwinkel ::> x[angle] :+= sonnenbewegung;
   ]

/********************** Methoden mit Zugriff auf die Struktur ***********/
// Beschattungsfunktion: prueft, ob es fuer einen Knoten "s"
// innerhalb eines Sichtkegels mit der Oeffnungsweite 40 Grad (kann man veraendern)
// ein Objekt f der Modul-Klasse Petiole (Blattstiel) gibt, wobei "s" (als 
// "Beobachter" sozusagen) nicht beruecksichtigt wird. Rueckgabewert ist false
// oder true
private static boolean isShaded(Node s)
   {
   return !empty( (* f:Petiole, (f != s && f in cone(s, zz, 40)) *) );
   } 



/********************** Methoden zur Berechnung ************************/

private int getLPI(int maxLeaf, int lnfb)
   {
   return maxLeaf - lnfb - 1; 
   }

private int getGrowthSection(int maxleaf, int lnfb)
   {
   int section = 0;
   if (lnfb <= BASESECTION_MAXLEAFNUMBER)
      section = BASESECTION;
   else
      if(getLPI(maxleaf, lnfb) <= TOPSECTION_MAXLPI)
         section = TOPSECTION;
      else
         section = MIDDLESECTION;
   return section;	
   }

/********************* Methoden fuer Serviveausgaben ********************/

private void newDataFile()
   {
   tmpfile = new PrintWriter(new FileWriter("c:\\tmp\\dataoutput3.txt"));
   tmpFileOpen = true;
   }

//--------  Spezielle Methode (public) zum Schliessen der Datei
public void closeDataFile()
   {
   tmpfile.close();	// Schliessen der Datei
   tmpFileOpen = false;
   }