tiny_mc.c

Drop-Dead Simple Monte Carlo Codes 数个文件啊 好东西

char   t1[80] = "Tiny Monte Carlo by Scott Prahl (http://omlc.ogi.edu)";

char   t2[80] = "1 W Point Source Heating in Infinite Isotropic Scattering Medium";


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define SHELL_MAX  
101

double mu_a = 2;	
		   /* Absorption Coefficient in 1/cm !!non-zero!! */
double mu_s = 20;	
		   /* Reduced Scattering Coefficient in 1/cm */
double microns_per_shell = 50; 
/* Thickness of spherical shells in microns */

long   i, shell, photons = 10000;

double x, y, z, u, v, w, weight;

double albedo, shells_per_mfp, xi1, xi2, t, heat[SHELL_MAX];


int main () 

{
	albedo = mu_s / (mu_s + mu_a);
	
shells_per_mfp = 1e4/microns_per_shell/(mu_a+mu_s);
	
	
for (i = 1; i <= photons; i++)

	{
		x = 0.0; y = 0.0; z = 0.0;			
		/*launch*/  
		u = 0.0; v = 0.0; w = 1.0;		

		weight = 1.0;
		
		
for (;;) 
{
			t = -log((rand()+1.0)/(RAND_MAX+1.0));	
/*move*/
			x += t * u;
		
	y += t * v;
			z += t * w;  

		
	shell=sqrt(x*x+y*y+z*z)*shells_per_mfp;	/*absorb*/
		
	if (shell > SHELL_MAX-1) shell = SHELL_MAX-1;	
		
	heat[shell] += (1.0-albedo)*weight;
			
weight *= albedo;
		
			
for(;;) 
{								/*new direction*/
		
		xi1=2.0*rand()/RAND_MAX - 1.0; 
				
xi2=2.0*rand()/RAND_MAX - 1.0; 
				
if ((t=xi1*xi1+xi2*xi2)<=1) break;
			}
		
	u = 2.0 * t - 1.0;
			
v = xi1 * sqrt((1-u*u)/t);
			
w = xi2 * sqrt((1-u*u)/t);
			
	
		if (weight < 0.001)
{ 					/*roulette*/
		
		if (rand() > 0.1 * RAND_MAX) break; 
		
		weight /= 0.1;
			}
	
	} 
	}	
	
	
printf("%s\n%s\n\nScattering = %8.3f/cm\nAbsorption = %8.3f/cm\n",t1,t2,mu_s,mu_a);
	
printf("Photons    = %8ld\n\n Radius         Heat\n[microns]     [W/cm^3]\n",photons);
	
t = 4*3.14159*pow(microns_per_shell,3)*photons/1e12;
	
for (i=0;i<SHELL_MAX-1;i++)
		
printf("%6.0f    %12.5f\n",i*microns_per_shell, heat[i]/t/(i*i+i+1.0/3.0));

	printf(" extra    %12.5f\n",heat[SHELL_MAX-1]/photons);
	
return 0;
}