program FDCCT_DR
%%#include <math.h>
%%#include <unistd.h>

double DR_CDATA[10000];
assign DR_CDATA to "BMD:FAST:DCCT:CDATA";

double DR_THR;
assign DR_THR to "BMD:FAST:DCCT:THR";

double DR_CURRENT;
assign DR_CURRENT to "BMD:FAST:DCCT:CURRENT";

double DR_RAW;

double DR_OFFSET;
assign DR_OFFSET to "BMD:FAST:DCCT:OFFSET";

double DR_INJRATE;
assign DR_INJRATE to "BMD:FAST:DCCT:INJRATE";

double DR_DECRATE;
assign DR_DECRATE to "BMD:FAST:DCCT:DECRATE";

int DR_EVENT;
assign DR_EVENT to "BMD:FAST:CDATA:EVENT";
int DR_EVENT_PRE;
monitor DR_EVENT;

double DR_MAX;
assign DR_MAX to "BMD:FAST:DCCT:MAX";
double DR_MIN;
assign DR_MIN to "BMD:FAST:DCCT:MIN";

double DR_LOW;
assign DR_LOW to "BMD:FAST:DCCT:LOW";
double DR_EXT;
assign DR_EXT to "BMD:FAST:DCCT:EXT";
int DR_EXTCOUNT;
assign DR_EXTCOUNT to "BMD:FAST:DCCT:EXTCOUNT";



int DR_NORD;
assign DR_NORD to "BMD:FAST:DATA.NORD";

int DR_GETOFF;
assign DR_GETOFF to "BMD:FAST:DCCT:GETOFF";
monitor DR_GETOFF;

int DR_INJ;
assign DR_INJ to "BMD:FAST:DCCT:INJCOUNT";

int DR_FIT;
assign DR_FIT to "BMD:FAST:DCCT:FITCOUNT";

double DR_fitarray[100000];

double DR_rate[10000];
assign DR_rate to "BMD:FAST:DCCT:RATE_A";
double DR_inj[10000];

double NG_ARRAY[10000];
assign NG_ARRAY to "BMD:FAST:DCCT:NG_ARRAY";
int    NG_NUM;


int DR_fitnum;
int DR_injnum;

int i;
int j;
int numelem;
int position;
int counter;
double diff;
double diff2;
double samp_time;
double Hfactor;

%{
static void lsqfit(double Data[], int FSTART, int FEND, double samp_time, double FITN, double *rate);
}%


ss state_slife
{

state init
{
when()
{
for (i=0;i<10000;i++){
DR_fitarray[i]=0.0;
DR_inj[i]=0.0;
DR_rate[i]=0.0;
}
samp_time = 1/5e3;


pvGet(DR_EVENT);
DR_EVENT_PRE=DR_EVENT;
} state calrate
}
state calrate
{
when(DR_EVENT != DR_EVENT_PRE)
{
for (i=0;i<10000;i++){
 NG_ARRAY[i]=0.0;
 }
NG_NUM = 0;
DR_EVENT_PRE = DR_EVENT;
pvGet(DR_CDATA);
pvGet(DR_NORD);
pvGet(DR_THR);
pvGet(DR_OFFSET);
pvGet(DR_LOW);
pvGet(DR_EXT);
DR_EXTCOUNT = 0;
if (DR_NORD > 1) samp_time = 1/(DR_NORD*1.0);
DR_CURRENT = 0;
for (i=0;i<DR_NORD;i++){
DR_inj[i]=0.0;
}
DR_MAX = 0.0;
DR_MIN = 100.0;
for (i=0;i<DR_NORD;i++){
DR_CURRENT = DR_CURRENT + DR_CDATA[i];
if (DR_MAX < DR_CDATA[i]) DR_MAX = DR_CDATA[i];
if (DR_MIN > DR_CDATA[i]) DR_MIN = DR_CDATA[i];
}
DR_CURRENT = DR_CURRENT/DR_NORD;
pvPut(DR_CURRENT);
pvPut(DR_MAX);
pvPut(DR_MIN);
DR_RAW = DR_CURRENT + DR_OFFSET;


position = 0;
numelem=0;
DR_fitnum = 0;
DR_injnum = 0;

counter =0;




%{
DR_fitarray[0]=DR_CDATA[0];
if (DR_fitarray[0]>DR_LOW) NG_ARRAY[0] = 1.0;
else NG_NUM = NG_NUM +1; 
do {
numelem++;
position++;
counter++;
DR_fitarray[counter]=DR_CDATA[position];


diff = DR_fitarray[counter]-DR_fitarray[counter-1];
if ((position+1)< DR_NORD) {
  diff2 = DR_CDATA[position+1] - DR_fitarray[counter-1];
  }
else diff2 = diff;
 
if ((diff >DR_THR) && (diff2 > DR_THR)) {
DR_inj[DR_injnum]=diff;
DR_injnum++;
if (numelem >1){


lsqfit(DR_fitarray,0,numelem,samp_time,numelem,&DR_rate[DR_fitnum]);


NG_ARRAY[DR_fitnum] = 1.0;DR_fitnum++;
numelem=0;
counter=0;
position++;
}
DR_fitarray[counter]=DR_CDATA[position];
}
else if (diff < DR_EXT) {
DR_EXTCOUNT = DR_EXTCOUNT + 1;
if (numelem >2){

numelem = numelem -1;
counter = counter -1;
lsqfit(DR_fitarray,0,numelem,samp_time,numelem,&DR_rate[DR_fitnum]);
DR_fitnum++;
numelem = 0;
counter =0;
if (DR_CDATA[position-1] > DR_LOW) NG_ARRAY[DR_fitnum] = 1.0;
else NG_NUM = NG_NUM+1;
}


DR_fitarray[counter]=DR_CDATA[position];



}
}
while (position < DR_NORD);


DR_INJRATE = 0.0;
DR_DECRATE = 0.0;


if (numelem>1) {
lsqfit(DR_fitarray,0,numelem,samp_time,numelem,&DR_rate[DR_fitnum]);
if (DR_fitarray[position]>DR_LOW) NG_ARRAY[DR_fitnum] = 1.0;
else NG_NUM = NG_NUM +1; 


DR_fitnum++;

}




for (i=0;i<DR_fitnum;i++)
{
if (DR_rate[i] > 0) {
 NG_ARRAY[i] = 0;
 NG_NUM = NG_NUM +1;
 }


DR_DECRATE = DR_DECRATE + DR_rate[i]*NG_ARRAY[i];
}
for (i=0;i<DR_injnum;i++)
{
DR_INJRATE=DR_INJRATE+DR_inj[i];
}


}%
DR_INJ=DR_injnum;
DR_FIT=DR_fitnum;
DR_INJRATE = DR_INJRATE;
if ((DR_fitnum - NG_NUM)<1) {DR_DECRATE= 0.0;}
else
{
DR_DECRATE=DR_DECRATE/(DR_fitnum-NG_NUM);
}

pvPut(DR_INJRATE);
pvPut(DR_DECRATE);
pvPut(DR_INJ);
pvPut(DR_FIT);
pvPut(DR_EXTCOUNT);
pvPut(NG_ARRAY);

pvPut(DR_rate);




} state calrate
when(delay(.1)){
} state calrate


}}

ss state_offsetget
{

state init
{
when()
{
DR_GETOFF = 0;
pvPut(DR_GETOFF);
} state get
}
state get
{
when(DR_GETOFF !=0)
 {
 DR_GETOFF = 0;
 pvPut(DR_GETOFF);
 DR_OFFSET = DR_RAW;
 pvPut(DR_OFFSET); 
 } state get
}

}%{
static void lsqfit(double Data[], int FSTART, int FEND, double samp_time, double FITN, double *rate)
{
double a1_L,a2_L,a3_L,a4_L,a5_L,a6_L, s_x_L,s_x2_L,s_y_L, s_xy_L;
int i;

   s_x_L =0;
   s_x2_L = 0;
   s_y_L = 0;
   s_xy_L = 0; 

  for (i=FSTART;i<FEND; i++)
  {
    s_x_L  += i*samp_time;
    s_x2_L += (i*samp_time)*(i*samp_time);
    s_y_L  += (Data[i]);
    s_xy_L += (i*samp_time)*(Data[i]);
  }
  a1_L = FITN*s_x2_L - s_x_L*s_x_L;
  a2_L = FITN*s_xy_L - s_x_L*s_y_L;
  a3_L = a2_L/a1_L;

  *rate= a3_L;
}
}%