//LLNL Maturity Prediction System
//ReWritten with C++ by ZHS(水南)
/* references:
(1)Vitrinite reflectance (%) using the LLNL kinetic model
of (Sweeney and Burnham, 1990)
(2)Basin 2
*/

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
using namespace std;

#define R 8.31432
#define A 1.0E13

int main()
{
double fm[20]={.03, .03, .04, .04, .05,
.05, .06, .04, .04, .07,
.06, .06, .06, .05, .05,
.04, .03, .02, .02, .01};
double Ea[20]={34.0, 36.0, 38.0, 40.0, 42.0,
44.0, 46.0, 48.0, 50.0, 52.0,
54.0, 56.0, 58.0, 60.0, 62.0,
64.0, 66.0, 68.0, 70.0, 72.0};
double tc[200]={0},tk[200],t[200]={0},ts[200],dtk=2,sumF[200]={0};
double dt,h[50],dc[200][20],F[200][20]={0},k[20];
ifstream getdata;
ofstream savefile;
int i,j,num,count=0;
getdata.open("./data.txt");
if (!getdata)
{
cout<<"It seems that you do not have the data file named data.txt"<<endl
<<"Please check it and then run the program again."<<endl;
cout<<"※※※※※※※※※※※※※※※※※※※"<<endl;
cout<<"Press any key to continue..."<<endl;
    getchar();
exit(1);
}
i=0;
while(!getdata.eof())
{
getdata>>t[i];
getdata>>tc[i];
i++;
}
num=i;
if(num>=200)
{
  cout<<"Error! You have too much data"<<endl;
  exit(1);
  }
  tk[0]=tc[0];
  ts[0]=t[0]; //ts表示反应时间，此处单位为Ma
/*
对时间和温度进行插值，使时间间隔为2Ma，时间范围为输入的时间范围
  */
  i=0;
  j=1;
  while(ts[i]<=t[num-1])
{
  if(ts[i]+dtk<=t[j])
   {
     ts[i+1]=ts[i]+dtk; //dtk 为时间间隔，默认为2K/Ma
       tk[i+1]=(ts[i+1]-t[j-1])*(tc[j]-tc[j-1])/(t[j]-t[j-1])+tc[j-1];
       i++;
     }
     else j++;
  }
  
for(i=0;i!=20;i++)
{
     Ea[i]=Ea[i]*4186.8;//转换活化能的单位Kcal->J
     F[0][i]=fm[i];
}
dtk=dtk*3.16E13; //时间间隔转换为S
for(i=1;ts[i]<=t[num-1];i++)
{
for (j=0;j!=20;j++)
k[j]=A*exp(-Ea[j]/R/(tk[i]+273));
for (j=0;j!=20;j++)
F[i][j]=F[i-1][j]*exp(-k[j]*dtk);
    count++;  
}
  savefile.open("./LLNL_VRo.txt",ios::app);
  for(i=0;i<=count;i++)  
  {
    for(j=0;j!=20;j++)
       sumF[i]=sumF[i]+F[i][j];
   sumF[i]=0.85-sumF[i];  
  }
savefile<<"Time(Ma)"<<"\t"<<"Temperature(C)"<<"\t"<<"VRo(%)"<<"\n";
  for(i=0;i<=count;i++)
    savefile<<ts[i]<<"\t"<<"\t"<<tk[i]<<"\t"<<"\t"<<exp(-1.6+3.7*sumF[i])<<"\n";
  cout<<"The result has been successfully written to ./LLNL_VRo.txt"<<endl
   <<"Press any key to continue"<<endl;
  getchar();

}