//
//  BEM.cpp
//  
//
//  Created by Arne Pommerening on 01.02.2014.
//  Copyright 2014 Philodendron International. All rights reserved.
//

#include <Rcpp.h>
#include <cmath>
#include <ctime>
using namespace Rcpp;
using namespace std;


/**
 * This function calculates the age transformation.
 */
double transformAge(int age, double trans) {
    return (age - trans) / trans;
}

/**
 * This function transforms a relative growth variable to a multiplier.
 */
double relGrowthToMultiplier(double p) {
    return 1 / (1 - p);
}

/**
 * This uses the Wenk function to calculate the relative increment for 10 years.
 */
double calcIncrementTenYears(double c1, double c2, double c3, double c4, int tage, double trans) {
    double t_trans = transformAge(tage, trans);
    return exp(-c1 / (c4 + 1) * exp((c4 + 1) * log(t_trans)) * (1 - exp(-c2 * t_trans * (1 - exp(-c3 * t_trans))))); // Wenk's original
}


/**
 * This function calculates the annual multiplier based on the Gerold & Roemisch paper.
 */
double calcAnnualMultiplier(double c1, double c2, double c3, double c4, int age, double trans) {
    double px;
    double product[2], mvx[2];
    // 0 und 1 im array legen einfach nur fest, was Zaehler und Nenner ist.
    for (int i = 0; i < 2; i++)
        product[i] = 1;
    int tage = age;
    do {
        px = calcIncrementTenYears(c1, c2, c3, c4, tage + 9, trans);
        mvx[0] = relGrowthToMultiplier(px);
        product[0] *= mvx[0];
        px = calcIncrementTenYears(c1, c2, c3, c4, tage + 10, trans);
        mvx[1] = relGrowthToMultiplier(px);
        product[1] *= mvx[1];
        tage += 10;
    }
    while (mvx[0] / mvx[1] > 1+1E-16);
    return product[0] / product[1];
}

/**
 * This function calculates the final multiplier for a given time period (interval), age1 is really the initial age, not the final age.
 */
// [[Rcpp::export]]
double interpolateMultiplier(int age1, double trans, int interval, double c1, double c2, double c3, double c4) {
	double mvy = 1.0;
	int age = age1;
	for (int i = 0; i < interval; i++) {
		age++;
		mvy *= calcAnnualMultiplier(c1, c2, c3, c4, age, trans);
	}
	return mvy;
}