// G. Hagopian doing the SIR modeling

// Pseudocode to produce data for the plot
// precondition: initial values s(0), i(0), r(0) and 
//	             parameters b and k.
// postcondition: Three vector<double> objects: susceptibles, infectious 
//             and recovered that are produce by the iterative processes
//
//                  s(n) = s(n-1) + change in s
//                  i(n) = i(n-1) + change in i
//                  r(n) = r(n-1) + change in r
//
//             following the model described above.

#include <iostream>
#include <vector>
#include <fstream>
#include <iomanip>
using namespace std;

void print2file(ofstream& , vector<double>, vector<double>, vector<double>);

int main() {
	ofstream ofs("SIR_data.dat");
	vector<double> susceptibles, infectious, recovered;
	double b{}, k{};
	//NYC parameters s(0) = 1, i(0)≈1.27e−6, r(0) = 0
	b = 0.5;
	k = 1.0 / 3.0;
	// vectors to contain population proportions
	susceptibles.push_back(1.);
	infectious.push_back(1.27e-6);
	recovered.push_back(0);
	double delta_s{}, delta_i{}, delta_r{};
	//∆s=−b·s(t)·i(t)∆t  
	// note: we take ∆t = 1 day
	delta_s = -b * susceptibles[0] * infectious[0];
	//∆r(t)=k·i(t)∆t
	delta_r = k * infectious[0];
	//∆i=(b·s(t)i(t)−k·i(t))∆t
	delta_i = -(delta_r + delta_s);
	susceptibles.push_back(susceptibles[0] + delta_s);
	infectious.push_back(infectious[0] + delta_i);
	recovered.push_back(recovered[0] + delta_r);
	// insert loop here...set flag when i starts decreasing, and double days from there
	print2file(ofs, susceptibles, infectious, recovered);
}


void print2file(ofstream& ofs, vector<double> s, vector<double> i, vector<double> r) {
	ofs << "Susceptibles\tInfections\tRecovered\n";
	for (int j = 0; j < s.size(); ++j) {
		ofs << s[j] << ", " << i[j] << ", " << r[j] << endl;
	}

}