//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop
#include <stdio.h>
#include <math.h>
#include "Sim.h"
#include "Data.h"
//#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TSimForm *SimForm;
//---------------------------------------------------------------------------
__fastcall TSimForm::TSimForm(TComponent* Owner)
  : TForm(Owner)
{
  SepUnits = true;
  Restricted = false;
  Lake = false;
  Puddle = false;
  BipModel = true;
  RateModel = false;
  DataForm->RateModelBtnClick(NULL);
  TwoStep = true;
  FiveStep = false;
  DataForm->TwoStepTrBtnClick(this);
  TauQAPQ  = 600;
  TauPQ    = 10000;
  TauQAQB  = 160;
  TauQAQB2 = 800;
  TauQB2PQ = 100;
  TauAss   = 600;
  NoRC     = 100000;
  PQPool   = 10;
  kp   = 75;    //per ns
  kl   = 23;
  kf   = 2;
  kto  = 20;
  kbo  = 5;
  ktc  = 4;
  kbc  = 4;
  klrc = 1;
  SimTime = 50000;
  SimStep = 1.0;
  SimRec  = 1.0;
  SimRuns = 10;
  TimeMax = 100000 - 1;
}
//---------------------------------------------------------------------------


void __fastcall TSimForm::SimRun(int NoRun)
{
  int rn, dest;
  bool  QA_Reoxidize_Pending = false;
  for(int i = 0; i<1000000; i++) ev_ptr[i] = 0;
  for(int i = 0; i<1000000; i++) ev_rcnext[i] = 0;
  evfirst = 1;
  evlast = 1;
  evcurr = 1;
  unsigned char op;
  float tsim, quantas, quantaleft, tsimnext, ll;
  sig = DataForm->sig;
  quantaleft = 0;
  float kallopen = kp + kl + kf;
  float kallclos = kl + kf;
  int kpopen = int(100000 * kp/kallopen);
  int kfopen = int(100000 * kf/kallopen);
  int klopen = int(100000 * kl/kallopen);
  int kfclos = int(100000 * kf/kallclos);
  int klclos = int(100000 * kl/kallclos);
  float fyield_open = kf/kallopen;
  float fyield_closed = kf/kallclos;
  int ranindex = random(100000);
  //if (ranindex > 99999) ranindex = 99999;
  int tauqadec = int(100000 * SimStep/TauQAPQ);
  redox = 0;
  pqredox = 0;
  float runweigth = (float(NoRun) - 1)/SimRuns;
  for(int i = 0; i < NoRC; i++)
  {
    QA_reduced[i] = false;
    QA_reoxidation_pending[i] = false;
    PQ_reduced[i] = 0;
  }
  delete(DataForm->FIF);
  DataForm->FIF = new TFastLineSeries(DataForm);
  DataForm->FIF->SeriesColor = clRed;
  DataForm->FIF->ParentChart = DataForm->SimChart;
  delete(DataForm->FIP);
  DataForm->FIP = new TFastLineSeries(DataForm);
  DataForm->FIP->SeriesColor = clBlack;
  DataForm->FIP->ParentChart = DataForm->SimChart;
  delete(DataForm->RED);
  DataForm->RED = new TFastLineSeries(DataForm);
  DataForm->RED->SeriesColor = clBlue;
  DataForm->RED->ParentChart = DataForm->SimChart;



  //DataForm->FIF->Clear();
  //DataForm->RED->Clear();
  DataForm->PQRED->Clear();

  int is, ir;
  evvalu[1] = 0;
  for(int iev = 0; iev < noev; iev ++)
  {
    float quantas = evvalu[iev] * 3.3 * 1.e-8 * sig * NoRC * SimStep;
    tsim = evtime[iev];
    tsimnext = evtime[iev + 1];
    do
    {
      is = int(tsim);
      quantaleft += quantas;
      while (quantaleft > 1)
      {
        exs[is] += 1;
        ir = random(NoRC);
        rn = random(100000);
        if(!QA_reduced[ir])                           // RC open
        {
          if(rn < kpopen)
          {
            QA_reduced[ir] = true;                    // reduce QA, close RC
            redox += 1;
            if(PQ_reduced[ir] < PQPool)
            {
              ll = 0.00001 * random(100000) + 0.00001;
              dest = is - int(TauQAPQ * log(ll));
              if(dest > TimeMax) dest = TimeMax;
              ev_rc[evlast] = ir;
              ev_op[evlast] = 1;
              ev_rcnext[evlast] = ev_ptr[dest];
              ev_ptr[dest] = evlast;
              evlast ++;
              QA_reoxidation_pending[ir] = false;
              if(evlast == 1000000) evlast = 1;
            }
            else QA_reoxidation_pending[ir] = true;   // PQ is fully reduced
          }
          else if((rn >= kpopen) && (rn < kpopen + kfopen)) fif[is] += 1;
          else fih[is] += 1;
        }
        else
        {
          if(rn < kfclos) fif[is] += 1;
          else fih[is] +=1;
        }
        quantaleft -= 1;
      }
//---- Service the scheduled events ------
      while(!(ev_ptr[is] == 0))
      {
        evcurr = ev_ptr[is];
        ir = ev_rc[evcurr];
        op = ev_op[evcurr];
        if(op == 2)                                  // Reoxidize PQ Pool
        {
          PQ_reduced[ir] --;
          pqredox --;
          fip[is] ++;
          if(QA_reoxidation_pending[ir])
          {
            QA_reoxidation_pending[ir] = false;
            ll = 0.00001 * random(100000) + 0.00001;
            dest = is - int(TauQAPQ * log(ll));     // Only now Shedule QA reoxidation
            if(dest > TimeMax) dest = TimeMax;      // PQ was fully reduced
            ev_rc[evlast] = ir;
            ev_op[evlast] = 1;                      // op for QA reoxidation
            ev_rcnext[evlast] = ev_ptr[dest];
            ev_ptr[dest] = evlast;
            evlast ++;
            if(evlast == 1000000) evlast = 1;
          }
        }
        if((op == 1) && (PQ_reduced[ir] < PQPool))  // Move electron from QA to PQ
        {
          QA_reduced[ir] = false;
          ll = 0.00001 * random(100000) + 0.00001;
          dest = is - int(TauPQ * log(ll));         // Shedule PQPool reoxidation
          if(dest > TimeMax) dest = TimeMax;
          ev_rc[evlast] = ir;
          ev_op[evlast] = 2;                        // op for PQPool reoxidation
          ev_rcnext[evlast] = ev_ptr[dest];
          ev_ptr[dest] = evlast;
          evlast ++;
          if(evlast == 1000000) evlast = 1;
          PQ_reduced[ir] ++;
          pqredox ++;
          redox -= 1;
        }
        ev_ptr[is] = ev_rcnext[evcurr];
      }
//--- Done with scheduled events --------
      red[is] = runweigth * red[is] + (1 - runweigth) * redox/NoRC;
      pqred[is] = runweigth * pqred[is] + (1 - runweigth) * pqredox/NoRC;
        //if(exs[is] > 0) DataForm->EXS->AddXY(is, exs[is]* 0.01, ' ', clBlue);
        //if(ems[is] > 0)

        //DataForm->FIF->AddXY(is, fif[is]/(exs[is] + 0.0001) * 1000, ' ', clRed);
        //DataForm->FIF->AddXY(is, fif[is]/(quantas * NoRun + 0.001) * 1000, ' ', clRed);
      if(fmod(is,10) == 0)
      {
        DataForm->RED->AddXY(is, red[is] * 100, ' ', clBlue);
        DataForm->PQRED->AddXY(is, pqred[is] * 100, ' ', clGreen);
      }
      tsim += SimStep;
    }while (tsim < tsimnext);
  }
  for(int idata = 0; idata < DataForm->nodata; idata ++)
  {
    float redox = red[int(tim[idata])];
    f[idata] = (1 - redox) * fyield_open + redox * fyield_closed;
    f[idata] *= 1000;
    DataForm->FIF->AddXY(DataForm->tim[idata], f[idata], ' ', clRed);
  }
  //float fmax = DataForm->FIF->MaxYValue();
  DataForm->SimChart->LeftAxis->SetMinMax(-2, 100);
  DataForm->SimChart->Update();
  AnsiString cc = "";
  char vv[10];
  sprintf(vv, "%10i", evlast);
  //MainForm->Memo->Lines->Add(vv);
}

void __fastcall TSimForm::Simulate()
{
  int i;
  for(i = 0; i < 100000; i++)
  {
    exs[i] = 0;
    fif[i] = 0;
    fip[i] = 0;
    fih[i] = 0;
    red[i] = 0;
    pqred[i] = 0;
  }
  randomize();
  DataForm->fo = 1000 * kf/(kf + kl + kp);
  DataForm->fm = 1000 * kf/(kf + kl);
  sprintf(text, "%10.2f",DataForm->fo);
  DataForm->FoEd->SetTextBuf(text);
  DataForm->FoEd->Update();
  sprintf(text, "%10.2f",DataForm->fm);
  DataForm->FmEd->SetTextBuf(text);
  DataForm->FmEd->Update();
  for(int NoRun = 0; NoRun < SimRuns; NoRun ++)
  {
    SimRun(NoRun + 1);
    sprintf(text, "%d", NoRun + 1);
    DataForm->RunEd->SetTextBuf(text);
    DataForm->RunEd->Update();
  }
  DataForm->F->Clear();
  for(int idata = 0; idata < DataForm->nodata; idata ++)
  {
    DataForm->f[idata] = f[idata];
    DataForm->F->AddXY(DataForm->tim[idata], DataForm->f[idata], ' ', clRed);
  }
  DataForm->FitChart->LeftAxis->SetMinMax(-2, 100);
}