//*** The code for the simulations of the article 
//*** "Antiangiogenic Therapy Efficacy Can Be Tumor Size Dependent, as Mathematical Modeling Suggests"
//*** by Maxim Kuznetsov and Andrey Kolobov
//*** submitted to Mathematics
//*** 2023
//*** This code works fine in Dev-C++, version 5.11
//*** Note, that some small parts of the code may need to be revised depending on the utilized software 

//***********************************************************************************************************************************************************
//********************************* INTRODUCTION ************************************************************************************************************
//***********************************************************************************************************************************************************

//********* Initializing needful things

#include <stdio.h>
#include <tchar.h>
#include <windows.h>
#include <ctime>
#include <iostream>
#include <cmath>
#include <conio.h>
#include <fstream>

using namespace std;

#define output_frequency 1. // how often the distributions of varibles will be recorded

//********* MODEL PARAMETERS
double xi = 5.; // radius of chemoterapeutic agent molecule
double chi = 0.05; //1./xi/xi/xi; // sensitivity to chemotherapeutic agent
//Chemotherapy
int AATonCTdose = 5; // on which dose of CT AAT is performed
int sizeCT = 6; // how many injections are in the scheme
double RadCTbegin = 4.; // radius of tumor, at whih chemotherapy starts (in mm)
double R0T = 40.; // initial radius of tissue region
#define inception_width 2. // the radius of inital tumor cell population


// Cells
double B = 0.01; // maximum rate of tumor cell proliferation
double sp = 15.; // critical stress for tumor cell proliferation
double eps = 500.; // smoothing parameter for Heaviside functions
double M = 0.01; //  maximum degradation rate of deadcells
double nu = 0.003; // rate of death by starvation
double gd = 0.001; // critical level of glucose for tumor cell survival

// Stress
double k = 500.; // solid stress coefficient 
double ss = 0.3; // minimum fraction of interacting cells
double s0 = 0.8; // normal fraction of cells

// Interstitial fluid 
double Ln = 0.1; // hydraulic conductivity of normal capillaries
double La = 0.22; // hydraulic conductivity of abnormal capillaries
double pc = 4.; // fluid pressure in capillaries
double Ll = 1300.; // hydraulic conductivity of lymphatic capillaries
double K = 0.1; // tissue hydraulic conductivity

// VEGF
double Sv = 1.; // secretion rate
double omega = 1.; // internalization rate
double Mv = 0.01; // degradation rate
double Dv = 21.; // diffusion coefficient

// Capillaries
double R = 0.008; // maximum rate of angiogenesis
double cmax = 5.; // maximum surface area density
double Mc = 0.03; // microvasculature degradation rate
double kM = 2.; // coefficient of degradation in the tumor core
double Vn = 0.1; // normalization rate
double Vd = 0.1; // denormalization rate
double mu = 0.002; // pruning rate
double vst = 0.001; // Michaelis constant for VEGF action
double Dc = 0.03; // coefficient of active motion

// Glucose
double gst = 0.01; // Michaelis constant for consumption rate
double Pgn = 4.; // permeability of normal capillaries
double Pga = 10.; // permeability of abnormal capillaries
double nug = 1200.; // coefficient of proliferating tumor cells consumption rate
double Qgh = 0.5; // glucose consumption rate by normal cells 
double Dg = 10.; // glucose diffusion coefficient

//Chemotherapeutic agent


double Du = 13.; // diffusion coefficient
double Cu = 0.0015; // clearance rate

double sun = 0.09; // fraction of available pore cross-section area of normal capillaries
double sua = 0.58; // fraction of available pore cross-section area of abnormal capillaries
double Pun = 0.007; // permeability of normal capillaries
double Pua = 0.25;  // permeability of abnormal capillaries

#define size 5000. // the radius of a simulation region, in 10^(-2) cm

# define PI 3.14159265358979323846  // pi

//********* SIMULATION PARAMETERS

int ddebug = 99999;

double tau = 0.00001; // initial time step in full model -- will be adjusted 
double tauAdjust;
double Is0NonMonIndexMin = 1.e-7; // lower value of this parameter...
double tauAdjustParameter = 50; // and greater value of this parameter provide more accurate calculations
#define h 0.05 // space step
double tauMax = tau;
double tauMin = tau;

#define N (int(size/h)) // number of points in space, minus one
#define tau_adjust_frequency 99999.//0.1 // how often time step will be adjusted
#define rsn_frequency 1. // how often the tumor radius, speed of tumor front and number of tumor cells will be recorded
FILE* results;
#define inception_width_int (int)(inception_width/h) // the size, in cell grids, of inital tumor cell population

//********* COUNTERS

double t = 0.; // main
double ta = 0.; // for adjustment of time step
double ti = 0.; // for output of the distributions of varibles
double tr = 0.; // for output of the tumor radius, speed of tumor front and number of tumor cells
int d = 0; // for serial numbers of the files with distributions of varibles

//********* MODEL VARIABLES

double np0[N + 1], npr[N + 1], np1[N + 1]; // proliferating tumor cells
double nq0[N + 1], nqr[N + 1], nq1[N + 1]; // quiescent tumor cells
double h0[N + 1], hr[N + 1], h1[N + 1]; // normal cells
double m0[N + 1], mr[N + 1], m1[N + 1]; // damaged cells
double v0[N + 1], vr[N + 1], v1[N + 1]; // VEGF
double V0[N + 1], VR[N + 1], V1[N + 1]; // VEGF -- for easy solution in spherically-symmetrical case
double cn0[N + 1], cnr[N + 1], cn1[N + 1]; // normal capillaries
double ca0[N + 1], car[N + 1], cac[N + 1], ca1[N + 1]; // abnormal capillaries
double CA0[N + 1], CAR[N + 1], CA1[N + 1]; // abnormal capillaries -- for easy solution in spherically-symmetrical case
double g0[N + 1], gr[N + 1], g1[N + 1]; // glucose
double G0[N + 1], GR[N + 1], G1[N + 1]; // glucose -- for easy solution in spherically-symmetrical case
double nh0[N + 1], nh1[N + 1]; // all cells -- for convenient solution of convection equation
double u0[N + 1], ur[N + 1], uc[N + 1], u1[N + 1]; // chemotherapeutic agent
double U0[N + 1], UR[N + 1], U1[N + 1]; // chemotherapeutic agent -- for easy solution in spherically-symmetrical case
double ubl; // chemotherapeutic agent in blood
double sigma[N + 1]; // stress

//********* PROCEDURES AND FUNCTIONS

void Begin(); // complex of actions for the beginning, ...
void Cycle_Ending(); // ...ending of main cycle, ...
void End(); // ...and ending of the program execution

void TauAdjust(); // for adjustment of time step

void Convection(); // solving the equation for cells convection
void ConvectionCA(); // solving the equation for chemotherapeutic agent convection
void Stress(double* nh); // computing stress
double StressSingle(double nh);
void ConvSpeed(double* ccn, double* cca, double* hh, double* si, double* II, double* IIf, double* nh); // computing convective speed
void CT(); // Chemotherapy

void Diffusion_glucose(); //  solving the equation for diffusion of glucose
void Diffusion_VEGF(); //  solving the equation for diffusion of VEGF
void Diffusion_cap(); //  solving the equation for diffusion of abnormal capillaries
void Diffusion_CA(); //  solving the equation for diffusion of chemotherapeutic agent

void TMA(double apr, double* bpr, double cpr, double* F, double* prog,
	double muleft, double hileft, double muright, double hiright); // tridiagonal matrix algorithm

void TMAc(double apr, double bpr, double cpr, double* F, double* prog,
	double muleft, double hileft, double muright, double hiright);

void Euler(); // Euler method
double Dq[N + 1], MMg[N + 1], MMv[N + 1], MMvst[N + 1], TG[N + 1], Pr[N + 1], kI[N + 1], kIf[N + 1], BUCA[N + 1]; // auxiliary variables for kinetic terms

double Tp(double ss); // function that defines rate of tumor cells proliferation depending on stress
double Ttr(double gg); // function that defines transitions of tumor cells to proliferation and back depending on the level of glucose
double Td(double gg); // function that defines rate of death of cells depending on the level of glucose
#define Nm 10000 // number of precalculated values of theta functions
double TtrPre[Nm + 1]; // values for Ttr(gg) are precalculated for optimization purposes
double Ttrg[N + 1];

void Profile_output(); // output of distributions of variables
void Radius_Speed_CellNumber_output(); // output of the tumor radius, whole tissue radius, speed of tumor front and number of tumor cells
void Radius_estimation(); // estimation of tumor radius and whole tissue radius
void Speed_estimation(); // estimation of tumor front speed

//********* CHEMOTHERAPY

double** CTscheme; // current CT injection schedule
int doseCTcnt = 1; // what fraction is to be done now
int CTisOn = 0; // if it is 1, then tumor has reached the radius, at which chemotherapy starts 
double t_preCT; // to remember the value before chemotherapy

//********* AUXILIARY VARIABLES

FILE* inputInitial;
double lenmmprev, lenmm, npIn, nqIn, hIn, mIn, vIn, cnIn, caIn, gIn, IIn, sigmaIn, uIn;

// smoothen speed field and h
double ksm = 0.3;
int ci = 0;
double I0sm[N + 2], hsm[N + 2];

////// for glucose diffusion
double kg = 2 * h * h / (tau * Dg);
double aprg = 1.;
double bprg = 2. + kg;
double cprg = 1.;
double muleftg, hileftg, murightg, hirightg;
double Fg[N + 1];
double progg[N + 1];
double pp[N + 1], qq[N + 1];

////// for VEGF diffusion
double kv = 2 * h * h / (tau * Dv);
double aprv = 1.;
double bprv = 2. + kv;
double cprv = 1.;
double muleftv, hileftv, murightv, hirightv;
double Fv[N + 1];
double progv[N + 1];

////// for capillaries diffusion
double kc = 2 * h * h / (tau * Dc);
double aprc = 1.;
double bprc = 2. + kc;
double cprc = 1.;
double muleftc, hileftc, murightc, hirightc;
double Fc[N + 1];
double progc[N + 1];

////// for chemotherapeutic agent diffusion
double ku = 2 * h * h / (tau * Du);
double apru = 1.;
double bpru = 2. + ku;
double cpru = 1.;
double muleftu, hileftu, murightu, hirightu;
double Fu[N + 1];
double progu[N + 1];

////// for convection
double Sgn(double x);
double Is0[N + 2], Is0_not_corr[N + 3];
double If0[N + 2], If0_not_corr[N + 3];
double Is0Max;
double Is0NonMonIndex;
double Q_plus[N + 2], Q_minus[N + 2];
double delta_nh[N + 2], delta_np[N + 2], delta_nq[N + 2], delta_m[N + 2], delta_cn[N + 2], delta_ca[N + 2];
double delta_uf[N + 2], delta_ub[N + 2];
double one_eighth[N + 2];
double f_c_nh[N + 2], f_c_np[N + 2], f_c_nq[N + 2], f_c_m[N + 2], f_c_cn[N + 2], f_c_ca[N + 2];
double f_c_uf[N + 2], f_c_ub[N + 2];

int NN = inception_width_int; // number of point in space where tumor cells are
int NH = (int(R0T / h)); // number of point in space where tumor and normal cells are
double hNN = 0.; // the variable size of the last grid element, NN, where tumor cells are
double hNH = 0.; // the variable size of the last grid element, NH, where normal cells are
double hNNnew;
double VNNnew;
double hNHnew;
double VNHnew;
double knNN, khNN; // coefficient of volume of NN-th grid element related to tumor and normal tissue

////// for tumor radius, speed of tumor front and number of tumor cells
double RadTis; // current whole tissue radius, in mm
double Rad; // current tumor radius, in mm
double Rad_prev; // previously measured tumor radius, in mm
double T_prev; // time of previous measurement of tumor radius
double V; // tumor front speed
double NTC, NTCsurv, NTCr, NTCpr, NTCcm1, NTCpm1; // relative number of tumor cells
double WDnoNP, WDNP, BEDh, WDNPprop; // total fraction of cells that would die with nanoparticels, without them, and normal tissue damage 
double WDnoNPopt, WDNPopt; // total fraction of cells that would die under optimized distiburtion of radiation with nanoparticels, without them
double HTC, HTCr, HTCcm1, HTCrLIN, HTCcm1LIN; // relative number of normal cells -- for solution of convective equation
double PIT; // relative number of particles in tumor
double T_est_fin = 0.; // auxiliary variable for program stop
double V_prev; // auxiliary variable for program stop

////// for work with files for distribution of variables
char File_name1[80], num[10], radius[8];

// volumes of spatial grids
double Vol[N + 1], Vact[N + 1];

// for spatial optimization
double EnTot;
double Dact[N + 1];
double OERa[N + 1];
double OERb[N + 1];
double EffGrid[N + 1], EffGridTot[N + 1];
int GridChart[N + 1], GridChartRev[N + 1];
int Stop, jj, xxx;
double aux;

double hstat = s0;

//***********************************************************************************************************************************************************
//********************************* MAIN CYCLE **************************************************************************************************************
//***********************************************************************************************************************************************************

int main()
{
	char buffer[MAX_PATH];
	GetModuleFileNameA(NULL, buffer, MAX_PATH);
	string str(buffer);


	Begin();

	// the program runs until tumor grows no more, that is checked in Speed_estimation()
	while (1 == 1)
	{
		t = t + tau;
		Radius_estimation();
		Speed_estimation();

		if (CTisOn == 1 && t - t_preCT - CTscheme[0][doseCTcnt] >= 0 && doseCTcnt <= sizeCT) // injection
		{
			if(AATonCTdose<=doseCTcnt)
				for (int i = 0; i <= NH; i++)
				{
					v0[i] = 0.;
					Sv=0.;
				}
	
			ubl = ubl + CTscheme[1][doseCTcnt];
			doseCTcnt++;
		}

		Euler();
		Convection();

		if (CTisOn == 1)
		{
			ConvectionCA();
			Diffusion_CA();
		}

		Diffusion_glucose();

		Diffusion_VEGF();
		Diffusion_cap();

		if (t >= tr - 0.00001 || d > ddebug)
			Radius_Speed_CellNumber_output();
		if (t >= ti - 0.00001 || d > ddebug) // a trifle is subtracted to overcome the inevitable minor machine errors during arithmetic operations
			Profile_output();

		if (t >= ta - 0.00001)
			TauAdjust();

		Cycle_Ending();

		if (CTisOn == 0 && Rad > RadCTbegin)
		{
			CTisOn = 1;
			t_preCT = t;
		}
	}
}

//***********************************************************************************************************************************************************
//********************************* BEGINNING OF SIMULATION *************************************************************************************************
//***********************************************************************************************************************************************************

void Begin()
{
	// standard "volumes" of grid elements
	for (int i = 0; i <= N; i++)
		Vol[i] = (i + 1)*h*(i + 1)*h*(i + 1)*h - i*h*i*h*i*h;
		
	//********* PRECALCULATING Theta
	for (int i = 0; i <= Nm; i++)
		TtrPre[i] = Ttr(1. * i / Nm);

	//********* CREATING FILE FOR TUMOR RADIUS, FRONT SPEED AND CELL NUMBER
	if ((results = fopen("results.txt", "wt")) == NULL)
		fprintf(stderr, "Cannot open file results.txt\n\n");

	fprintf(results, "t(h)	Rad(mm)	Speed 	TCells 		RTi(mm)	ubl\n");
	fflush(results);

	//********* INITIAL CONDITIONS
	
	for (int i = 1; i <= 1000; i++)
		hstat = hstat - 0.00001 * (Ln * (pc + StressSingle(hstat)) + Ll * StressSingle(hstat));

	if ((inputInitial = fopen("Initial.txt", "rt")) != NULL)
	{
		cout << "\nReading initial conditions \n";
		NN = 0;
		NH = 0;

		npIn = 1.;
		char* word;
		int iScan;
		word = new char[512];
		fgets(word, 512, inputInitial);
		fgets(word, 512, inputInitial);
		do {
			fgets(word, 512, inputInitial);
			lenmmprev = lenmm;
			iScan = sscanf(word, "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", &lenmm, &npIn, &nqIn, &hIn, &mIn, &vIn, &cnIn, &caIn, &gIn, &IIn, &sigmaIn, &uIn);
			np0[NN] = npIn;
			nq0[NN] = nqIn;
			h0[NN] = hIn;
			m0[NN] = mIn;
			v0[NN] = vIn;
			cn0[NN] = cnIn;
			ca0[NN] = caIn;
			g0[NN] = gIn;
			u0[NN] = uIn;
			NN++;
		} while (npIn + nqIn > 0.);

		NN--;
		NN--;
		hNN = 10. * lenmmprev - (1. * NN) * h;
		NH = NN + 1;
		do {
			fgets(word, 512, inputInitial);
			lenmmprev = lenmm;
			iScan = sscanf(word, "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", &lenmm, &npIn, &nqIn, &hIn, &mIn, &vIn, &cnIn, &caIn, &gIn, &IIn, &sigmaIn, &uIn);
			np0[NH] = npIn;
			nq0[NH] = nqIn;
			h0[NH] = hIn;
			m0[NH] = mIn;
			v0[NH] = vIn;
			cn0[NH] = cnIn;
			ca0[NH] = caIn;
			g0[NH] = gIn;
			u0[NH] = uIn;
			NH++;
		} while (!feof(inputInitial));

		NH--;
		NH--;
		hNH = 10. * lenmmprev - (1. * NH) * h;
		fclose(inputInitial);
	}
	else
	{
		for (int i = 0; i <= NH; i++)
		{
			np0[i] = 0.;
			nq0[i] = 0.;
			h0[i] = hstat;
			m0[i] = 0.;
			v0[i] = 0.;
			cn0[i] = 1.;
			ca0[i] = 0.;
			g0[i] = 1.;
			u0[i] = 0.;
		}

		for (int i = 0; i <= NN - 1; i++)
		{
			np0[i] = hstat;
			h0[i] = 0.;
			cn0[i] = 0.;
		}

		np0[NN] = hstat;
		h0[NN] = hstat;
	}
	
	for (int i = 0; i <= N; i++) // this is made for easy output of the initial conditions
	{
		np1[i] = np0[i];
		nq1[i] = nq0[i];
		h1[i] = h0[i];
		m1[i] = m0[i];
		v1[i] = v0[i];
		cn1[i] = cn0[i];
		ca1[i] = ca0[i];
		g1[i] = g0[i];
		u1[i] = u0[i];
	}

	Radius_estimation();

	CTscheme = new double* [2]; // 0 -- time of injection, 1 -- dose

	for (int i = 0; i <= 1; i++)
		CTscheme[i] = new double[sizeCT + 1];

	CTscheme[0][0] = 0.; // they will be ignored
	CTscheme[1][0] = 0.;

	for (int i = 1; i <= sizeCT; i++)
	{
		CTscheme[0][i] = 24.*7.*3. * (i - 1)+1.;
		CTscheme[1][i] = 1.;
	}

	//********* COUNTING TUMOR CELLS
	NTC = 0.;

	for (int i = 0; i <= NN - 1; i++)
		NTC = NTC + Vol[i] * (np1[i] + nq1[i]);
	//hNN has zero width yet

	NTC = 4. * PI * 300. * NTC / 3.;

	Radius_Speed_CellNumber_output();
	Profile_output();
}

//***********************************************************************************************************************************************************
//********************************* OUTPUT OF DISTRIBUTIONS OF VARIABLES ************************************************************************************
//***********************************************************************************************************************************************************

void Profile_output()
{
	//********* CREATING A FILENAME WITH SERIAL NUMBER

	_itoa(d, num, 10);

	if (d < 10)
		strcpy(File_name1, "res-d0000");
	else if (d < 100)
		strcpy(File_name1, "res-d000");
	else if (d < 1000)
		strcpy(File_name1, "res-d00");
	else if (d < 10000)
		strcpy(File_name1, "res-d0");
	else
		strcpy(File_name1, "res-d");

	strcat(File_name1, num);

	if (RadTis < 1.)
	{
		_gcvt(RadTis, 1, radius);
		if (radius[1] == '\0')
			radius[1] = '0';
	}
	else if (RadTis < 10.)
	{
		_gcvt(RadTis, 2, radius);
		if (radius[2] == '\0')
			radius[2] = '0';
	}
	else
	{
		_gcvt(RadTis, 3, radius);// in mm
		if (radius[3] == '\0')
			radius[3] = '0';
	}

	strcat(File_name1, "_RadTis(mm)=");

	strcat(File_name1, radius);

	strcat(File_name1, ".txt");

	//********* PRINTING THE SERIAL NUMBER, RADIUS, SPEED OF TUMOR FRONT AND NUMBER OF TUMOR CELLS ON A SCREEN
	printf("#%i, radius = %.3f mm, speed = %.3f mm/week, tissue radius = %.3f mm, cell number = %.3e\n", d, Rad, V, RadTis, NTC);

	//********* OUTPUT OF DISTRIBUTIONS OF VARIABLES INTO A FILE
	FILE* out;
	if ((out = fopen(File_name1, "wt")) == NULL)
		fprintf(stderr, "Cannot open file\n");

	fprintf(out, "Len(mm)\t\tnp\t\tnq\t\th\t\tm\t\tv\t\tcn\t\tca\t\tg\t\tI\t\tsigma\t\tu\n\n");

	for (int i = 0; i <= NN - 1; i++)
		fprintf(out, "%.6f\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.5f\t%.7e\n", (1. * i + 0.5) * h / 10., np1[i], nq1[i], h1[i], m1[i], v1[i], cn1[i], ca1[i], g1[i], Is0[i], sigma[i], u1[i]);

		fprintf(out, "%.6f\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.5f\t%.7e\n", ((1. * NN) * h + hNN) / 10., np1[NN], nq1[NN], h1[NN], m1[NN], v1[NN - 1] + (v1[NN + 1] - v1[NN - 1]) * ((h / 2.) + hNN) / (2. * h), cn1[NN - 1] + (cn1[NN + 1] - cn1[NN - 1]) * ((h / 2.) + hNN) / (2. * h), ca1[NN - 1] + (ca1[NN + 1] - ca1[NN - 1]) * ((h / 2.) + hNN) / (2. * h), g1[NN - 1] + (g1[NN + 1] - g1[NN - 1]) * ((h / 2.) + hNN) / (2. * h), Is0[NN], sigma[NN], u1[NN]);

	for (int i = NN + 1; i <= NH - 1; i++)
		fprintf(out, "%.6f\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.5f\t%.7e\n", (1. * i + 0.5) * h / 10., np1[i], nq1[i], h1[i], m1[i], v1[i], cn1[i], ca1[i], g1[i], Is0[i], sigma[i], u1[i]);

		fprintf(out, "%.6f\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.5f\t%.7e\n", ((1. * NH) * h + hNH) / 10., np1[NH], nq1[NH], h1[NH], m1[NH], v1[NH], cn1[NH], ca1[NH], g1[NH], Is0[NH], sigma[NH], u1[NH]);

	fclose(out);

	//********* SWITCHING THE COUNTER FOR OUTPUT OF DISTRIBUTION OF VARIABLES 
	ti = ti + output_frequency;
	d++;
}

//***********************************************************************************************************************************************************
//********************************* ESTIMATION OF RADIUS AND SPEED ******************************************************************************************
//***********************************************************************************************************************************************************

//********* ESTIMATION OF RADIUS

void Radius_estimation()
{
	Rad = (NN * h + hNN) / 10.;
	RadTis = (NH * h + hNH) / 10.;
}

//********* ESTIMATION OF TUMOR FRONT SPEED

void Speed_estimation()
{
	if (t - T_prev > 5. || (abs(Rad - Rad_prev) > 0.2 && 10. * Rad > R0T))
	{
		V = 7. * 24. * (Rad - Rad_prev) / (t - T_prev); // in mm/week
		T_prev = t;
		Rad_prev = Rad;
	}

	if (t - T_est_fin > 10. * 24.)
	{
	/*	if (abs(V) < 1.e-5 && abs(V_prev) < 1.e-5) // growth stop
		{
			cout << "\n\nEnding.";
			_getch();
			End();
		}
		else
			T_est_fin = t;
*/
		V_prev = V;
	}
}

//***********************************************************************************************************************************************************
//********************************* OUTPUT OF THE TUMOR RADIUS, SPEED OF TUMOR FRONT AND NUMBER OF TUMOR CELLS **********************************************
//***********************************************************************************************************************************************************

void Radius_Speed_CellNumber_output()
{
	tr = tr + rsn_frequency;

	//********* COUNTING TUMOR CELLS

	NTC = 0.;

	for (int i = 0; i <= NN - 1; i++)
		NTC = NTC + Vol[i] * (np1[i] + nq1[i]);

	NTC = NTC + ((NN*h + hNN)*(NN*h + hNN)*(NN*h + hNN) - NN*h*NN*h*NN*h) * (np1[NN] + nq1[NN]);

	NTC = 4. * PI * 300. * NTC / 3.;
	
	//********* OUTPUT
	fprintf(results, "%.1f\t%.5f\t%.5f\t%.3e\t%.5f\t%.5f\n", t, Rad, V, NTC, RadTis, ubl);
	fflush(results);
}

//***********************************************************************************************************************************************************
//********************************* SOLUTION OF KINETIC TERMS ***********************************************************************************************
//***********************************************************************************************************************************************************

//********* function that defines rate of tumor cells proliferation depending on stress
double Tp(double ss)
{
	return (0.5 * (1. + tanh(eps * (sp - ss))));
}

//********* function that defines transitions of tumor cells to proliferation and back depending on the level of glucose
double Ttr(double gg)
{
	return (0.5 * (1. + tanh(eps * (gg - gst))));
}


//********* function that defines death rate depending on the level of glucose
double Td(double gg)
{
	return (0.5 * (1. + tanh(eps * (gd - gg))));
}

//********* EULER PROCEDURE

void Euler()
{
	// redefinition of variables for convenience of computing convective speed
	for (int i = 0; i <= NN - 1; i++)
		nh0[i] = np0[i] + nq0[i] + m0[i];

	nh0[NN] = ((np0[NN] + nq0[NN] + m0[NN]) * ((NN * h + hNN) * (NN * h + hNN) * (NN * h + hNN) - (NN * h) * (NN * h) * (NN * h)) + h0[NN] * (((NN + 1) * h) * ((NN + 1) * h) * ((NN + 1) * h) - (NN * h + hNN) * (NN * h + hNN) * (NN * h + hNN))) / (((NN + 1) * h) * ((NN + 1) * h) * ((NN + 1) * h) - (NN * h) * (NN * h) * (NN * h));

	for (int i = NN + 1; i <= NH; i++)
		nh0[i] = h0[i];

	Stress(nh0);
	ConvSpeed(cn0, ca0, h0, sigma, Is0, If0, nh0);

	for (int i = 0; i <= NN; i++)
		Ttrg[i] = TtrPre[(int)(g0[i] * Nm + 0.5)];

	// reaction terms
	for (int i = 0; i <= NH; i++)
	{
		MMv[i] = v0[i] / (v0[i] + vst);
		MMvst[i] = vst / (v0[i] + vst);
		kI[i] = 2. * Is0[i] / ((1. * i + 0.5) * h);
	}

	for (int i = 0; i <= NN - 1; i++)
	{
		MMg[i] = g0[i] / (g0[i] + gst);
		TG[i] = B * (1 - Ttrg[i]) * np0[i] - B * Ttrg[i] * nq0[i];
		Pr[i] = B * np0[i] * Tp(sigma[i]);

		npr[i] = np0[i] + tau * (Pr[i] * MMg[i] - TG[i] - kI[i] * np0[i] - chi * u0[i]*np0[i]);
		nqr[i] = nq0[i] + tau * (TG[i] -  nu * nq0[i]* Td(g0[i]) - kI[i] * nq0[i]);
		mr[i] = m0[i] + tau * (- M * m0[i] - kI[i] * m0[i] + chi * u0[i]*np0[i] + nu * nq0[i]* Td(g0[i]));
		vr[i] = v0[i] + tau * (Sv * nq0[i] - omega * (cn0[i] + ca0[i]) * v0[i] - Mv * v0[i]);
		cnr[i] = cn0[i] + tau * (-Mc * kM *(nq0[i] + m0[i]) * cn0[i] + Vn * MMvst[i] * ca0[i] - Vd * MMv[i] * cn0[i] - kI[i] * cn0[i]);
		if (cn0[i] > 1.)
			cnr[i] = cnr[i] - tau * (mu * (cn0[i] - 1.));
		car[i] = ca0[i] + tau * (-Mc * (np0[i] + kM * (nq0[i] + m0[i])) * ca0[i] + R * MMv[i] * (cn0[i] + ca0[i]) * (1. - (cn0[i] + ca0[i]) / cmax) - Vn * MMvst[i] * ca0[i] + Vd * MMv[i] * cn0[i] - kI[i] * ca0[i]);
		gr[i] = g0[i] + tau * ((Pgn * cn0[i] + Pga * ca0[i]) * (1. - g0[i]) - (nug * Pr[i] + Qgh * (nq0[i] + np0[i] * (1. - Tp(sigma[i])))) * MMg[i]);
	}

	knNN = ((NN * h + hNN) * (NN * h + hNN) * (NN * h + hNN) - (NN * h) * (NN * h) * (NN * h)) / (((NN + 1) * h) * ((NN + 1) * h) * ((NN + 1) * h) - (NN * h) * (NN * h) * (NN * h));
	khNN = (((NN + 1) * h) * ((NN + 1) * h) * ((NN + 1) * h) - (NN * h + hNN) * (NN * h + hNN) * (NN * h + hNN)) / (((NN + 1) * h) * ((NN + 1) * h) * ((NN + 1) * h) - (NN * h) * (NN * h) * (NN * h));

	npr[NN] = np0[NN] + tau * (Pr[NN] * MMg[NN] - TG[NN] - kI[NN] * np0[NN] - chi * u0[NN]*np0[NN]);
	nqr[NN] = nq0[NN] + tau * (TG[NN] -  nu * nq0[NN]* Td(g0[NN]) - kI[NN] * nq0[NN]);
	hr[NN] = h0[NN] - tau * (kI[NN] * h0[NN]);
	mr[NN] = m0[NN] + tau * (-M * m0[NN] - kI[NN] * m0[NN] + chi * u0[NN]*np0[NN] + nu * nq0[NN]* Td(g0[NN]));
	vr[NN] = v0[NN] + tau * (-omega * (cn0[NN] + ca0[NN]) * v0[NN] - Mv * v0[NN]);

	cnr[NN] = cn0[NN] + tau * (-Mc * kM *(nq0[NN] + m0[NN]) * knNN * cn0[NN] + Vn * MMvst[NN] * ca0[NN] - Vd * MMv[NN] * cn0[NN] - kI[NN] * cn0[NN]);
	if (cn0[NN] > 1.)
		cnr[NN] = cnr[NN] - tau * (mu * (cn0[NN] - 1.));
	car[NN] = ca0[NN] + tau * (-Mc * (np0[NN] + kM * (nq0[NN] + m0[NN])) * knNN * ca0[NN] + R * MMv[NN] * (cn0[NN] + ca0[NN]) * (1. - (cn0[NN] + ca0[NN]) / cmax) - Vn * MMvst[NN] * ca0[NN] + Vd * MMv[NN] * cn0[NN] - kI[NN] * ca0[NN]);

	gr[NN] = g0[NN] + tau * ((Pgn * cn0[NN] + Pga * ca0[NN]) * khNN * (1. - g0[NN]) - (nug * Pr[NN] * knNN + Qgh * (h0[NN] * khNN + nq0[NN] * knNN + np0[NN] * knNN * (1. - Tp(sigma[NN])))) * MMg[NN]);

	for (int i = NN + 1; i <= NH - 1; i++)
	{
		hr[i] = h0[i] - tau * (kI[i] * h0[i]);
		gr[i] = g0[i] + tau * ((Pgn * cn0[i] + Pga * ca0[i]) * (1. - g0[i]) - Qgh * h0[i] * MMg[i]);
		vr[i] = v0[i] + tau * (-omega * (cn0[i] + ca0[i]) * v0[i] - Mv * v0[i]);
		cnr[i] = cn0[i] + tau * (Vn * MMvst[i] * ca0[i] - Vd * MMv[i] * cn0[i] - kI[i] * cn0[i]);
		if (cn0[i] > 1.)
			cnr[i] = cnr[i] - tau * (mu * (cn0[i] - 1.));
		car[i] = ca0[i] + tau * (R * MMv[i] * (cn0[i] + ca0[i]) * (1. - (cn0[i] + ca0[i]) / cmax) - Vn * MMvst[i] * ca0[i] + Vd * MMv[i] * cn0[i] - kI[i] * ca0[i]);
	}

	hr[NH] = h0[NH] + tau * ((-2.) * h0[NH] * Is0[NH] / (NH * h + hNH / 2.));
	vr[NH] = v0[NH] + tau * (-omega * (cn0[NH] + ca0[NH]) * v0[NH] - Mv * v0[NH]);
	cnr[NH] = cn0[NH] + tau * (Vn * MMvst[NH] * ca0[NH] - Vd * MMv[NH] * cn0[NH] - 2. * cn0[NH] * Is0[NH] / ((1. * NH + 0.5) * h));
	if (cn0[NH] > 1.)
		cnr[NH] = cnr[NH] - tau * (mu * (cn0[NH] - 1.));
	car[NH] = ca0[NH] + tau * (R * MMv[NH] * (cn0[NH] + ca0[NH]) * (1. - (cn0[NH] + ca0[NH]) / cmax) - Vn * MMvst[NH] * ca0[NH] + Vd * MMv[NH] * cn0[NH] - 2. * ca0[NH] * Is0[NH] / ((1. * NH + 0.5) * h));
	gr[NH] = g0[NH] + tau * ((Pgn * cn0[NH] + Pga * ca0[NH]) * (1. - g0[NH]) - Qgh * h0[NH] * MMg[NH]);

	if (CTisOn == 1)
	{
		for (int i = 0; i <= NH; i++)
			kIf[i] = 2. * If0[i] / ((1. * i + 0.5) * h);

		ubl = ubl - tau * Cu * ubl;

		for (int i = 0; i <= NH; i++)
		{
			if (pc + sigma[i] >= 0.)
				ur[i] = u0[i] + tau * (((Ln * cn0[i] * sun + La * ca0[i] * sua) * (pc + sigma[i])) * ubl + (Pun * cn0[i] + Pua * ca0[i]) * (ubl - u0[i]) + Ll * h0[i] * sigma[i] * u0[i] - kIf[i] * u0[i]);
			else
				ur[i] = u0[i] + tau * (((Ln * cn0[i] * sun + La * ca0[i] * sua) * (pc + sigma[i])) * u0[i] + (Pun * cn0[i] + Pua * ca0[i]) * (ubl - u0[i]) + Ll * h0[i] * sigma[i] * u0[i]  - kIf[i] * u0[i]);
		}
	}
}

//********* TRIDIAGONAL MATRIX ALGORITHM

void TMA(double apr, double* bpr, double cpr, double* F, double* prog,
	double muleft, double hileft, double muright, double hiright)
{
	pp[1] = hileft;
	qq[1] = muleft;

	for (int i = 1; i <= NH - 1; i++)
	{
		pp[i + 1] = cpr / (bpr[i] - pp[i] * apr);
		qq[i + 1] = (apr * qq[i] + F[i]) / (bpr[i] - pp[i] * apr);
	}

	prog[NH] = (muright + hiright * qq[NH]) / (1. - pp[NH] * hiright);

	for (int i = NH - 1; i >= 0; i--)
		prog[i] = pp[i + 1] * prog[i + 1] + qq[i + 1];
}

void TMAc(double apr, double bpr, double cpr, double* F, double* prog,
	double muleft, double hileft, double muright, double hiright)
{
	pp[1] = hileft;
	qq[1] = muleft;

	for (int i = 1; i <= NH - 1; i++)
	{
		pp[i + 1] = cpr / (bpr - pp[i] * apr);
		qq[i + 1] = (apr * qq[i] + F[i]) / (bpr - pp[i] * apr);
	}

	prog[NH] = (muright + hiright * qq[NH]) / (1. - pp[NH] * hiright);

	for (int i = NH - 1; i >= 0; i--)
		prog[i] = pp[i + 1] * prog[i + 1] + qq[i + 1];
}

//***********************************************************************************************************************************************************
//********************************* CONVECTION OF CELLS *****************************************************************************************************
//***********************************************************************************************************************************************************

//********* PROCEDURE FOR COMPUTING CONVECTIVE SPEED

void Stress(double* nh)
{
	for (int i = 0; i <= NH; i++)
	{
		if (nh[i] < ss)
			sigma[i] = 0.;
		else if (nh[i] < 1)
			sigma[i] = (min(15000., k * (nh[i] - s0) * (nh[i] - ss) * (nh[i] - ss) / pow(1. - nh[i], 0.1)));
		else
			sigma[i] = (15000.); // to avoid technical difficulties, here just any large number should be
	}
}

double StressSingle(double nh)
{
	if (nh < ss)
		return (0.);
	else if (nh < 1.)
		return (min(15000., k * (nh - s0) * (nh - ss) * (nh - ss) / pow(1. - nh, 0.1)));
	else
		return (15000.);
}

void ConvSpeed(double* ccn, double* cca, double* hh, double* si, double* II, double* IIf, double* nh)
{
	II[0] = ((Ln * ccn[0] + La * cca[0]) * (pc + si[0]) + Ll * hh[0] * si[0]) * 0.25 * h * 0.25 * h * 0.5 * h;

	for (int i = 1; i <= NH - 1; i++)
		II[i] = II[i - 1] + ((Ln * ccn[i - 1] + La * cca[i - 1]) * (pc + si[i - 1]) + Ll * hh[i - 1] * si[i - 1]) * ((1. * i - 0.25) * h) * ((1. * i - 0.25) * h) * 0.5 * h + ((Ln * ccn[i] + La * cca[i]) * (pc + si[i]) + Ll * hh[i] * si[i]) * ((1. * i + 0.25) * h) * ((1. * i + 0.25) * h) * 0.5 * h;

	II[NH] = II[NH - 1] + ((Ln * ccn[NH - 1] + La * cca[NH - 1]) * (pc + si[NH - 1]) + Ll * hh[NH - 1] * si[NH - 1]) * ((1. * NH - 0.25) * h) * ((1. * NH - 0.25) * h) * 0.5 * h + ((Ln * ccn[NH] + La * cca[NH]) * (pc + si[NH]) + Ll * hh[NH] * si[NH]) * (NH * h + 0.5 * hNH) * (NH * h + 0.5 * hNH) * hNH;

	for (int i = 0; i <= NH - 1; i++)
		II[i] = II[i] / (((1. * i + 0.5) * h) * ((1. * i + 0.5) * h));

	II[NH] = II[NH] / ((NH * h + hNH) * (NH * h + hNH));

	if (CTisOn == 1)
	{
#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
		for (int i = 1; i <= NH - 2; i++)
		{
			IIf[i] = II[i] + (1. - 1. / (1 - nh[i])) * K * (sigma[i - 1] - sigma[i + 1]) / (2. * h);
			II[i] = II[i] + K * (sigma[i - 1] - sigma[i + 1]) / (2. * h);
		}

		IIf[NH - 1] = II[NH - 1] + (1. - 1. / (1 - nh[NH - 1])) * K * (sigma[NH - 2] - (sigma[NH - 1] + h * (sigma[NH] - sigma[NH - 1]) / (hNH + h / 2.))) / (2. * h);
		II[NH - 1] = II[NH - 1] + K * (sigma[NH - 2] - (sigma[NH - 1] + h * (sigma[NH] - sigma[NH - 1]) / (hNH + h / 2.))) / (2. * h);

		IIf[NH] = II[NH];
	}
	else
	{
		II[0] = II[0] + K * (sigma[0] - sigma[1]) / (2. * h);

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
		for (int i = 1; i <= NH - 2; i++)
			II[i] = II[i] + K * (sigma[i - 1] - sigma[i + 1]) / (2. * h);

		II[NH - 1] = II[NH - 1] + K * (sigma[NH - 2] - (sigma[NH - 1] + h * (sigma[NH] - sigma[NH - 1]) / (hNH + h / 2.))) / (2. * h);
	}
}

//********* SIGN FUNCTION

double Sgn(double x)
{
	if (x > 0.)
		return (1.);
	else if (x == 0.)
		return (0.);
	else
		return (-1.);
}

//********* PROCEDURE FOR CELLS' CONVECTION (more precisely, its "linear" part)

void Convection()
{
	for (int i = 0; i <= NN - 1; i++)
		nh0[i] = npr[i] + nqr[i] + mr[i];

	nh0[NN] = (npr[NN] + nqr[NN] + mr[NN]) * (hNN / h) + hr[NN] * ((h - hNN) / h);

	for (int i = NN + 1; i <= NH; i++)
		nh0[i] = hr[i];

	Stress(nh0);
	ConvSpeed(cn0, ca0, h0, sigma, Is0_not_corr, If0_not_corr, nh0);

	// adjusting speed
	if (Is0_not_corr[0] >= 0.)
		Is0[0] = (1. - (Is0_not_corr[0] * tau / (2. * h))) * Is0_not_corr[0] + (Is0_not_corr[0] * tau / (2. * h)) * Is0_not_corr[1];
	else
		Is0[0] = (1. + (Is0_not_corr[0] * tau / (2. * h))) * Is0_not_corr[0] + (Is0_not_corr[0] * tau / (2. * h)) * Is0_not_corr[0];

	for (int i = 1; i <= NH - 1; i++)
	{
		if (Is0_not_corr[i] >= 0.)
			Is0[i] = (1. - (Is0_not_corr[i] * tau / (2. * h))) * Is0_not_corr[i] + (Is0_not_corr[i] * tau / (2. * h)) * Is0_not_corr[i + 1];
		else
			Is0[i] = (1. + (Is0_not_corr[i] * tau / (2. * h))) * Is0_not_corr[i] - (Is0_not_corr[i] * tau / (2. * h)) * Is0_not_corr[i - 1];
	}

	if (Is0_not_corr[NH] < 0.)
		Is0[NH] = (1. + (Is0_not_corr[NH] * tau / (h + hNH))) * Is0_not_corr[NH] - (Is0_not_corr[NH] * tau / (h + hNH)) * Is0_not_corr[NH - 1];
	else
		Is0[NH] = Is0_not_corr[NH];

	// defining auxiliary variables
	for (int i = 0; i <= NH - 2; i++)
		Q_plus[i] = (0.5 - Is0[i] * tau / h) / (1. + (Is0[i + 1] - Is0[i]) * tau / h);
	for (int i = 1; i <= NH - 1; i++)
		Q_minus[i] = (0.5 + Is0[i] * tau / h) / (1. - (Is0[i - 1] - Is0[i]) * tau / h);

	Q_plus[NH - 1] = (0.5 * h - Is0[NH - 1] * tau) / (0.5 * (h + hNH) + (Is0[NH] - Is0[NH - 1]) * tau);

	Q_minus[NH] = (0.5 * hNH + Is0[NH] * tau) / (0.5 * (h + hNH) + (Is0[NH] - Is0[NH - 1]) * tau);

	// calculating the result of the transport stage
	nh1[0] = 0.5 * Q_plus[0] * Q_plus[0] * (nh0[1] - nh0[0]) + Q_plus[0] * nh0[0] + 0.5 * nh0[0];

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NH - 2; i++)
		nh1[i] = 0.5 * Q_minus[i] * Q_minus[i] * (nh0[i - 1] - nh0[i]) + 0.5 * Q_plus[i] * Q_plus[i] * (nh0[i + 1] - nh0[i]) + (Q_plus[i] + Q_minus[i]) * nh0[i];

	nh1[NH - 1] = 0.5 * Q_minus[NH - 1] * Q_minus[NH - 1] * (nh0[NH - 2] - nh0[NH - 1]) + 0.5 * Q_plus[NH - 1] * Q_plus[NH - 1] * (nh0[NH] - nh0[NH - 1]) * (h + hNH) / (2. * h) + (Q_plus[NH - 1] * (h + hNH) / (2. * h) + Q_minus[NH - 1]) * nh0[NH - 1];

	nh1[NH] = hstat;

	// quiescent cells separately 
	nq1[0] = 0.5 * Q_plus[0] * Q_plus[0] * (nqr[1] - nqr[0]) + Q_plus[0] * nqr[0] + 0.5 * nqr[0];

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NN - 1; i++)
		nq1[i] = 0.5 * Q_minus[i] * Q_minus[i] * (nqr[i - 1] - nqr[i]) + 0.5 * Q_plus[i] * Q_plus[i] * (nqr[i + 1] - nqr[i]) + (Q_plus[i] + Q_minus[i]) * nqr[i];

	nq1[NN] = nq1[NN - 1];

	// damaged cells separately 
	m1[0] = 0.5 * Q_plus[0] * Q_plus[0] * (mr[1] - mr[0]) + Q_plus[0] * mr[0] + 0.5 * mr[0];

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NN - 1; i++)
		m1[i] = 0.5 * Q_minus[i] * Q_minus[i] * (mr[i - 1] - mr[i]) + 0.5 * Q_plus[i] * Q_plus[i] * (mr[i + 1] - mr[i]) + (Q_plus[i] + Q_minus[i]) * mr[i];

	m1[NN] = m1[NN - 1];

	// capillaries
	cn1[0] = 0.5 * Q_plus[0] * Q_plus[0] * (cnr[1] - cnr[0]) + Q_plus[0] * cnr[0] + 0.5 * cnr[0];

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NH - 2; i++)
		cn1[i] = 0.5 * Q_minus[i] * Q_minus[i] * (cnr[i - 1] - cnr[i]) + 0.5 * Q_plus[i] * Q_plus[i] * (cnr[i + 1] - cnr[i]) + (Q_plus[i] + Q_minus[i]) * cnr[i];

	cn1[NH - 1] = 0.5 * Q_minus[NH - 1] * Q_minus[NH - 1] * (cnr[NH - 2] - cnr[NH - 1]) + 0.5 * Q_plus[NH - 1] * Q_plus[NH - 1] * (cnr[NH] - cnr[NH - 1]) * (h + hNH) / (2. * h) + (Q_plus[NH - 1] * (h + hNH) / (2. * h) + Q_minus[NH - 1]) * cnr[NH - 1];
	cn1[NH] = 1.;


	cac[0] = 0.5 * Q_plus[0] * Q_plus[0] * (car[1] - car[0]) + Q_plus[0] * car[0] + 0.5 * car[0];

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NH - 2; i++)
		cac[i] = 0.5 * Q_minus[i] * Q_minus[i] * (car[i - 1] - car[i]) + 0.5 * Q_plus[i] * Q_plus[i] * (car[i + 1] - car[i]) + (Q_plus[i] + Q_minus[i]) * car[i];

	cac[NH - 1] = 0.5 * Q_minus[NH - 1] * Q_minus[NH - 1] * (car[NH - 2] - car[NH - 1]) + 0.5 * Q_plus[NH - 1] * Q_plus[NH - 1] * (car[NH] - car[NH - 1]) * (h + hNH) / (2. * h) + (Q_plus[NH - 1] * (h + hNH) / (2. * h) + Q_minus[NH - 1]) * car[NH - 1];
	cac[NH] = 0.;

	//********* ANTIDIFFUSION STAGE
	for (int i = 0; i <= NH - 1; i++)
	{
		delta_nh[i] = nh1[i + 1] - nh1[i];
		delta_cn[i] = cn1[i + 1] - cn1[i];
		delta_ca[i] = cac[i + 1] - cac[i];
	}

	delta_nh[NH] = 0. - nh1[NH];
	delta_cn[NH] = 0. - cn1[NH];
	delta_ca[NH] = 0. - cac[NH];

	for (int i = 0; i <= NH + 1; i++)
		one_eighth[i] = 0.125 + 0.5 * (Is0[i] * tau / h) * (Is0[i] * tau / h);

	f_c_nh[0] = Sgn(delta_nh[0]) * max(0., min((-1.) * delta_nh[0] * Sgn(delta_nh[0]), min(one_eighth[0] * abs(delta_nh[0]), delta_nh[1] * Sgn(delta_nh[0]))));
	f_c_cn[0] = Sgn(delta_cn[0]) * max(0., min((-1.) * delta_cn[0] * Sgn(delta_cn[0]), min(one_eighth[0] * abs(delta_cn[0]), delta_cn[1] * Sgn(delta_cn[0]))));
	f_c_ca[0] = Sgn(delta_ca[0]) * max(0., min((-1.) * delta_ca[0] * Sgn(delta_ca[0]), min(one_eighth[0] * abs(delta_ca[0]), delta_ca[1] * Sgn(delta_ca[0]))));

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NH - 1; i++)
	{
		f_c_nh[i] = Sgn(delta_nh[i]) * max(0., min(delta_nh[i - 1] * Sgn(delta_nh[i]), min(one_eighth[i] * abs(delta_nh[i]), delta_nh[i + 1] * Sgn(delta_nh[i]))));
		f_c_cn[i] = Sgn(delta_cn[i]) * max(0., min(delta_cn[i - 1] * Sgn(delta_cn[i]), min(one_eighth[i] * abs(delta_cn[i]), delta_cn[i + 1] * Sgn(delta_cn[i]))));
		f_c_ca[i] = Sgn(delta_ca[i]) * max(0., min(delta_ca[i - 1] * Sgn(delta_ca[i]), min(one_eighth[i] * abs(delta_ca[i]), delta_ca[i + 1] * Sgn(delta_ca[i]))));
	}

	nh1[0] = nh1[0] - f_c_nh[0] - Sgn(delta_nh[0]) * max(0., min(delta_nh[1] * Sgn(delta_nh[0]), min((-1.) * one_eighth[1] * abs(delta_nh[0]), (-1.) * delta_nh[0] * Sgn(delta_nh[0]))));
	cn1[0] = cn1[0] - f_c_cn[0] - Sgn(delta_cn[0]) * max(0., min(delta_cn[1] * Sgn(delta_cn[0]), min((-1.) * one_eighth[1] * abs(delta_cn[0]), (-1.) * delta_cn[0] * Sgn(delta_cn[0]))));
	cac[0] = cac[0] - f_c_ca[0] - Sgn(delta_ca[0]) * max(0., min(delta_ca[1] * Sgn(delta_ca[0]), min((-1.) * one_eighth[1] * abs(delta_ca[0]), (-1.) * delta_ca[0] * Sgn(delta_ca[0]))));

	for (int i = 1; i <= NH - 1; i++)
	{
		nh1[i] = nh1[i] - f_c_nh[i] + f_c_nh[i - 1];
		cn1[i] = cn1[i] - f_c_cn[i] + f_c_cn[i - 1];
		cac[i] = cac[i] - f_c_ca[i] + f_c_ca[i - 1];
	}

	nh1[NH] = nh1[NH] + f_c_nh[NH - 1];
	cn1[NH] = cn1[NH] + f_c_cn[NH - 1];
	cac[NH] = cac[NH] + f_c_ca[NH - 1];

	// quiescent cells separately 
	for (int i = 0; i <= NN - 1; i++)
		delta_nq[i] = nq1[i + 1] - nq1[i];

	delta_nq[NN] = 0. - nq1[NN];

	f_c_nq[0] = Sgn(delta_nq[0]) * max(0., min((-1.) * delta_nq[0] * Sgn(delta_nq[0]), min(one_eighth[0] * abs(delta_nq[0]), delta_nq[1] * Sgn(delta_nq[0]))));

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NN - 1; i++)
		f_c_nq[i] = Sgn(delta_nq[i]) * max(0., min(delta_nq[i - 1] * Sgn(delta_nq[i]), min(one_eighth[i] * abs(delta_nq[i]), delta_nq[i + 1] * Sgn(delta_nq[i]))));

	nq1[0] = nq1[0] - f_c_nq[0] - Sgn(delta_nq[0]) * max(0., min(delta_nq[1] * Sgn(delta_nq[0]), min((-1.) * one_eighth[1] * abs(delta_nq[0]), (-1.) * delta_nq[0] * Sgn(delta_nq[0]))));

	for (int i = 1; i <= NN - 1; i++)
		nq1[i] = nq1[i] - f_c_nq[i] + f_c_nq[i - 1];

	nq1[NN] = nq1[NN] + f_c_nq[NN - 1];

	// damaged cells separately 
	for (int i = 0; i <= NN - 1; i++)
		delta_m[i] = m1[i + 1] - m1[i];

	delta_m[NN] = 0. - m1[NN];

	f_c_m[0] = Sgn(delta_m[0]) * max(0., min((-1.) * delta_m[0] * Sgn(delta_m[0]), min(one_eighth[0] * abs(delta_m[0]), delta_m[1] * Sgn(delta_m[0]))));

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NN - 1; i++)
		f_c_m[i] = Sgn(delta_m[i]) * max(0., min(delta_m[i - 1] * Sgn(delta_m[i]), min(one_eighth[i] * abs(delta_m[i]), delta_m[i + 1] * Sgn(delta_m[i]))));

	m1[0] = m1[0] - f_c_m[0] - Sgn(delta_m[0]) * max(0., min(delta_m[1] * Sgn(delta_m[0]), min((-1.) * one_eighth[1] * abs(delta_m[0]), (-1.) * delta_m[0] * Sgn(delta_m[0]))));

	for (int i = 1; i <= NN - 1; i++)
		m1[i] = m1[i] - f_c_m[i] + f_c_m[i - 1];

	m1[NN] = m1[NN] + f_c_m[NN - 1];

	//********* COUNTING "NORMAL CELLS" BEFORE CONVECTIVE STAGE (since this part of equation is solved in linear regime and this variable should be conserved)
	// it is made now since further hNN will probably change and NN may change

	if (R0T != size)
	{
		HTCrLIN = h0[NN] * (h - hNN);

		for (int i = NN + 1; i <= NH - 1; i++)
			HTCrLIN = HTCrLIN + h0[i] * h;

		HTCrLIN = HTCrLIN + h0[NH] * hNH;


		HTCr = h0[NN] * (h * NN + hNN / 2.) * (h * NN + hNN / 2.) * (h - hNN);

		for (int i = NN + 1; i <= NH - 1; i++)
			HTCr = HTCr + h0[i] * ((1. * i + 0.5) * h) * ((1. * i + 0.5) * h) * h;

		HTCr = HTCr + h0[NH] * (NH * h + hNH / 2.) * (NH * h + hNH / 2.) * hNH;
	}

	//********* COUNTING "PROLIFERATING TUMOR CELLS" BEFORE CONVECTIVE STAGE (since this part of equation is solved in linear regime and this variable should be conserved)

	NTCpr = 0.;
	for (int i = 0; i <= NN - 1; i++)
		NTCpr = NTCpr + npr[i] * h;

	NTCpr = NTCpr + npr[NN] * hNN;

	//********* COUNTING "TUMOR CELLS" BEFORE CONVECTIVE STAGE (since this part of equation is solved in linear regime and this variable should be conserved)

	NTCr = 0.;
	for (int i = 0; i <= NN - 1; i++)
		NTCr = NTCr + nh0[i] * h;

	NTCr = NTCr + nh0[NN] * hNN;

	//********* COUNTING "TUMOR CELLS" AFTER CONVECTIVE STAGE UP TO NN-1
	NTCcm1 = 0.;
	for (int i = 0; i <= NN - 1; i++)
		NTCcm1 = NTCcm1 + nh1[i] * h;

	hNNnew = (NTCr - NTCcm1) / nh1[NN];

	//********* COUNTING "PROLIFERATING TUMOR CELLS" AFTER CONVECTIVE STAGE UP TO NN-1
	NTCpm1 = 0.;
	for (int i = 0; i <= NN - 1; i++)
		NTCpm1 = NTCpm1 + np1[i] * h;

	// REDISTRIBUTING n AND h

	for (int i = 0; i <= NN - 1; i++)
	{
		np1[i] = nh1[i] - nq1[i] - m1[i];
		h1[i] = 0.;
	}
	np1[NN] = nh1[NN] - nq1[NN] - m1[NN];

	if (hNNnew < -1. * h)
	{
		cout << "\nError! hNNnew <-h";
		system("Pause");
		exit(0);
	}
	else if (hNNnew < 0.)
	{
		np1[NN] = 0.;
		nq1[NN] = 0.;
		m1[NN] = 0.;
		h1[NN] = nh1[NN];
		NN--;
		h1[NN] = nh1[NN];
		hNN = h + hNNnew;
	}
	else if (hNNnew < h)
	{
		h1[NN] = nh1[NN];
		hNN = hNNnew;
	}
	else if (hNNnew < 2. * h)
	{
		h1[NN] = 0.;
		NN++;
		np1[NN] = np1[NN - 1];
		nq1[NN] = nq1[NN - 1];
		m1[NN] = m1[NN - 1];
		h1[NN] = nh1[NN];
		hNN = (NTCr - NTCcm1 - nh1[NN - 1] * h) / nh1[NN];
	}
	else
	{
		cout << "\nError! hNNnew >= 2.*h";
		system("Pause");
		exit(0);
	}

	if (R0T == size)
	{
		cout << "\n R0T == size \n";
		system("PAUSE");
		exit(0);
	}
	else
	{
		//********* COUNTING "NORMAL CELLS" AFTER CONVECTIVE STAGE UP TO NN-1

		HTCcm1LIN = nh1[NN] * (h - hNN);

		for (int i = NN + 1; i <= NH - 1; i++)
			HTCcm1LIN = HTCcm1LIN + nh1[i] * h;

		hNHnew = (HTCrLIN - HTCcm1LIN) / nh1[NH];

		HTCcm1 = nh1[NN] * (h * NN + hNN / 2.) * (h * NN + hNN / 2.) * (h - hNN);

		for (int i = NN + 1; i <= NH - 1; i++)
			HTCcm1 = HTCcm1 + nh1[i] * ((1. * i + 0.5) * h) * ((1. * i + 0.5) * h) * h;

		hNHnew = (HTCr - HTCcm1) / (nh1[NH] * (NH * h + hNHnew / 2.) * (NH * h + hNHnew / 2.));
		hNHnew = (HTCr - HTCcm1) / (nh1[NH] * (NH * h + hNHnew / 2.) * (NH * h + hNHnew / 2.));
		hNHnew = (HTCr - HTCcm1) / (nh1[NH] * (NH * h + hNHnew / 2.) * (NH * h + hNHnew / 2.));

		if (hNHnew < -1. * h)
		{
			cout << "\nError! hNHnew <-h";
			system("Pause");
			exit(0);
		}
		else if (hNHnew < 0.)
		{
			nh1[NH] = 0.;
			u1[NH] = 0.;
			NH--;
			hNH = h + hNHnew;
		}
		else if (hNHnew < h)
			hNH = hNHnew;
		else if (hNHnew < 2. * h)
		{
			NH++;
			nh1[NH] = nh1[NH - 1];
			g0[NH] = g0[NH - 1];
			u1[NH] = u1[NH - 1];
			hNH = hNHnew - h;
		}
		else if (hNHnew < 3. * h)
		{
			NH++;
			nh1[NH] = nh1[NH - 1];
			g0[NH] = g0[NH - 1];
			u1[NH] = u1[NH - 1];
			NH++;
			nh1[NH] = nh1[NH - 1];
			g0[NH] = g0[NH - 1];
			u1[NH] = u1[NH - 1];
			hNH = hNHnew - 2. * h;
		}
		else
		{
			cout << "\nError! hNHnew >= 3.*h";
			system("Pause");
			exit(0);
		}
	}

	// REDISTRIBUTING n AND h
	for (int i = NN + 1; i <= NH; i++)
	{
		np1[i] = 0.;
		nq1[i] = 0.;
		m1[i] = 0.;
		h1[i] = nh1[i];
	}
}

void ConvectionCA()
{
	// adjusting speed
	if (If0_not_corr[0] >= 0.)
		If0[0] = (1. - (If0_not_corr[0] * tau / (2. * h))) * If0_not_corr[0] + (If0_not_corr[0] * tau / (2. * h)) * If0_not_corr[1];
	else
		If0[0] = (1. + (If0_not_corr[0] * tau / (2. * h))) * If0_not_corr[0] + (If0_not_corr[0] * tau / (2. * h)) * If0_not_corr[0];

	for (int i = 1; i <= NH - 1; i++)
	{
		if (If0_not_corr[i] >= 0.)
			If0[i] = (1. - (If0_not_corr[i] * tau / (2. * h))) * If0_not_corr[i] + (If0_not_corr[i] * tau / (2. * h)) * If0_not_corr[i + 1];
		else
			If0[i] = (1. + (If0_not_corr[i] * tau / (2. * h))) * If0_not_corr[i] - (If0_not_corr[i] * tau / (2. * h)) * If0_not_corr[i - 1];
	}

	if (If0_not_corr[NH] < 0.)
		If0[NH] = (1. + (If0_not_corr[NH] * tau / (h + hNH))) * If0_not_corr[NH] - (If0_not_corr[NH] * tau / (h + hNH)) * If0_not_corr[NH - 1];
	else
		If0[NH] = If0_not_corr[NH];

	// defining auxiliary variables
	for (int i = 0; i <= NH - 2; i++)
		Q_plus[i] = (0.5 - If0[i] * tau / h) / (1. + (If0[i + 1] - If0[i]) * tau / h);
	for (int i = 1; i <= NH - 1; i++)
		Q_minus[i] = (0.5 + If0[i] * tau / h) / (1. - (If0[i - 1] - If0[i]) * tau / h);

	Q_plus[NH - 1] = (0.5 * h - If0[NH - 1] * tau) / (0.5 * (h + hNH) + (If0[NH] - If0[NH - 1]) * tau);

	Q_minus[NH] = (0.5 * hNH + If0[NH] * tau) / (0.5 * (h + hNH) + (If0[NH] - If0[NH - 1]) * tau);

	// calculating the result of the transport stage
	uc[0] = 0.5 * Q_plus[0] * Q_plus[0] * (ur[1] - ur[0]) + Q_plus[0] * ur[0] + 0.5 * ur[0];

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NH - 2; i++)
		uc[i] = 0.5 * Q_minus[i] * Q_minus[i] * (ur[i - 1] - ur[i]) + 0.5 * Q_plus[i] * Q_plus[i] * (ur[i + 1] - ur[i]) + (Q_plus[i] + Q_minus[i]) * ur[i];

	uc[NH - 1] = 0.5 * Q_minus[NH - 1] * Q_minus[NH - 1] * (ur[NH - 2] - ur[NH - 1]) + 0.5 * Q_plus[NH - 1] * Q_plus[NH - 1] * (ur[NH] - ur[NH - 1]) * (h + hNH) / (2. * h) + (Q_plus[NH - 1] * (h + hNH) / (2. * h) + Q_minus[NH - 1]) * ur[NH - 1];
	uc[NH] = uc[NH - 1];

	//********* ANTIDIFFUSION STAGE
	for (int i = 0; i <= NH - 1; i++)
		delta_uf[i] = uc[i + 1] - uc[i];

	delta_uf[NH] = 0. - uc[NH];

	for (int i = 0; i <= NH + 1; i++)
		one_eighth[i] = 0.125 + 0.5 * (If0[i] * tau / h) * (If0[i] * tau / h);

	f_c_uf[0] = Sgn(delta_uf[0]) * max(0., min((-1.) * delta_uf[0] * Sgn(delta_uf[0]), min(one_eighth[0] * abs(delta_uf[0]), delta_uf[1] * Sgn(delta_uf[0]))));

#pragma loop(hint_parallel(6))  
#pragma loop(ivdep) 
	for (int i = 1; i <= NH - 1; i++)
		f_c_uf[i] = Sgn(delta_uf[i]) * max(0., min(delta_uf[i - 1] * Sgn(delta_uf[i]), min(one_eighth[i] * abs(delta_uf[i]), delta_uf[i + 1] * Sgn(delta_uf[i]))));

	uc[0] = uc[0] - f_c_uf[0] - Sgn(delta_uf[0]) * max(0., min(delta_uf[1] * Sgn(delta_uf[0]), min((-1.) * one_eighth[1] * abs(delta_uf[0]), (-1.) * delta_uf[0] * Sgn(delta_uf[0]))));

	for (int i = 1; i <= NH - 1; i++)
		uc[i] = uc[i] - f_c_uf[i] + f_c_uf[i - 1];

	uc[NH] = uc[NH] + f_c_uf[NH - 1];
}

//***********************************************************************************************************************************************************
//********************************* DIFFUSION ***************************************************************************************************************
//***********************************************************************************************************************************************************

void Diffusion_glucose()
{
	for (int i = 0; i <= NH; i++)
	{
		G0[i] = g0[i] * (1. * i + 0.5) * h; // for easy solution in spherical symmetrical case
		GR[i] = gr[i] * (1. * i + 0.5) * h;
	}

	for (int i = 1; i <= NH - 1; i++)
		Fg[i] = (G0[i + 1] - 2. * G0[i] + G0[i - 1]) + kg * GR[i];

	muleftg = 0.;
	hileftg = 0.5 / 1.5;

	murightg = 0.;
	hirightg = (1. * NH + 0.5) / (1. * NH - 0.5);

	TMAc(aprg, bprg, cprg, Fg, progg, muleftg, hileftg, murightg, hirightg);

	for (int i = 0; i <= NH; i++)
		G1[i] = progg[i];

	for (int i = 0; i <= NH; i++)
		g1[i] = G1[i] / ((1. * i + 0.5) * h);
}

void Diffusion_VEGF()
{
	for (int i = 0; i <= NH; i++)
	{
		V0[i] = v0[i] * (1. * i + 0.5) * h; // for easy solution in spherical symmetrical case
		VR[i] = vr[i] * (1. * i + 0.5) * h;
	}

	for (int i = 1; i <= NH - 1; i++)
		Fv[i] = (V0[i + 1] - 2. * V0[i] + V0[i - 1]) + kv * VR[i];

	muleftv = 0.;
	hileftv = 0.5 / 1.5;

	murightv = 0.;
	hirightv = (1. * NH + 0.5) / (1. * NH - 0.5);

	TMAc(aprv, bprv, cprv, Fv, progv, muleftv, hileftv, murightv, hirightv);

	for (int i = 0; i <= NH; i++)
		V1[i] = progv[i];

	for (int i = 0; i <= NH; i++)
		v1[i] = V1[i] / ((1. * i + 0.5) * h);
}

void Diffusion_cap()
{
	for (int i = 0; i <= NH; i++)
	{
		CA0[i] = ca0[i] * (1. * i + 0.5) * h; // for easy solution in spherical symmetrical case
		CAR[i] = cac[i] * (1. * i + 0.5) * h;
	}

	for (int i = 1; i <= NH - 1; i++)
		Fc[i] = (CA0[i + 1] - 2. * CA0[i] + CA0[i - 1]) + kc * CAR[i];

	muleftc = 0.;
	hileftc = 0.5 / 1.5;

	murightc = 0.;
	hirightc = (1. * NH + 0.5) / (1. * NH - 0.5);

	TMAc(aprc, bprc, cprc, Fc, progc, muleftc, hileftc, murightc, hirightc);

	for (int i = 0; i <= NH; i++)
		CA1[i] = progc[i];

	for (int i = 0; i <= NH; i++)
		ca1[i] = CA1[i] / ((1. * i + 0.5) * h);
}

void Diffusion_CA()
{
	for (int i = 0; i <= NH; i++)
	{
		U0[i] = u0[i] * (1. * i + 0.5) * h; // for easy solution in spherical symmetrical case
		UR[i] = uc[i] * (1. * i + 0.5) * h;
	}

	for (int i = 1; i <= NH - 1; i++)
		Fu[i] = (U0[i + 1] - 2. * U0[i] + U0[i - 1]) + ku * UR[i];

	muleftu = 0.;
	hileftu = 0.5 / 1.5;

	murightu = 0.;
	hirightu = (1. * NH + 0.5) / (1. * NH - 0.5);

	TMAc(apru, bpru, cpru, Fu, progu, muleftu, hileftu, murightu, hirightu);

	for (int i = 0; i <= NH; i++)
		U1[i] = progu[i];

	for (int i = 0; i <= NH; i++)
		u1[i] = U1[i] / ((1. * i + 0.5) * h);
}
//***********************************************************************************************************************************************************
//********************************* TIME STEP ADJUSTMENT ****************************************************************************************************
//***********************************************************************************************************************************************************

void TauAdjust()
{
	ta = ta + tau_adjust_frequency;

	// adjusting tau
	Is0Max = 0.;
	Is0NonMonIndex = 0;

	for (int i = NN + 1; i <= NH; i++)
		Is0Max = max(Is0Max, h0[i]);

	for (int i = NN + 1; i <= NH - 1; i++)
		if ((h0[i] - h0[i - 1]) * (h0[i + 1] - h0[i]) < -1.e-10)
			Is0NonMonIndex = Is0NonMonIndex - (h0[i] - h0[i - 1]) * (h0[i + 1] - h0[i]);

	tauAdjust = h / (2. * Is0Max) / tauAdjustParameter;

	if (Is0NonMonIndex > Is0NonMonIndexMin)
		tau = max(tauMin, tau / 2.);
	else if (tau < tauAdjust)
		tau = min(tauMax, 1.02 * tau);
	else
		tau = max(tauMin, tauAdjust);

	tau = 1.e-8 * (int)(1.e+8 * tau);

	kg = 2 * h * h / (tau * Dg);
	bprg = 2. + kg;
}

//***********************************************************************************************************************************************************
//********************************* MAIN CYCLE ENDING *******************************************************************************************************
//***********************************************************************************************************************************************************

void Cycle_Ending()
{
	for (int i = 0; i <= NH; i++)
	{
		np0[i] = np1[i];
		nq0[i] = nq1[i];
		h0[i] = h1[i];
		m0[i] = m1[i];
		v0[i] = v1[i];
		cn0[i] = cn1[i];
		ca0[i] = ca1[i];
		g0[i] = g1[i];
		u0[i] = u1[i];
	}
}

//***********************************************************************************************************************************************************
//********************************* PROGRAM ENDING **********************************************************************************************************
//***********************************************************************************************************************************************************

void End()
{
	if (clock() / CLOCKS_PER_SEC >= 3600)
		printf("\n%s%ld%s%ld%s%ld%s", "The end. Total time: ", (clock() / CLOCKS_PER_SEC) / 3600, " h ", (clock() / CLOCKS_PER_SEC) / 60 - 60 * ((clock() / CLOCKS_PER_SEC) / 3600), " min ", (clock() / CLOCKS_PER_SEC) - 3600 * ((clock() / CLOCKS_PER_SEC) / 3600) - 60 * ((clock() / CLOCKS_PER_SEC) / 60 - 60 * ((clock() / CLOCKS_PER_SEC) / 3600)), " s.");
	else if (clock() / CLOCKS_PER_SEC >= 60)
		printf("\n%s%ld%s%ld%s", "The end. Total time: ", (clock() / CLOCKS_PER_SEC) / 60 - 60 * ((clock() / CLOCKS_PER_SEC) / 3600), " min ", (clock() / CLOCKS_PER_SEC) - 3600 * ((clock() / CLOCKS_PER_SEC) / 3600) - 60 * ((clock() / CLOCKS_PER_SEC) / 60 - 60 * ((clock() / CLOCKS_PER_SEC) / 3600)), " s.");
	else
		printf("\n%s%ld%s", "The end. Total time: ", (clock() / CLOCKS_PER_SEC) - 3600 * ((clock() / CLOCKS_PER_SEC) / 3600) - 60 * ((clock() / CLOCKS_PER_SEC) / 60 - 60 * ((clock() / CLOCKS_PER_SEC) / 3600)), " s.");
	//_getch();
	exit(0);
}

//********************************************************
//********* THE END
//********************************************************