ratsphere.txt

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#define PI 3.1415926
#define dbl double
dbl fn[100005];
dbl save[100005];
int main(){
	dbl rad;
	int i,j;
	dbl count = 0.00,x,y,z;
	dbl rats[3];
	dbl cats[3],fl=0;
rep:
	fprintf(stderr,"R??");
	scanf("%lf",&rad);
	dbl range = rad/10;
	srand(time(0));
	for(i=0;i<10000;i++){
		rats[0] = -rad;
		rats[1] = rats[2] = cats[1] = 0;	
		cats[0] = rad;
		cats[2] = 1;	
		fprintf(stderr,"%d\n",i);
		for(j=1;j<=10000;j++){
			int flag=0;
			do{
				dbl theta = (dbl)rand()/(dbl)RAND_MAX;
				theta*=PI;
				dbl phi = (dbl)rand()/(dbl)RAND_MAX;
				phi*=2*PI;
				x = cos(theta);
				y = sin(theta) * cos(phi);
				z = sin(theta) * sin(phi);
				if(sqrt(pow(rats[0]+x,2) + pow(rats[1]+y,2) + pow(rats[2]+z,2)) <= rad){
					flag=1;
					rats[0]+=x;
					rats[1]+=y;
					rats[2]+=z;
				}
			}
			while(flag==0);
			flag=0;
			do{
				dbl theta = (dbl)rand()/(dbl)RAND_MAX;
				theta*=PI;
				dbl phi = (dbl)rand()/(dbl)RAND_MAX;
				phi*=2*PI;
				x = cos(theta);
				y = sin(theta) * cos(phi);
				z = sin(theta) * sin(phi);
				if(sqrt(pow(cats[0]+x,2) + pow(cats[1]+y,2) + pow(cats[2]+z,2)) <= rad){
					flag=1;
					cats[0]+=x;
					cats[1]+=y;
					cats[2]+=z;
				}
			}
			while(flag==0);
			if(sqrt(pow(cats[0]-rats[0],2) + pow(cats[1]-rats[1],2) +pow(cats[3]-rats[3],2)) < range) break; 
			else if(fl==0) fn[j]+=1.0;
			else save[j]+=1.0;
		}
	}
	if(fl==0){
		fl=1;
		goto rep;
	}
	for(i=1;i<=10000;i++){
		printf("%d %lf %lf\n",i,(dbl)(fn[i]/10000.00),(dbl)(save[i]/10000.00));
	}
	return 0;
}