paper

 float d=0.5;

 float t0=0.0;

 const double angle00=0.0;
 const double angle01=0.0;
 const double angle10=PI/6;
 const double angle11=PI/3;

 b0.x=cos(angle00);
 b0.y=sin(angle00) * cos(angle01);
 b0.z=sin(angle00) * sin(angle01);

 b2.x=cos(angle10);
 b2.y=sin(angle10) * cos(angle11);
 b2.z=sin(angle10) * sin(angle11);

 cout<<b0.x<<" "<<b0.y<<"  "<<b0.z<<"  "<<endl;
 cout<<b2.x<<" "<<b2.y<<"  "<<b2.z<<"  "<<endl;

 float ccos,ssin,ssin2,ttg;
 ccos=vector3::Magnitude(b0+b2)/2;
 ssin=sqrt(1.0-ccos);
 ssin2=2.0*ccos*ssin;
 cout<<ccos<<endl;
 ttg=ssin/ccos;

 float k=0;

 float t=3*PI/2;
 k=ttg*sin(t)/cos(t);

 vector3 b11,b22,b33;
 b11=(b0+b2)/(2.0*ccos*ccos);
 b22=vector3::cross(b0,b2);
 b33=b22*k;
 b1=b11+b33/(ssin2);
 float w;
// w=ssin/(sqrt(k*k+ttg*ttg));
 w=ccos*cos(t);

 float xt,yt,zt;

 xt=(b0.x*(1-t0)*(1-t0)+b1.x*w*2.0*t0*(1-t0)+b2.x*t0*t0)/(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0));
 yt=(b0.y*(1-t0)*(1-t0)+b1.y*w*2.0*t0*(1-t0)+b2.y*t0*t0)/(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0));
 zt=(b0.z*(1-t0)*(1-t0)+b1.z*w*2.0*t0*(1-t0)+b2.z*t0*t0)/(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0));

 float xtDao,ytDao,ztDao;

 xtDao=((t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0)) * (-2.0*(1-t0)*b0.x+2.0*b1.x*w*(1.0-2*t)+2.0*t0*b2.x)-(b0.x*(1-t0)*(1-t0)+b1.x*w*2.0*t0*(1-t0)+b2.x*t0*t0)*(-2.0*(1-t0)+2.0*w*(1.0-2.0*t0)+2.0*t0))/(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0))*(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0));
 ytDao=((t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0)) * (-2.0*(1-t0)*b0.y+2.0*b1.y*w*(1.0-2*t)+2.0*t0*b2.y)-(b0.y*(1-t0)*(1-t0)+b1.y*w*2.0*t0*(1-t0)+b2.y*t0*t0)*(-2.0*(1-t0)+2.0*w*(1.0-2.0*t0)+2.0*t0))/(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0))*(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0));
 ztDao=((t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0)) * (-2.0*(1-t0)*b0.z+2.0*b1.z*w*(1.0-2*t)+2.0*t0*b2.z)-(b0.z*(1-t0)*(1-t0)+b1.z*w*2.0*t0*(1-t0)+b2.z*t0*t0)*(-2.0*(1-t0)+2.0*w*(1.0-2.0*t0)+2.0*t0))/(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0))*(t0*t0+w*2.0*t0*(1-t0)+(1-t0)*(1-t0));

 float a000,b000,c000;

 a000=yt*xtDao-xt*ytDao;
 b000=zt*xtDao-xt*ztDao;
 c000=xtDao*cos(d);

 float a111,b111,c111;

 a111=xtDao;
 b111=(c000/a000)*ytDao;
 c111=ztDao-(b000/a000)*ytDao;

 float a222,b222,c222;

 a222=(c111*c111)/(a111*a111);
 b222=(2.0*b111*c111)/(a111*a111)-(2.0*b000*c000)/(a000*a000);
 c222=(b111*b111)/(a111*a111)+(c000*c000)/(a000*a000)-1.0;

 
 float xx,yy,zz;

 zz=(-b222+(sqrt(b222*b222-4.0*a222*c222)))/(2.0*a222);
 xx=(-b111-zz*c111)/(a111);
 yy=(c000-zz*b000)/(a000);