translating code from C++ to matlab
显示 更早的评论
Hi. I seem to be having a problem translating a code from c++ to matlab. Its a simple monte carlo simulation and there seems to be an undefined variable/function in 'heat[]'. I tried translating everything as is and know that heat is undefined and yet the program runs swimmingly in c++. Can anyone help me out with this? I am trying to turn heat[] into an array of collected data points resulting and Id like matlab to print these but juts somehow cant. I even tried defining heat in the beginning and see if results would come out, but it kept giving me the error '??? Subscript indices must either be real positive integers or logicals'...
I am first going to post the C++ code and then I will post my translated version that cant run. PLEASE SOMEONE HELP
C++ version that does run:
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;
}
my crappy rendition
intro = 'Tiny Monte Carlo retranslated into Matlab by Jon';
intro_2= ' 1 W Point Source Heating in Infinite Isotropic Scattering Medium';
shell_max = 101;
mu_a = 2; mu_s = 20;
microns_per_shell = 50;
photons = 10000;
j = 1; h = 1;
max = 32767;
albedo = mu_s/ (mu_s + mu_a);
shells_per_mfp = 1e4/microns_per_shell/ (mu_a+mu_s);
for i = 1:(photons- 1)
i= i+1;
x = 0.0; y = 0.0; z = 0.0;
u = 0.0; v = 0.0; w = 1.0;
weight = 1.0;
while j
t = -log((randi(max)+1.0)/(max+1));
x = x + t*u;
y = y + t*v;
z = z + t*w;
shell = ((abs(x)+abs(y)+abs(z))*shells_per_mfp);
if shell > (shell_max -1.0)
shell = shell_max - 1.0; end
heat(shell)= heat(shell) + (1.0-albedo)*weight;
weight = weight* albedo;
while h
xi1 = (2.0 * randi(max))/ (max - 1.0);
xi2 = (2.0 * randi(max))/ (max - 1.0);
t = (xi1 * xi1)+ (xi2 *xi2);
if t <=1
h=0;
end
end % while h line 33
u = 2.0 * t- 1.0;
v = xi1 * sqrt ((1 - u*u)/t);
w = xi2 * sqrt ((1 - u*u)/t);
if (weight < 0.0001)
if randi > 0.1 * max
end
weight= weight/0.1;
j= 0;
end
end
end
sprintf('%s\n%s\n\nScattering = %8.3f/cm\nAbsorption = %8.3f/cm\n', intro, intro_2, mu_s, mu_a)
sprintf('Photons = %81d\n\n Radius Heat\n[microns] [W/cm^3]\n', photons)
t = 4* 3.141592654 * microns_per_shell^3 * photons/1e12;
for i=0:(shell_max -1)
sprintf ('%6.0f %12.5f\n', i * microns_per_shell, heat(i)/t/((i*i + i + 1.0)/3));
end;
sprintf( ' extra %12.5f\n', heat(shell_max-1)/photons);
return
1 个评论
Geoff
2012-5-24
Please edit your question and format code correctly. There's a button for it.
回答（1 个）
Walter Roberson
2012-5-24
0 个投票
The variable "shell" is declared as being of type "long" in the C code. Any floating point expression assigned to it would automatically be converted to integer as if MATLAB's fix() was called (C doesn't have a specific routine with the same functionality.)
You need to put in fix() in various places.
类别
在 帮助中心 和 File Exchange 中查找有关 MATLAB 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
