Boolean HitTumBoundary(PhotonStruct * Photon_Ptr,
InputStruct * In_Ptr)
{
short layer = Photon_Ptr->layer;
double ux = Photon_Ptr->ux;
double uy = Photon_Ptr->uy;
double uz = Photon_Ptr->uz;
double x0 = Photon_Ptr->x;
double y0 = Photon_Ptr->y;
double z0 = Photon_Ptr->z;
double x1 = 0;
double y1 = 0;
double z1 = Photon_Ptr->z;
//double z1 = In_Ptr->tz;
double rt = In_Ptr->rtum;
double a, b, c, mut;
double dl_b;
double L1,L2;
double t, lr, zz1, zz2;
double zb1 = In_Ptr->layerspecs[In_Ptr->ntum].z0;
double zb2 = In_Ptr->tz;
double rr = sqrt(x0*x0+y0*y0);
//double mut;
Boolean hit;
/* Distance to the boundary. */
t=((x1-x0)*ux+(y1-y0)*uy); //the length of vector(x1-x0),(y1-y0),(z1-z0) projects onto the photon packet trajectory
lr=sqrt((x0-x1)*(x0-x1)+2*(x0-x1)*ux*t+ux*ux*t*t //distance from the tumor center to the photon packet trajectory
+(y0-y1)*(y0-y1)+2*(y0-y1)*uy*t+uy*uy*t*t
+(z0-z1)*(z0-z1)+2*(z0-z1)*uz*t+uz*uz*t*t);
a=ux*ux+uy*uy; //(x+L*ux)^2+(y+L*uy)^2=rsb^2
b=2*x0*ux+2*y0*uy;
c=x0*x0+y0*y0-rt*rt;
//u1if(a!=1.0)printf("a=%lf\n",a);//yk 20130225
if(a!=0 && b*b-4*a*c>=0){
L1=(-b+sqrt(b*b-4*a*c))/(2*a); //length of the photon position to the boundary
L2=(-b-sqrt(b*b-4*a*c))/(2*a); //length of the photon position to the boundary
}else{
L1=0; //length of the photon position to the boundary
L2=0; //length of the photon position to the boundary
}
zz1 = z0+L1*uz;
zz2 = z0+L2*uz;
if(Photon_Ptr->layer == In_Ptr->num_layers+1 && lr<rt && lr>0 && uz != 0.0 && Photon_Ptr->s > L1 && zz1>In_Ptr->layerspecs[In_Ptr->ntum].z0 && L1>0){//在tumor內
mut = In_Ptr->tmua //yk 20130223
+ In_Ptr->tmus;
Photon_Ptr->sleft = (Photon_Ptr->s - L1)*mut;
Photon_Ptr->s = L1;
if(L1<0) printf("L1=%lf\n",L1);//yk 20130225
//if(Photon_Ptr->s<0) printf("s1=%lf\n",Photon_Ptr->s);//yk 20130225
Photon_Ptr->hitsideboundary = 1;
hit = 1;
//printf("hit1=%d\n",hit);//yk 20130305
}
else if(Photon_Ptr->layer == In_Ptr->ntum && lr<rt && lr>0 &&uz != 0.0 && Photon_Ptr->s > L2 && zz2<In_Ptr->tz && zz2>In_Ptr->layerspecs[In_Ptr->ntum].z0 && L2>0){//在tumor外
mut = In_Ptr->layerspecs[layer].mua //yk 20130223
+ In_Ptr->layerspecs[layer].mus;
Photon_Ptr->sleft = (Photon_Ptr->s - L2)*mut;
Photon_Ptr->s = L2;
if(L2<0) printf("L2=%lf\n",L2);//yk 20130225
//if(Photon_Ptr->s<0) printf("s2=%lf\n",Photon_Ptr->s);//yk 20130225
Photon_Ptr->hitsideboundary = 1;
hit = 1;
}
else{
hit = 0;
Photon_Ptr->hitsideboundary = 0;
}
return(hit);
}