// Copyleft 2012 Chris Korda
// This program is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the Free
// Software Foundation; either version 2 of the License, or any later version.
/*
		chris korda

		rev		date		comments
		00		16oct12		initial version

		convert atmospheric CO2 concentrations to relative temperatures

*/

// This application computes average global temperature changes resulting from
// historical and projected atmospheric CO2 concentrations. Thermal inertia is
// modeled using Climate Response Functions given in the 2011 NASA/GISS paper
// "Earth's energy imbalance and implications" by Hansen et al.
//
// The application takes as input six files containing RCP concentration data,
// in ASCII format. The input files are named RCPfoo_MIDYEAR_CONCENTRATIONS.DAT
// where foo is the scenario name. The files are available from IPCC AR5 RCP
// Scenario data group at the following URL (see the link 'Download all data'):
// http://www.pik-potsdam.de/~mmalte/rcps/index.htm#Download
//
// The application creates a single output file in tab-delimited ASCII format,
// containing four temperature series for each of the six RCP scenarios: one
// temperature series for each of the three Climate Response Functions (slow,
// intermediate, fast) and a fourth series for instantaneous (naive) response.

// CO2ConcToTemps.cpp
//

#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "memory.h"
#include "assert.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

// naive function that converts an atmospheric CO2 concentration (in ppm) to a
// temperature change (in degrees C), assuming instantaneous climate response,
// given a base concentation (in ppm) and a climate sensitivity (in degrees C)
//
double PPM2DTNaive(double Concentration, double BaseConcentration, double Sensitivity)
{
	return(log(Concentration / BaseConcentration) / log(2.0) * Sensitivity);
}

void TestPPM2DTNaive()	// reality check for PPM2DTNaive
{
	double	base = 280;
	double	sens = 3;
	printf("%f\n", PPM2DTNaive(280, base, sens));	// result = 0
	printf("%f\n", PPM2DTNaive(560, base, sens));	// result = 3
	printf("%f\n", PPM2DTNaive(1120, base, sens));	// result = 6
	printf("%f\n", PPM2DTNaive(2240, base, sens));	// result = 9
}

enum {
	FUNC_SEGMENTS = 4	// number of line segments in a climate response function
};

struct RESPONSE_FUNC {	// climate response function
	const char	*Name;	// function name
	double	Pct[FUNC_SEGMENTS];	// final percentage of each function segment
};

enum {	// enumerate climate response functions; must match gRespFunc
	RF_SLOW,
	RF_INTERMEDIATE,
	RF_FAST,
	RF_INSTANTANEOUS,
	RESPONSE_FUNCS,	// total number of climate response functions
};

// Replicate the climate response functions in Hansen et al 2012.
// Each function consists of a polyline in log space. The polyline
// x values are segment final years, defined by gSegmentEndYear;
// the y values are climate response percentages, defined below.
//
static const RESPONSE_FUNC gRespFunc[RESPONSE_FUNCS] = {
//								segment final year 
//	function name			1		10		100		2000
	{"slow",			{	15,		45,		60,		100	}},
	{"intermediate",	{	15,		55,		75,		100	}},
	{"fast",			{	15,		65,		90,		100	}},
	{"instantaneous",	{	100,	100,	100,	100	}},
};

// Define the ending year for each climate response function segment;
// they're all powers of ten EXCEPT the last one, so consequently the
// last segment requires a hard-coded case in CalcResponsePercentage.
//
static const double	gSegmentEndYear[FUNC_SEGMENTS] = {1, 10, 100, 2000};

// Compute the climate response percentage for a given year, by interpolating
// the climate response function specified by RespFuncIdx. The valid range for
// year is [1..2000]. Note that the interpolation is done in log space.
//
double CalcResponsePercentage(int RespFuncIdx, int Year)
{
	assert(RespFuncIdx >= 0 && RespFuncIdx < RESPONSE_FUNCS);
	if (Year <= 0)
		return(0);	// year can't be before end of first segment
	// determine which of the response function's segments contains the given year
	int	iSegment;
	for (iSegment = 0; iSegment < FUNC_SEGMENTS; iSegment++) {
		if (Year <= gSegmentEndYear[iSegment])
			break;
	}
	if (!iSegment)	// if first segment
		return(gRespFunc[RespFuncIdx].Pct[0]);	// no interpolation needed
	if (iSegment >= FUNC_SEGMENTS)	// if beyond last segment
		return(0);	// no response remains to be distributed
	// convert year to a position within segment, normalized between 0 and 1
	double	SegPos;
	if (iSegment < FUNC_SEGMENTS - 1)	// if any segment but last one
		SegPos = log10(double(Year)) - (iSegment - 1);	// generic case for powers of ten
	else	// last segment is special because it doesn't end on a power of ten
		SegPos = (log10(double(Year)) - 2) / (log10(2.0) + 1);	// hard-coded for 2000
	// calculate percentage delta and offset for this segment
	const RESPONSE_FUNC& rc = gRespFunc[RespFuncIdx];
	double	PctDelta = rc.Pct[iSegment] - rc.Pct[iSegment - 1];
	double	PctOffset = rc.Pct[iSegment - 1];
	double	r = SegPos * PctDelta + PctOffset;	// interpolate response percentage
	return(r);
}

struct RESPONSE_SAMPLE {	// climate response percentages for a given year
	double	Pct[RESPONSE_FUNCS];	// one percentage for each climate response function
};

enum {
	RESPONSE_FUNC_YEARS = 2000,	// timespan of a climate response function, in years
};

RESPONSE_SAMPLE	gRespPct[RESPONSE_FUNC_YEARS];	// yearly climate response percentages

// Generate yearly climate response percentages for each of the climate
// response functions, and store them in the gRespPct array.
//
void MakeResponseData()
{
	for (int iYear = 0; iYear < RESPONSE_FUNC_YEARS; iYear++) {
		for (int iFunc = 0; iFunc < RESPONSE_FUNCS; iFunc++) {
			double	Pct = CalcResponsePercentage(iFunc, iYear + 1);
			gRespPct[iYear].Pct[iFunc] = Pct;
		}
	}
}

// Output the yearly climate response percentages to a tab-delimited text file;
// the percentages are generated by MakeResponseData and taken from gRespPct.
//
bool WriteResponseData(const char *Path)
{
	FILE	*fp = fopen(Path, "w");
	if (fp == NULL) {
		printf("can't create %s\n", Path);
		return(false);
	}
	// write header
	fprintf(fp, "year");
	for (int iFunc = 0; iFunc < RESPONSE_FUNCS; iFunc++)
		fprintf(fp, "\t%s", gRespFunc[iFunc].Name);
	fprintf(fp, "\n");
	// write data rows
	for (int iYear = 0; iYear < RESPONSE_FUNC_YEARS; iYear++) {
		fprintf(fp, "%d", iYear + 1);	// write year
		// write percentages for each climate response function
		for (int iFunc = 0; iFunc < RESPONSE_FUNCS; iFunc++)
			fprintf(fp, "\t%.15g", gRespPct[iYear].Pct[iFunc]);
		fprintf(fp, "\n");
	}
	fclose(fp);	// close file
	return(true);
}

enum {
	RCP_SCENARIOS = 6,			// number of RCP scenarios; must match gRCPName
	RCP_FIRST_YEAR = 1765,		// first year of RCP data
	RCP_LAST_YEAR = 2500,		// last year of RCP data
	RCP_YEARS = RCP_LAST_YEAR - RCP_FIRST_YEAR + 1,	// number of years in RCP data
	MODEL_YEARS = RCP_YEARS + RESPONSE_FUNC_YEARS,	// timespan of model, in years
	RCP_HEADER_READS = 34,		// number of reads to skip RCP header
	RCP_DATA_COLS = 35,			// number of RCP data columns
};

static const char *gRCPName[RCP_SCENARIOS] = {	// RCP scenario names
	"3PD",			// RCP3-PD (Peak & Decline)
	"45",			// RCP4.5
	"45SCP45TO3PD",	// RCP4.5 Stablized Concentration Pathway 4.5 to 3PD
	"6",			// RCP6
	"6SCP6TO45",	// RCP6 Stablized Concentration Pathway 6 to 4.5
	"85"			// RCP8.5
};

double	gClimateSensitivity = 3;						// climate sensitivity, in degrees C

// Read yearly absolute CO2 concentrations from an IPCC AR5 RCP concentration
// dataset in ASCII format, and store them in the Concentration array.
//
bool ReadRCPConcentrationData(const char *Path, double Concentration[RCP_YEARS])
{
	FILE	*fp = fopen(Path, "r");
	if (fp == NULL) {
		printf("can't open %s\n", Path);
		return(false);
	}
	// skip RCP data file's header
	for (int iSkip = 0; iSkip < RCP_HEADER_READS; iSkip++)
		fscanf(fp, "\r\n%*[^\r\n]");
	for (int iSamp = 0; iSamp < RCP_YEARS; iSamp++) {
		int	year;
		double	samp;
		// read year and CO2 concentration (first and fourth columns)
		int	retc = fscanf(fp, "%d %*f %*f %lf", &year, &samp);
		// verify that read worked and year is as expected
		if (retc != 2 || year != iSamp + RCP_FIRST_YEAR) {
			printf("parse error\n");
			return(false);
		}
		Concentration[iSamp] = samp;	// store the sample
		for (int iSkip = 0; iSkip < RCP_DATA_COLS - 3; iSkip++)	// skip remaining columns
			fscanf(fp, "%*f");
	}
	return(true);
}

// Compute a series of yearly temperatures given a climate response function,
// based on absolute CO2 concentrations in the Concentration array, and store
// the output in the Temperature array. The unit is average global surface
// temperature, relative to zero at the start of the data, in degrees C.
//
void MakeTemperatureData(int RespFuncIdx, const double Concentration[RCP_YEARS], double Temperature[MODEL_YEARS])
{
	// build table of yearly concentration deltas taking climate response into account
	double	ConcentrationDelta[MODEL_YEARS];			// yearly concentration deltas, in ppm
	memset(ConcentrationDelta, 0, sizeof(ConcentrationDelta));	// zero concentration deltas
	// note iSamp starts at one; first year's concentration delta is zero by default
	for (int iSamp = 1; iSamp < RCP_YEARS; iSamp++) {
		double	ConcDelta = Concentration[iSamp] - Concentration[iSamp - 1];
		double	PrevRespPct = 0;
		// distribute this year's concentration delta over the next two millenia,
		// according to yearly percentages defined by the climate response function
		for (int iYear = 0; iYear < RESPONSE_FUNC_YEARS; iYear++) {
			double	RespPct = gRespPct[iYear].Pct[RespFuncIdx];
			double	ConcDeltaFrac = ConcDelta * (RespPct - PrevRespPct) / 100;
			ConcentrationDelta[iSamp + iYear] += ConcDeltaFrac;
			PrevRespPct = RespPct;
		}
	}
	// now create yearly relative temperatures from the RCP CO2 concentration data
	double	sum = Concentration[0];	// init sum to first year's absolute concentration
	for (int iYear = 0; iYear < MODEL_YEARS; iYear++) {
		sum += ConcentrationDelta[iYear];	// add concentration delta to sum; order matters
		// compute temperature from CO2 concentration given a climate sensitivity
		Temperature[iYear] = PPM2DTNaive(sum, Concentration[0], gClimateSensitivity);
	}
}

double	gTemperature[RCP_SCENARIOS][RESPONSE_FUNCS][MODEL_YEARS];	// temperatures, in degrees C

// Output the gTemperature array to a tab-delimited text file.
//
bool WriteTemperatureData(const char *Path)
{
	FILE	*fp = fopen(Path, "w");
	if (fp == NULL) {
		printf("can't create %s\n", Path);
		return(false);
	}
	// write first header row
	int	iScen;
	for (iScen = 0; iScen < RCP_SCENARIOS; iScen++) {	// for each RCP scenario
		for (int iFunc = 0; iFunc < RESPONSE_FUNCS; iFunc++) {	// for each response function
			fprintf(fp, "\tRCP%s", gRCPName[iScen]);
		}
	}
	fprintf(fp, "\n");
	// write second header row
	fprintf(fp, "year");
	for (iScen = 0; iScen < RCP_SCENARIOS; iScen++) {	// for each RCP scenario
		for (int iFunc = 0; iFunc < RESPONSE_FUNCS; iFunc++) {	// for each response function
			fprintf(fp, "\t%s", gRespFunc[iFunc].Name);
		}
	}
	fprintf(fp, "\n");
	// write data rows
	for (int iYear = 0; iYear < MODEL_YEARS; iYear++) {
		fprintf(fp, "%d", iYear + RCP_FIRST_YEAR);	// write year
		for (int iScen = 0; iScen < RCP_SCENARIOS; iScen++) {	// for each RCP scenario
			for (int iFunc = 0; iFunc < RESPONSE_FUNCS; iFunc++) {	// for each response function
				fprintf(fp, "\t%.15g", gTemperature[iScen][iFunc][iYear]);	// write temperature
			}
		}
		fprintf(fp, "\n");
	}
	fclose(fp);	// close file
	return(true);
}

bool CO2ConcToTemps()
{
	MakeResponseData();	// make climate response data
	if (!WriteResponseData("ClimateResponseFuncs.txt"))	// output climate response data
		return(false);
	for (int iScen = 0; iScen < RCP_SCENARIOS; iScen++) {	// for each RCP scenario
		char	path[256];
		sprintf(path, "RCP%s_MIDYEAR_CONCENTRATIONS.DAT", gRCPName[iScen]);
		double	RCPConcentration[RCP_YEARS];	// RCP absolute CO2 concentration values in ppm
		if (!ReadRCPConcentrationData(path, RCPConcentration))	// read RCP CO2 concentrations
			return(false);
		for (int iFunc = 0; iFunc < RESPONSE_FUNCS; iFunc++) {	// for each response function
			// create temperature data
			MakeTemperatureData(iFunc, RCPConcentration, gTemperature[iScen][iFunc]);
		}
	}
	if (!WriteTemperatureData("RCPTemperatures.txt"))	// output temperature data
		return(false);
	return(true);
}

int main()
{
	//TestPPM2DTNaive();
	return(!CO2ConcToTemps());
}