//源岩降解动力学过程在C++的初步实现 
 
//--------------------------------------------------------- 
 
//                !!CORPYRIGHT DECLARATION!! 
 
//This program was wirrten by Zeng H.S using std C++ language. 
 
//You can use the program under GNU licence. 
 
//                                2006/08/11(YY/MM/DD) 
 
//------------------------------------------------------------ 
 
#include <iostream> 
 
#include <math.h> 
 
#include<conio.h> 
 
#include <stdlib.h> 
 
#define R 8.31441 
 
#define T2 2 
 
#define T_initial 25 
 
 
 
using namespace std; 
 
 
 
int main() 
 
{ 
 
//--------------------------------------------------------- 
 
  //begin varialbe declaration 
 
  long double delta_t,time1,time2,A,Ea;   
 
                            //delta_t是累积时间，time是时间 
 
                        //A是指前因子，Ea是活化能 
 
  long double T[200],TOC_initial,C_initial,Tx,alpha_x=2,x[200],k; 
 
  //end varialble declaration 
 
  //--------------------------------------------------------- 
 
  //begin data initiation  
 
  T[0]=T_initial;                          //initial temperature 
 
  Tx=alpha_x+15*log10(alpha_x/2);         //heating rate 
 
  Ea=259.1*1E3;                           //A type kerogen  Ea 
 
  A=41.23E17;                              //frequency 
 
  //---------------------------------------------------------- 
 
  //user input 
 
L1: cout<<"Please input Heating Rate(C/Ma):"<<endl 
 
      <<"Heating Rate(C/Ma)="; 
 
  cin>>alpha_x; 
 
  cout<<endl; 
 
  cout<<"Please input Ea(KJ/Mol):"<<endl 
 
      <<"Ea(KJ/Mol)="; 
 
  cin>>Ea; 
 
  Ea=Ea*1000; 
 
  cout<<endl<<"Please input A(frequency factor):"<<endl 
 
      <<"A="; 
 
  cin>>A; 
 
  cout<<endl;       
 
         
 
   
 
  //---------------------------------------------------------- 
 
  //end data initiation 
 
  //---------------------------------------------------------- 
 
  //begin calculation 
 
  int i=0; 
 
  time1=0; 
 
  x[i]=1; 
 
  while(time1<=100*(1E6)*365*24*3600) 
 
  { 
 
    i++; 
 
    time2=time1+(1E6)*365*24*3600; 
 
    T[i]=T[i-1]+Tx; 
 
    k=A*exp(-Ea/R/(T[i]+273)); 
 
    x[i]=x[i-1]*exp(-k*(time2-time1)); 
 
    time1=time2; 
 
  } 
 
   
 
   
 
  //end calculation 
 
  //---------------------------------------------------------- 
 
  //begin output the results 
 
    cout<<"Temperature(C)                "<<"X(%)"<<endl; 
 
  for(i=0;i<=250;i++) 
 
    { 
 
     cout<<"   "<<T[i]<<"                 "<<x[i]<<endl; 
 
       if(x[i]<=1E-10) break; 
 
    } 
 
  cout<<"Calculate another source rock? Enter Y to continue or N to end the program"<<endl; 
 
  char loop; 
 
  cin>>loop; 
 
  if(loop=='y'||loop=='Y') 
 
  goto L1; 
 
  else 
 
  { 
 
    cout<<"Press any key to quit out"; 
 
    getch(); 
 
  } 
 
   
 
  //end output the results                           
 
                  
 
}