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

[Geoff]

Please edit your question and format code correctly. There's a button for it.

Answer 1

[Walter Roberson]

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.