From f2e450db8fd435139eb9599745e94119fbeeb76c Mon Sep 17 00:00:00 2001
From: Christian Fries <email@christian-fries.de>
Date: Thu, 31 Aug 2023 14:03:12 +0200
Subject: [PATCH] Added old implementation of Emission Intensity Function

---
 .../EmissionIndustrialIntensityFunction.java  | 78 +++++++++++++++++++
 1 file changed, 78 insertions(+)
 create mode 100644 src/main/java/net/finmath/climate/models/dice/submodels/EmissionIndustrialIntensityFunction.java

diff --git a/src/main/java/net/finmath/climate/models/dice/submodels/EmissionIndustrialIntensityFunction.java b/src/main/java/net/finmath/climate/models/dice/submodels/EmissionIndustrialIntensityFunction.java
new file mode 100644
index 0000000000..2273f4d0b4
--- /dev/null
+++ b/src/main/java/net/finmath/climate/models/dice/submodels/EmissionIndustrialIntensityFunction.java
@@ -0,0 +1,78 @@
+package net.finmath.climate.models.dice.submodels;
+
+import java.util.function.BiFunction;
+
+import net.finmath.time.TimeDiscretization;
+
+/**
+ * The function that maps \(i, \sigma(t_{i})) \) to \sigma(t_{i+1})), where \( \sigma(t) \) is the emission intensity (in kgCO2 / USD = GtCO2 / (10^12 USD)).
+ *
+ * The emission intensity is the factor that is applied to the GDP to get the corresponding emissions.
+ *
+ * The function is modelled as an exponential decay, where the decay rate decays exponentially (double exponential).
+ *
+ * Note: This is a continuous time version of the function \( \sigma(t) \) from the original model, except that the division by \( (1-\mu(0)) \) is missing here.
+ *
+ * @author Christian Fries
+ */
+public class EmissionIndustrialIntensityFunction implements BiFunction<Integer, Double, Double> {
+
+	private static double e0 = 35.85;					// Initial emissions
+	private static double q0 = 105.5;					// Initial global output
+	private static double sigma0 = e0/q0;				// Calculated initial emissions intensity, the 1/(1-mu0) is outside
+
+	//	private static double mu0 = 0.03;					// Initial mitigation rate
+	//	private static double sigma0 = e0/(q0*(1-mu0));		// Calculated initial emissions intensity
+
+	private final TimeDiscretization timeDiscretization;
+	private final double emissionIntensityInitial;		// sigma0;
+	private final double emissionIntensityRateInitial;	// = 0.0152;		// -g	// per year
+	private final double emissionIntensityRateDecay;	// exp decay rate corresponding to annual -0.001;			// -d	// per year
+
+	/**
+	 * The evolution of the emission intensity
+	 * @param timeDiscretization The time discretization used.
+	 * @param emissionIntensityInitial The initial emission intensity. Unit: GtCO2 / (10^12 USD)
+	 * @param emissionIntensityRateInitial Initial value for the exponential decay rate.
+	 * @param emissionIntensityRateDecay Exponential decay rate for the decay rate.
+	 */
+	public EmissionIndustrialIntensityFunction(TimeDiscretization timeDiscretization, double emissionIntensityInitial,
+	                                           double emissionIntensityRateInitial, double emissionIntensityRateDecay) {
+		super();
+		this.timeDiscretization = timeDiscretization;
+		this.emissionIntensityInitial = emissionIntensityInitial;
+		this.emissionIntensityRateInitial = emissionIntensityRateInitial;
+		this.emissionIntensityRateDecay = emissionIntensityRateDecay;
+	}
+
+	public EmissionIndustrialIntensityFunction(TimeDiscretization timeDiscretization) {
+		// Parameters from original model
+		this(timeDiscretization, sigma0, 0.0152, -Math.log(1-0.001));
+	}
+
+	@Override
+	public Double apply(Integer timeIndex, Double _emissionIntensity) {
+		double time = timeDiscretization.getTime(timeIndex);
+		final double emissionIntensityRate = emissionIntensityRateInitial * Math.exp(-emissionIntensityRateDecay * time);
+		
+		final double emissionIntensity = emissionIntensityInitial * Math.exp((emissionIntensityRate-emissionIntensityRateInitial) / emissionIntensityRateDecay);
+
+		return emissionIntensity;
+	}
+
+	public TimeDiscretization getTimeDiscretization() {
+		return timeDiscretization;
+	}
+
+	public double getEmissionIntensityInitial() {
+		return emissionIntensityInitial;
+	}
+
+	public double getEmissionIntensityRateInitial() {
+		return emissionIntensityRateInitial;
+	}
+
+	public double getEmissionIntensityRateDecay() {
+		return emissionIntensityRateDecay;
+	}
+}