program ardcct_l2
%{extern double log();}%
%{extern double fabs();
extern double exp();}%
float FB;
float LIFE;
float RATE;
%{
float ib[500];
float s_xy,s_x,s_y,s_x2,a1,a2,a3,a4,a5,a6;
float samp_time = 1.0;
int l_10, l_60, l_240;
int llflag1,llflag2,llflag3,llflag4;
}%
int i;
int MAXFIT;
int FITSTART;
int FITN;
float P_MON01;
float P_MON01P;
float MIN;
float MIN2;
float ZERO;
assign FB to "fb_ar_dcct:dcct";
assign LIFE to "fb_ar_dcct:life";
assign RATE to "fb_ar_dcct:rate";
assign FITSTART to "fb_ar_dcct:fitstart";
assign P_MON01 to "fb_ar_dcct:dcct";
monitor P_MON01;
ss state_main1
{
state init
{
when()
{
printf("init state\n");
ZERO = 0.0;
MAXFIT = 500;
LIFE = 0.001;
FITSTART = MAXFIT -10;
FITN = MAXFIT -FITSTART;
MIN=0.2;
MIN2=0.003;
P_MON01P=0;
samp_time = 1.0;
llflag1 = 0;
llflag2 = 0;
llflag3 = 0;
llflag4 = 1;
l_10 = 0;
l_60 = 0;
l_240 = 0;
for (i=0; i<200; i++)
{
ib[i]=0.01;
}
}state sta1
}
state sta1
{
when(P_MON01!=P_MON01P)
{
s_x = 0.0;
s_x2 = 0.0;
s_y = 0.0;
s_xy = 0.0;
pvGet(FB);
/* printf("FB= %f\n",FB); */
if (FB<MIN)
ib[MAXFIT-1]=0.01;
else
ib[MAXFIT-1] = FB;
if ((LIFE > 30.0) && (llflag3))
{
if ((l_240 > 239) && (!llflag4))
{
FITSTART=MAXFIT-240;
pvPut(FITSTART);
llflag4=1;
l_240 = 240;
}
else
{
l_240++;
}
}
else if (LIFE > 15.0)
{
if ((l_60 > 59) && (!llflag3) || (llflag4))
{
FITSTART=MAXFIT-60;
pvPut(FITSTART);
llflag3 = 1;
llflag4 = 0;
l_60 = 60;
l_240 = 60;
}
else
{
l_60++;
}
}
else if (LIFE> 1.0)
{
llflag3 = 0;
llflag4 = 0;
l_60 = 0;
l_240 = 0;
FITSTART=MAXFIT-10;
pvPut(FITSTART);
}
else if (LIFE> 0.0)
{
llflag3 = 0;
llflag4 = 0;
FITSTART = MAXFIT - 5;
pvPut(FITSTART);
l_60 = 0;
l_240 = 0;
}
else
{
llflag3 = 0;
llflag4 = 0;
l_60 = 0;
l_240 = 0;
FITSTART = MAXFIT - 3;
pvPut(FITSTART);
}
for (i=FITSTART;i<MAXFIT; i++)
{
s_x += i*samp_time;
s_x2 += (i*samp_time)*(i*samp_time);
s_y += log(ib[i]);
s_xy += (i*samp_time)*log(ib[i]);
}
FITN=MAXFIT-FITSTART;
a1 = FITN*s_x2 - s_x*s_x;
a2 = FITN*s_xy - s_x*s_y;
a3 = a2/a1;
RATE = a3*ib[FITSTART];
pvPut(RATE);
a3 = -a3*60;
/* printf(" a3 and ABS =%f %f\n",a3,fabs(a3)); */
if (fabs(a3)<MIN2)
a4 = 0;
else
a4 =1.0/a3;
LIFE = a4;
pvPut(LIFE);
P_MON01P=P_MON01;
for (i=0;i<MAXFIT-1; i++)
{
ib[i] = ib[i+1];
}
}state sta1
}
}