Don't know why my code crashes when j=317. I thought was the matrix dimensions, so tried some numbers. Nothing particular happened.
Every help would be very usefull.
#include <stdio.h>
#include <math.h>
int main(){
int T,i,j,t,c,r[T],Jmax,Z,k;
float v[10],R[10],EBNo[10],trapserv[10],a[10],w[10],Pblock[10];
float L[j+1][j],B[T+1][j+1],b[T+1][1000],P[T+1][j+1],q[j+1],qnorm[j+1];
float sigma,mtimi,mtimic,CV;
float Nmax,g,W,EIother,C,Nown,No,totofcell,sumq;
float Gamma,bita,QA,LA,BA;
float x;
//Data input//
` T=3;`
v[1]=0.45;v[2]=0.3;v[3]=0.8;
R[1]=12.2;R[2]=64.0;R[3]=144.0;
trapserv[1]=0.5;trapserv[2]=0.25;trapserv[3]=0.25;
EBNo[1]=5.0;EBNo[2]=4.0;EBNo[3]=3.0;
Nmax=0.8;g=0.005;
totofcell=0.1;W=3.84;CV=1.0;
EIother=2.0;No=3.98;
//Processing section//
//Calculate c, w[i], r[i] and a[i]//
C=Nmax/g;
for (i=1;i<T+1;i++){
w[i]=(pow(10,EBNo[i]/10)*R[i])/(W*1000+pow(10,EBNo[i]/10)*R[i]);
r[i]=floor((w[i]/g)+0.5);
a[i]=(0.1*trapserv[i])/(v[i]*w[i]);
printf("w[%d]= %f\n",i,w[i]);
printf("r[%d]= %f\n",i,floor(r[i]));
printf("a[%d]= %f\n",i,a[i]);
}
//Initialize b, B and q//
L[0][0]=1;
q[0]=1;
sigma=sqrt(log(1+pow(CV,2)));
mtimi=log(EIother*pow(10,-18))-0.5*log(1+pow(CV,2));
mtimic=mtimi+log(1-Nmax)-log(No*pow(10,-18));
for (t=1;t<T+1;t++){
x=0;
x=Nmax-w[t];
if (x>0){
Gamma=0;
Gamma=0.5+0.5*erf((log(x)-mtimic)/(sigma*sqrt(2)));
b[t][0]=1-Gamma;
}
else b[t][0]=1.0;
B[t][0]=b[t][0]; printf("b[%d][0]= %f\n",t,b[t][0]);
}
//Algorithm for b[t][j] calculation//
for (t=1;t<T+1;t++){
for (k=0;k<180;k++){
x=0;
x=Nmax-k*g-w[t];
if (x>0){
Gamma=0.0;
Gamma=0.5+0.5*erf((log(x)-mtimic)/(sigma*sqrt(2)));
b[t][k]=1.0-Gamma;
}
else b[t][k]=1;
B[t][k]=b[t][k];
}
}
//Basic algorithm for B[t][j], P[t][j], L[c][j], q[j] calculation//
j=0;
sumq=0;
while(B[1][j]<=0.999){
j=j+1;
QA=0;
for (t=1;t<T+1;t++){ //q[j] calculation//
r[t]=floor((w[t]/g)+0.5);
if (j-r[t]>=0) QA=QA+q[j-r[t]]*(1.0-B[t][j-r[t]])*a[t]*r[t];
else QA=QA+0;
}
q[j]=QA/j;
sumq=sumq+q[j];
for (t=1;t<T+1;t++){ //P[t][j] calculation//
if ((j-r[t]>=0)&&(q[j]>0)) P[t][j]=(q[j-r[t]]*(1.0-B[t][j-r[t]])*a[t]*r[t]) / (j*q[j]);
else P[t][j]=0;
}
for (c=0;c<j+1;c++){ //L[c][j] calculation//
LA=0;
for (t=1;t<T+1;t++){
if (c>t) L[c][t]=0;
if (j-r[t]>=0){
if (c-r[t]>=0) LA=LA+P[t][j]*(v[t]*L[c-r[t]][j-r[t]]+(1-v[t])*L[c][j-r[t]]);
else LA=LA+P[t][j]*((1-v[t])*L[c][j-r[t]]);
}
else LA=LA+0;
}
L[c][j]=LA; printf("L[%d][%d]= %f\n",c,j,L[c][j]);
}
for (t=1;t<T+1;t++){ //B[t][j] calculation//
BA=0;
for (c=0;c<j+1;c++){
BA=BA+L[c][j]*b[t][c];
}
B[t][j]=BA;
}
Jmax=j;
}
printf("Jmax= %d\n",Jmax);
//Normalize and print q[j]//
sumq=1+sumq;
for (i=0;i<Jmax+1;i++){
qnorm[i]=q[i]/sumq;
printf("q'[%d]= %f\n",i,qnorm[i]);printf("q[%d]= %f\n",i,q[i]);
}
//Print Pblock[t]//
for (t=1;t<T+1;t++){
for (i=0;i<Jmax+1;i++){
Pblock[t]=Pblock[t]+B[i][t]*qnorm[i];
}
printf("Pblock[%d]= %f\n",t,Pblock[t]);
}
scanf("\n");
}