| 1 | package dk.deepthought.sidious.greenhouse; |
| 2 | |
| 3 | /** |
| 4 | * @author Michael Rasmussen |
| 5 | */ |
| 6 | public class LeafPhotosynthesisModel |
| 7 | { |
| 8 | private static final double E_J = 37000; |
| 9 | private static final double E_C = 59356; |
| 10 | private static final double E_O = 35948; |
| 11 | private static final double E_Rd = 66405; |
| 12 | private static final double E_VC = 58520; |
| 13 | @SuppressWarnings("unused") |
| 14 | private static final double Q_10_Rd = 2.21; |
| 15 | private static final double r_b_H2O = 100; |
| 16 | private static final double S = 710; |
| 17 | private static final double H = 220000; |
| 18 | private static final double R_d_25 = 1.1; |
| 19 | private static final double Theta = 0.7; |
| 20 | private static final double R = 8.314; |
| 21 | private static final double a_0 = 0.0176715; |
| 22 | private static final double V_c_max_25 = 97.875; |
| 23 | private static final double J_max_25 = 210.15; |
| 24 | private static final double K_O_25 = 155; |
| 25 | private static final double K_C_25 = 310; |
| 26 | private static final double p_O2i = 210; |
| 27 | private static final double r_s_H2O = 50; |
| 28 | private static final double T_0 = 273.15; |
| 29 | private static final double T_25 = 298.15; |
| 30 | private static final double V_O_C = 0.21; |
| 31 | private static final double rho_CO2_T_0 = 1.98; |
| 32 | private static final double M_CO2 = 44e-3; |
| 33 | |
| 34 | /** |
| 35 | * This method calculates the photosynthesis. |
| 36 | * |
| 37 | * @param temp |
| 38 | * inside temperature (UC: 1) [° Celsius] |
| 39 | * @param co2 |
| 40 | * concentration of CO2 (UC: 14) [ppm] |
| 41 | * @param sun |
| 42 | * irradiance (UC: 56) [Wm^-2] |
| 43 | * @param glass_factor |
| 44 | * the factor of light that passes through the glass panes (UC: |
| 45 | * 158) [%] |
| 46 | * @param shade_factor |
| 47 | * the factor of light that passes through the screens (UC: 567) |
| 48 | * [%] |
| 49 | * @param shade |
| 50 | * the position of the screens (UC: 590) [%] |
| 51 | * @return the calculated photosynthesis |
| 52 | */ |
| 53 | public double calculate(double temp, double co2, double sun, double glass_factor, double shade_factor, double shade) |
| 54 | { |
| 55 | if (Double.isNaN(temp) || Double.isNaN(co2) || Double.isNaN(sun)) |
| 56 | return Double.NaN; |
| 57 | |
| 58 | if (Double.isNaN(glass_factor)) |
| 59 | glass_factor = 0.6; // default: 60% |
| 60 | else |
| 61 | glass_factor /= 100; |
| 62 | |
| 63 | if (Double.isNaN(shade_factor)) |
| 64 | shade_factor = 0.4; // default: 40% |
| 65 | else |
| 66 | shade_factor /= 100; |
| 67 | |
| 68 | if (Double.isNaN(shade)) |
| 69 | shade = 0.5; // default: 50% |
| 70 | else |
| 71 | shade /= 100; |
| 72 | |
| 73 | double light = 2.3 * sun * glass_factor * ((1 - shade) + (shade * shade_factor)); |
| 74 | |
| 75 | double T_l = temp + T_0; |
| 76 | double CO2 = co2; |
| 77 | double I_a = light / 4.6; // divide by 4.6, to convert from [µmol m^-2 s^-1] to [J m^-2 s^-1] |
| 78 | |
| 79 | double P_g; |
| 80 | double P_g_max; |
| 81 | double a_l; |
| 82 | double Gamma; |
| 83 | double P_n_max; |
| 84 | double P_n_CO2; |
| 85 | double J_max; |
| 86 | double rho_CO2_T_l; |
| 87 | double r_CO2; |
| 88 | double V_c_max; |
| 89 | double r_c_CO2; |
| 90 | double R_d; |
| 91 | |
| 92 | double P_n; |
| 93 | |
| 94 | double temp_diff = calculate_temp_diff(T_l); |
| 95 | |
| 96 | Gamma = calculate_Gamma(temp_diff); |
| 97 | rho_CO2_T_l = calculate_rho_CO2_T_l(T_l); |
| 98 | J_max = calculate_J_max(T_l, temp_diff); |
| 99 | 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 | } |