#include <iostream>
//#include <fstream>
#include <cmath>
#include <cstdlib>
using namespace std;
int main(int argc, char *argv[])
{
  //double radius[7]={5.05,8.85,12.25,29.9,37.1,44.3,51.4};  //in cm
  const double thickness[7]={0.025,0.025,0.025, 0.0285,0.0285,0.0285,0.0285}; //in cm
  const double initialVd[7]={65.,65.,65.,70.,70.,70.,70};
  const double beamOnTemp[7]={-13.,-13.,-13.,-3.,-3.,-3.,5.};
  const double beamOffTemp[7]={20.,20.,20.,20.,20.,20.,20.};
  const double neqfluence[7]={21.5e11, 8.91e11, 5.57e11, 1.66e11, 1.30e11, 1.07e11, 0.90e11};//per fb-1

  const double priodLength[19]={3.,8.,4.,8.,4.,8.,28.,8.,4.,8.,4.,8.,16.,8.,4.,8.,4.,8.,1.}; // in month
  const double coolingOn[19] = {0. ,1.,0. ,1.,0., 1.,0.,1.,0. ,1.,0., 1., 1.,1.,1.,1.,1.,1.,1.};
  const double        fac14[19]={0.,1.,0.,1.,0.,1.,0.,1.24,0.,1.24,0.,1.24,0.,1.24,0.,1.24,0.,1.24,0.};
  const double lumi[19]={0,0.05,0,4.95,0,12.,0,18.,0,55.,0,66.,0,64.,0,80.,0.,100.,0}; 
  const double VSCT[19]={150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.};
  const double VPix[19]={150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.,150.};
  // ofstream fout(argv[3]);  // 出力ファイルを作成する
  const double TR=20.; const double Eg=1.21; const double EI=1.09;
  const double alpha=6.90e-18;
  const double kB=8.617e-5;
  const double A[5]={0.42, 0.10, 0.23, 0.21, 0.04};
  const double tau[5]={ 833.33, 28.47, 2.57, 0.09, 0.01}; // in days
  double TemperCorrection[7];
  double year=0.;
  //if(!fout) {
  //  cout << "Cannot open output file.\n";
  //  return 1;
  //}
// log plot
  //double ion=1.e-6*16*pow((double)iz2,0.9);
  for(int i=0; i < 7; i++){
    double TOp=beamOnTemp[i]+273.15;
    TemperCorrection[i]=pow(TOp/266.15,2.)*exp(0.5*Eg/kB*(1/266.15-1/TOp));
    //cout << "Temperture Correction" << i <<":" << TemperCorrection[i] << endl;
  }
  // 19 periods
  const double Vol=1.;//Volume!
  const int stepnum=1000;
  double step=1./(float)stepnum;
  double Ijlast[7][5]={0.};
  for(int i=0; i<19; i++){
    double Ij[7][5],Isum[7]={0.};
    for(int ndev=0; ndev<7; ndev++){
      for(int j=0; j< 5; j++){
        double temp;
        if(lumi[i] == 0.) temp=beamOffTemp[ndev]+273.15;
        else temp=beamOnTemp[ndev]+273.15;
        double Fluence=neqfluence[ndev];
        if(fac14[i] != 1.) Fluence = Fluence*fac14[i];
        double Theta=exp(EI/kB*(1/293.15-1/temp));
        //cout << "Theta" << i <<":" << "device" << ndev<< ":" << Theta << endl;
        //cout << "DelI" <<":" ; 
        for(int istep=0; istep < stepnum; istep++){
          double DelI;
          //if(ndev < 3) Vop=VPix[i];
          //else Vop=VSCT[i];
          DelI=alpha*1.*A[j]*tau[j]*Fluence*lumi[i]*step;
          //cout << "alpha=" << alpha << " Vop=" << Vop << " tau=" << tau[j] << " Fluence=" << Fluence << " lumi=" << lumi[i] << " ";
          //cout << DelI*1.e3< " "; 
          double dump=exp(-Theta*30.5*priodLength[i]*step/tau[j]);
          DelI=DelI*(1.-dump)/Theta/(30.5*priodLength[i]*step); 
          Ij[ndev][j]=Ijlast[ndev][j]*dump+DelI;
          Ijlast[ndev][j]=Ij[ndev][j];
        }
        //cout << endl;
        Isum[ndev]+=Ij[ndev][j];
      }
    }
    year+=priodLength[i]/12.;
    cout << year << " ";
    for(int ndev=0; ndev<7; ndev++)
        cout << TemperCorrection[ndev]*Isum[ndev]*1.e3 << " ";
    cout << endl;
  }
  //ifstream fin(argv[1]);   // 入力ファイルを開く
  //int degree;
  //double dvalue;
  
  //while(!fin.eof()) {
  //  fin >> degree >> dvalue;
    
  //  cout << degree <<" " << dvalue*dvalue << endl;
  //}

  //fin.close();

  return 0;
}