French fry Length distribution Monte Carlo

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mymath.h"
#define RAD1 13.0
#define RAD2 10.0
#define RAD3 10.0

double ax1[3], ax2[3], ax3[3];

double co2;
void i_axis()
{

double th1=rand()*1.0/RAND_MAX*3.14159*2;  // 0 to 2pi
double z=rand()*2.0/RAND_MAX-1.0; //-1 to 1
ax1[2]=z;
double r=sqrt(1-z*z);
double c1=cos(th1);
double s1=sin(th1);
ax1[0]=c1*r;
ax1[1]=s1*r;

double n1[3]={ -s1,c1,0};
double n2[3]={ -z*c1,-z*s1,r};


#if 0
printf("Rtest %f\n", V_IPROD3(n1, ax1));
printf("Rtest %f\n", V_IPROD3(n2, ax1));
printf("Rtest %f\n", V_IPROD3(n2, n1));
#endif

double th2=rand()*1.0/RAND_MAX*3.14159*2;  // 0 to 2pi
double c2=cos(th2);
double s2=sin(th2);

ax2[0]= c2*n1[0]+s2*n2[0];
ax2[1]= c2*n1[1]+s2*n2[1];
ax2[2]= c2*n1[2]+s2*n2[2];

//printf("Rtest %f\n", V_IPROD3(ax1,ax2));

V_OPROD3(ax1, ax2,ax3);

//printf("%f %f %f\n", ax3[0],ax3[1],ax3[2]);
// t^2
 co2= 
+ax3[0]*ax3[0]/RAD1/RAD1
+ax3[1]*ax3[1]/RAD2/RAD2
+ax3[2]*ax3[2]/RAD3/RAD3;
}

double slen(double *pos)
{
// x=x0+t ax[0];
// y=y0+t ax[1];
// z=z0+t ax[2];
// x^2/rx^2 + .. -1 =0
//(x0^2 + t^2 ax3[0]^2 +2x0 t ax[0])/Rx^2
//(y0^2 + t^2 ax3[1]^2 +2y0 t ax[1])/Ry^2
//(z0^2 + t^2 ax3[2]^2 +2z0 t ax[2])/Rz^2

//t
double c1=
+2*pos[0]*ax3[0]/RAD1/RAD1
+2*pos[1]*ax3[1]/RAD2/RAD2
+2*pos[2]*ax3[2]/RAD3/RAD3;
//const
double c0= -1
+pos[0]*pos[0]/RAD1/RAD1
+pos[1]*pos[1]/RAD2/RAD2
+pos[2]*pos[2]/RAD3/RAD3;
double det=c1*c1-4*co2*c0;
//printf("C %f %f %f   det %f\n", co2,c1,c0, det);
if(det<0.0) return -1;
else return sqrt(det)/(2*co2);
}

double slen2(int ix, int iy, double dx, double dy)
{
double p0[3];

double smax=-1;
int iix,iiy,k;

for(iix=0;iix<2;iix++){
for(iiy=0;iiy<2;iiy++){
for(k=0;k<3;k++)
  p0[k]=ax1[k]*(ix+iix + dx) + ax2[k]*(iy+iiy+dy);
double s=slen(p0);
if(s>smax)smax=s;
}}
return smax;

}

int hist[20];
int main(int argc, char **argv)
{
int i,s;
for(i=0;i<20;i++)hist[i]=0;

for(s=0;s<100;s++){

i_axis();

double dx=rand()*1.0/RAND_MAX;
double dy=rand()*1.0/RAND_MAX;

int rmax=15;

int ix,iy;

for(ix=-rmax;ix<rmax;ix++){
for(iy=-rmax;iy<rmax;iy++){
double s=slen2(ix,iy,dx,dy);
if(s>0.0)
{
int ix=(int)(s+0.5);
hist[ix]++;
//printf("%f\n", s);
}
}}
}//s
for(i=0;i<20;i++)
if(hist[i]!=0)printf("%f %d\n", i+0.5,hist[i]);

}