001 package dk.deepthought.sidious.greenhouse;
002
003 /**
004 * @author Michael Rasmussen
005 */
006 public class LeafPhotosynthesisModel
007 {
008 private static final double E_J = 37000;
009 private static final double E_C = 59356;
010 private static final double E_O = 35948;
011 private static final double E_Rd = 66405;
012 private static final double E_VC = 58520;
013 @SuppressWarnings("unused")
014 private static final double Q_10_Rd = 2.21;
015 private static final double r_b_H2O = 100;
016 private static final double S = 710;
017 private static final double H = 220000;
018 private static final double R_d_25 = 1.1;
019 private static final double Theta = 0.7;
020 private static final double R = 8.314;
021 private static final double a_0 = 0.0176715;
022 private static final double V_c_max_25 = 97.875;
023 private static final double J_max_25 = 210.15;
024 private static final double K_O_25 = 155;
025 private static final double K_C_25 = 310;
026 private static final double p_O2i = 210;
027 private static final double r_s_H2O = 50;
028 private static final double T_0 = 273.15;
029 private static final double T_25 = 298.15;
030 private static final double V_O_C = 0.21;
031 private static final double rho_CO2_T_0 = 1.98;
032 private static final double M_CO2 = 44e-3;
033
034 /**
035 * This method calculates the photosynthesis.
036 *
037 * @param temp
038 * inside temperature (UC: 1) [° Celsius]
039 * @param co2
040 * concentration of CO2 (UC: 14) [ppm]
041 * @param sun
042 * irradiance (UC: 56) [Wm^-2]
043 * @param glass_factor
044 * the factor of light that passes through the glass panes (UC:
045 * 158) [%]
046 * @param shade_factor
047 * the factor of light that passes through the screens (UC: 567)
048 * [%]
049 * @param shade
050 * the position of the screens (UC: 590) [%]
051 * @return the calculated photosynthesis
052 */
053 public double calculate(double temp, double co2, double sun, double glass_factor, double shade_factor, double shade)
054 {
055 if (Double.isNaN(temp) || Double.isNaN(co2) || Double.isNaN(sun))
056 return Double.NaN;
057
058 if (Double.isNaN(glass_factor))
059 glass_factor = 0.6; // default: 60%
060 else
061 glass_factor /= 100;
062
063 if (Double.isNaN(shade_factor))
064 shade_factor = 0.4; // default: 40%
065 else
066 shade_factor /= 100;
067
068 if (Double.isNaN(shade))
069 shade = 0.5; // default: 50%
070 else
071 shade /= 100;
072
073 double light = 2.3 * sun * glass_factor * ((1 - shade) + (shade * shade_factor));
074
075 double T_l = temp + T_0;
076 double CO2 = co2;
077 double I_a = light / 4.6; // divide by 4.6, to convert from [µmol m^-2 s^-1] to [J m^-2 s^-1]
078
079 double P_g;
080 double P_g_max;
081 double a_l;
082 double Gamma;
083 double P_n_max;
084 double P_n_CO2;
085 double J_max;
086 double rho_CO2_T_l;
087 double r_CO2;
088 double V_c_max;
089 double r_c_CO2;
090 double R_d;
091
092 double P_n;
093
094 double temp_diff = calculate_temp_diff(T_l);
095
096 Gamma = calculate_Gamma(temp_diff);
097 rho_CO2_T_l = calculate_rho_CO2_T_l(T_l);
098 J_max = calculate_J_max(T_l, temp_diff);
099 R_d = calculate_R_d(T_l, temp_diff);
100 V_c_max = calculate_V_c_max(temp_diff);
101
102 r_c_CO2 = calculate_r_c_CO2(temp_diff, rho_CO2_T_l, V_c_max);
103
104 r_CO2 = calculate_r_CO2(r_c_CO2);
105 P_n_CO2 = calculate_P_n_CO2(CO2, Gamma, rho_CO2_T_l, r_CO2);
106 a_l = calculate_a_l(CO2, Gamma);
107
108 P_n_max = calculate_P_n_max(P_n_CO2, J_max);
109 P_g_max = calculate_P_g_max(P_n_max, R_d);
110
111 P_g = calculate_P_g(P_g_max, a_l, I_a);
112 P_n = calculate_P_n(P_g, R_d);
113
114 P_n *= 22.7222; // multiply by 22.7222 to convert from [mg m^-2 s^-1] to [µmol m^-2 s^-1]
115
116 // multiply by 60, to convert from pr second to pr minute
117 return P_n * 60;// * 0.2133700901; // multiply by 0.21337, for Justicia plant type
118 }
119
120 private double calculate_temp_diff(double T_l) // -- [mol J^-1]
121 {
122 return (T_l - T_25) / (T_l * R * T_25);
123 }
124
125 private double calculate_P_n(double P_g, double R_d) // -- [mg m^-2 s^-1]
126 {
127 return P_g - R_d;
128 }
129
130 private double calculate_P_g_max(double P_n_max, double R_d) // -- [mg m^-2 s^-1]
131 {
132 return P_n_max + R_d;
133 }
134
135 private double calculate_P_g(double P_g_max, double a_l, double I_a) // [1] -- [mg m^-2 s^-1]
136 {
137 double part1 = 1 - Math.exp((-1 * a_l * I_a) / P_g_max);
138 return P_g_max * part1;
139 }
140
141 private double calculate_a_l(double CO2, double Gamma) // [2] -- [mg J^-1]
142 {
143 double max = Math.max(CO2, Gamma);
144 return a_0 * ((max - Gamma) / (max + 2 * Gamma));
145 }
146
147 private double calculate_Gamma(double temp_diff) // [3] -- [µbar] = [µmol mol^-1] = [ppm]
148 {
149 double part1 = (p_O2i * V_O_C) / 2;
150 double part2_t = K_C_25 * Math.exp(E_C * temp_diff);
151 double part2_n = K_O_25 * Math.exp(E_O * temp_diff);
152
153 return part1 * (part2_t / part2_n);
154 }
155
156 private double calculate_P_n_max(double P_n_CO2, double J_max) // [4] -- [mg m^-2 s^-1]
157 {
158 double part1 = Math.sqrt(Math.pow(P_n_CO2 + J_max, 2) - (4 * Theta * P_n_CO2 * J_max));
159 return (P_n_CO2 + J_max - part1) / (2 * Theta);
160 }
161
162 private double calculate_P_n_CO2(double CO2, double Gamma, double rho_CO2_T_l, double r_CO2) // [5] -- [mg m^-2 s^-1]
163 {
164 return ((Math.max(CO2, Gamma) - Gamma) * rho_CO2_T_l) / r_CO2;
165 }
166
167 private double calculate_rho_CO2_T_l(double T_l) // [6] -- [kg m^-3]
168 {
169 return rho_CO2_T_0 * (T_0 / T_l);
170 }
171
172 private double calculate_r_c_CO2(double temp_diff, double rho_CO2_T_l, double V_c_max) // [7] -- [bar s m^-1]
173 {
174 double part1 = K_C_25 * Math.exp(E_C * temp_diff);
175 double part2_n = K_O_25 * Math.exp(E_O * temp_diff);
176 double part2 = 1 + (p_O2i / part2_n);
177
178 return ((part1 * part2 * rho_CO2_T_l) / V_c_max) / M_CO2;
179 }
180
181 private double calculate_r_CO2(double r_c_CO2) // [8] -- [s m^-1] -- OBS: r_c_CO2 ^^
182 {
183 return 1.6 * r_s_H2O + 1.37 * r_b_H2O + r_c_CO2;
184 }
185
186 private double calculate_J_max(double T_l, double temp_diff) // [9] -- [µmol m^-2 s^-1]
187 {
188 double part1 = Math.exp(E_J * temp_diff);
189 double part2 = 1 + Math.exp((S - H / T_25) / R);
190 double part3 = 1 + Math.exp((S - H / T_l) / R);
191
192 return (M_CO2 / 4) * J_max_25 * part1 * (part2 / part3);
193 }
194
195 private double calculate_R_d(double T_l, double temp_diff) // [10] -- [mg m^-2 s^-1]
196 {
197 //return M_CO2 * R_d_25 * Math.pow(Q_10_Rd, (T_l - T_25) / 10);
198 return M_CO2 * R_d_25 * Math.exp(E_Rd * temp_diff);
199 }
200
201 private double calculate_V_c_max(double temp_diff) // [11] -- [µmol m^-2 s^-1]
202 {
203 return V_c_max_25 * Math.exp(E_VC * temp_diff);
204 }
205 }
|