Commits

patrik  committed 018ed52

Initial revision of the Radiative Transfer code

  • Participants
  • Parent commits e71a17c

Comments (0)

Files changed (10)

File src/Makefile.in

+# Makefile for the utilities library
+
+####### BEGIN
+build_error:
+	@echo "This section should have been replaced by a config file."
+####### END
+
+### Utility directory files ###
+UTIL = ../utilities
+
+INT = $(UTIL)/integrator.h
+COU = $(UTIL)/counter.h
+MAT = $(UTIL)/matrix.h
+LSQ = $(UTIL)/lsqfitter.h
+LS2 = $(UTIL)/lsqfitter2.h
+BIO = $(UTIL)/biniostream.h
+LIS = $(UTIL)/lists.h
+SOL = $(UTIL)/solver.h
+INP = $(UTIL)/interpolator.h
+
+
+####### Compiler, tools and options
+
+CC	=	$(SYSCONF_CC)
+CFLAGS	=	$(SYSCONF_CFLAGS) -I$(UTIL)
+LFLAGS	=	$(SYSCONF_LFLAGS) -L$(UTIL)
+LIBS	=	$(SYSCONF_LIBS)
+MKLIB	= 	$(SYSCONF_MKLIB)
+
+SRCS = *.cc
+
+OBJS = xfer.o optical.o test.o emergence.o
+
+### RULES ###
+%.o : %.cc
+	g++ -c $(CFLAGS) $< -o $@
+
+### TARGETS
+
+depend:
+	makedepend -- $(CFLAGS) -- $(SRCS)
+
+clean:
+	rm *.o
+
+test: test.o $(OBJS)
+	g++ $(LFLAGS) $(LIBS) -lPJutil -o test $(OBJS)
+
+### DEPS
+
+# DO NOT DELETE THIS LINE -- make depend depends on it

File src/emergence.cc

+// emergence.h
+// The class for storing the emerging rays info
+// $Id$
+// $Log$
+
+#include "emergence.h"
+#include "matrix.h"
+#include "counter.h"
+#include <math.h>
+
+const et pi=atan(1)*4;
+
+emergence::emergence()
+{
+  radius=9;
+  n_pt=20;
+  n_pp=20;
+  n_dt=20;
+  n_dp=20;
+
+  // Create histogram
+  h=new histogram(4,-1.,1.,n_pt,-pi,pi,n_pp,-1.,1.,n_dt,-pi,pi,n_dp);
+
+
+  // Temporary projection plane
+  et theta=pi/2;
+  et phi=0;
+  et fov=15*pi/180;
+  xs=100;
+  ys=100;
+
+  lensradius=100; // Distance from center to the lens
+  f=2; // focal length of the lens
+  s=lensradius*f/(lensradius-f); // focus lens
+  R=100; // lens aperture - the lens is so fucking bad at large apertures
+  scale=xs/(s*2*tan(fov/2)); // image scale in pixels/image plane coordinate
+  maxix=xs/(2*scale); // max image plane coords in the image
+  maxiy=ys/(2*scale);
+
+  // Get the lens plane normal and point
+  lpn=matrix<et>(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta));
+  lpp=lpn*lensradius;
+
+  // Define unit axis vectors in the lens plane
+  nx=matrix<et>(-sin(theta)*sin(phi),sin(theta)*cos(phi),0.);
+  ny=matrix<et>(-cos(theta)*cos(phi),-cos(theta)*sin(phi),sin(theta));
+  nx=norm(nx);
+  ny=norm(ny);
+  image=0;
+  
+  // Allocate image
+  image= new double[xs*ys];
+  for(int i=0;i<xs*ys;i++)
+    image[i]=0;
+  maxval=0;
+  
+};
+
+emergence::~emergence()
+{
+  if(h)
+    delete h;
+}
+
+// Adds a ray to the emerging ensemble
+void emergence::add(ray* r)
+{
+  /*
+  // Extrapolate the position of the ray to a sphere exterior
+  // to the simulation volume
+  
+  double a=dot(r->pos,r->dir);
+  double b=dot(r->pos,r->pos)-radius*radius;
+  double l=-a+sqrt(a*a-b); // Positive solution is the one we want by constr.
+  matrix<double> p=r->pos+l*r->dir;
+
+  // Now get spherical angles of this position
+  double p_phi=atan2(p[1],p[0]);
+
+  // And spherical angles of the direction
+  double d_phi=atan2(r->dir[1],r->dir[0]);
+
+  // And add the intoensity of the ray to the bin
+  // cos(theta) is already in the z component of the vectors so we save some
+  // trig there
+  h->addvalue(r->intensity, p[2]/radius, p_phi, double(r->dir[2]), d_phi);
+  */
+
+  // Do the projection shit directly
+  // Find the intercept with the lens plane
+  et l=dot(lpn,lpp-r->pos)/dot(lpn,r->dir);
+  if(l>0) {
+    //cout << "nonzero bin\n";
+    // If it's heading away from the lens we don't want it
+    r->pos+=l*r->dir;
+    
+    // Get the components of position and angle of the ray
+    // in the lens plane
+    et lx=dot(r->pos,nx);
+    et ly=dot(r->pos,ny);
+    if(lx*lx+ly*ly<R*R) {
+      et cosalpha_x=dot(r->dir,nx);
+      et cosalpha_y=dot(r->dir,ny);
+      
+      // Now we are working in cartesian 2-space, puh!
+      // Use the thin lens formula to get the position in the image plane
+      et ix=lx*(1-s/f)+s*cosalpha_x/sqrt(1-cosalpha_x*cosalpha_x);
+      et iy=ly*(1-s/f)+s*cosalpha_y/sqrt(1-cosalpha_y*cosalpha_y);
+					      
+      // Now get the pixel numbers and increase that pixel value	    
+      if( (fabs(ix)<maxix) && (fabs(iy)<maxiy) ) {
+	// It's within the image
+	int px=int(ix*scale+xs/2);
+	int py=int(iy*scale+ys/2);
+	
+	//cout << lx << '\t' << ly << '\t' << cosalpha_x << '\t' << cosalpha_y << '\t' << ix << '\t' << iy << '\t' << px << '\t' << py <<'\n';
+	int i=px+py*xs;
+	image[i]+=r->intensity;
+	
+	if(image[i]>maxval) maxval=image[i];
+      }
+    }
+  }
+  
+};
+
+
+// Creates an image of the emerging rays by projecting them onto a plane
+// normal to (theta,phi). The field of view is fov and the image has
+// xs*ys pixels.
+void emergence::project(et theta, et phi, et fov, int xs, int ys)
+{
+  /*
+  et lensradius=100; // Distance from center to the lens
+  et f=2; // focal length of the lens
+  et s=lensradius*f/(lensradius-f); // focus lens
+  et R=100; // lens aperture - the lens is so fucking bad at large apertures
+  et scale=xs/(s*2*tan(fov/2)); // image scale in pixels/image plane coordinate
+  et maxix=xs/(2*scale); // max image plane coords in the image
+  et maxiy=ys/(2*scale);
+
+  // Get the lens plane normal and point
+  matrix<et> lpn(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta));
+  matrix<et> lpp=lpn*lensradius;
+
+  // Define unit axis vectors in the lens plane
+  matrix<et> nx(-sin(theta)*sin(phi),sin(theta)*cos(phi),0.);
+  matrix<et> ny(-cos(theta)*cos(phi),-cos(theta)*sin(phi),sin(theta));
+  nx=norm(nx);
+  ny=norm(ny);
+
+  // Allocate image
+  double* image= new double[xs*ys];
+  for(int i=0;i<xs*ys;i++)
+    image[i]=0;
+  et maxval=0;
+
+  counter cc(1000);
+  // Loop over all the bins (actually we only need to loop over half of them)
+  for(int a=0;a<n_pt;a++)
+    for(int b=0;b<n_pp;b++)
+      for(int c=0;c<n_dt;c++)
+	for(int d=0;d<n_dp;d++) {
+	  cc.incr();
+	  
+	  // Get the bin value
+	  bin bb = h->getbin(a,b,c,d);
+	  //	  bb.value=1;
+	  if(bb.value>0) {
+
+	    // Now we need to project these onto the image plane
+	    // Make up a ray
+	    // bb->coords[0] contains cos(theta) of the position
+	    // bb->coords[1] phi of the position
+	    ray r;
+	    et sintheta=sin(acos(bb.coords[0]));
+	    r.pos[0]=radius*sintheta*cos(bb.coords[1]);
+	    r.pos[1]=radius*sintheta*sin(bb.coords[1]);
+	    r.pos[2]=radius*bb.coords[0];
+	    // same goes for the direction
+	    sintheta=sin(acos(bb.coords[2]));
+	    r.dir[0]=sintheta*cos(bb.coords[3]);
+	    r.dir[1]=sintheta*sin(bb.coords[3]);
+	    r.dir[2]=bb.coords[2];
+	    
+	    // Find the intercept with the lens plane
+	    et l=dot(lpn,lpp-r.pos)/dot(lpn,r.dir);
+	    if(l>0) {
+	      //cout << "nonzero bin\n";
+	      // If it's heading away from the lens we don't want it
+	      r.pos+=l*r.dir;
+	      
+	      // Get the components of position and angle of the ray
+	      // in the lens plane
+	      et lx=dot(r.pos,nx);
+	      et ly=dot(r.pos,ny);
+	      if(lx*lx+ly*ly<R*R) {
+		et cosalpha_x=dot(r.dir,nx);
+		et cosalpha_y=dot(r.dir,ny);
+
+		// Now we are working in cartesian 2-space, puh!
+		// Use the thin lens formula to get the position in the image plane
+		et ix=lx*(1-s/f)+s/tan(acos(cosalpha_x));
+		et iy=ly*(1-s/f)+s/tan(acos(cosalpha_y));
+		
+		// Now get the pixel numbers and increase that pixel value	    
+		if( (fabs(ix)<maxix) && (fabs(iy)<maxiy) ) {
+		  // It's within the image
+		  int px=int(ix*scale+xs/2);
+		  int py=int(iy*scale+ys/2);
+		  
+		  //cout << lx << '\t' << ly << '\t' << cosalpha_x << '\t' << cosalpha_y << '\t' << ix << '\t' << iy << '\t' << px << '\t' << py <<'\n';
+		  int i=px+py*xs;
+		  image[i]+=bb.value;
+		  
+		  if(image[i]>maxval) maxval=image[i];
+		}
+	      }
+	    }
+	  }
+	}
+  */
+  // Ok, that's that. no worse eh! :)
+  // Now save the image in a ppm file
+
+  ofstream ppmfile("test.ppm");
+  ppmfile << "P2 " << xs << ' ' << ys << ' ' << 255 << '\n';
+
+  // Loop over pixels and output values
+  for(int y=0;y<ys;y++)
+    for(int x=0;x<xs;x++)
+      ppmfile << int( image[x+y*xs]/maxval*255 ) << '\n';
+}

File src/emergence.h

+// emergence.h
+// The class for storing the emerging rays info
+// $Id$
+// $Log$
+
+#ifndef __emergence__
+#define __emergence__
+
+#include "histogram.h"
+#include "ray.h"
+
+typedef double et;
+
+
+// The outgoing rays are propagated to a sphere exterior to the simulation
+// volume. The histogram axes are position cos(theta),phi and direction
+// cos(theta),phi
+// Using cos(theta) instead of straight theta gives the bins uniform surface
+// angle
+
+class emergence {
+private:
+  histogram* h; // This keeps the info about emerging rays
+  int n_pt,n_pp,n_dt,n_dp; // Number of bins
+  et radius; // Radius of recorded sphere
+
+  // projection shit
+  et lensradius;
+  et f;
+  et s;
+  et R;
+  et scale;
+  et maxix;
+  et maxiy;
+  matrix<et> lpn;
+  matrix<et> lpp;
+  matrix<et> nx;
+  matrix<et> ny;
+  double* image;
+  et maxval;
+
+  int xs,ys;
+
+public:
+  emergence();
+  ~emergence();
+
+  void add(ray*); // Adds a ray to the emerging ensemble
+
+  save(char*); // Saves itself to a file
+  load(char*); // Loads itself from a file
+
+  // Projects an image of the object from a specified angle
+  void project(et,et,et,int,int); 
+};
+
+#endif
+// grid.h
+// A class for a grid of varying refinement
+// $Id$
+// $Log$
+
+#ifndef __grid__
+#define __grid__
+
+#include "matrix.h"
+#include <math.h>
+
+typedef double gt;
+
+
+template<class T> class grid_cell;
+  
+template<class T> class grid {
+
+  friend class grid_cell<T>;
+  friend main();
+private:
+  // Grid dimensions
+  matrix<gt> min;
+  matrix<gt> max;
+  matrix<gt> size;
+  matrix<int> n;
+  int ntot;
+
+  // Array of the grid cells in this grid
+  grid_cell<T>* cells;
+  // Parent grid cell, the grid cell in which this grid is located
+  // (for refined grids only)
+  grid_cell<T>* parent;
+
+  // Returns the index of the cells array of a cell at pos x,y,z
+  int index(int,int,int) const;
+  int index(const matrix<int>&) const;
+
+  // enumeration counter
+  int counter;
+
+public:
+  grid();
+  grid(gt,gt,int,gt,gt,int,gt,gt,int);
+  ~grid();
+
+  // Finds cell containing the given point
+  grid_cell<T>* locate(matrix<gt>); 
+
+  // enumerates the grid cells
+  void init();
+  grid_cell<T>* enumerate();
+};
+
+// *** class grid_cell ***
+
+template<class T> class grid_cell {
+
+  friend class grid<T>;
+  friend main();
+private:
+  // Pointer to the grid this is a cell in
+  grid<T>* host;
+
+  // Pointer to a possible refined grid in this cell
+  grid<T>* sub;
+  // Pointer to data for this cell (these two are mutually exclusive)
+  T* celldata;
+
+  // Coordinates for grid corner with lowest coordinates
+  matrix<gt> min;
+
+  // Returns the xyz indices in the grid of this grid cell
+  matrix<int> indices();
+
+public:
+  grid_cell();
+  ~grid_cell();
+
+  T* data();
+
+  // Finds neighbouring cell
+  grid_cell<T>* neighbour(matrix<gt>,int,int);
+
+  // Refine this grid cell with a subgrid
+  void refine(matrix<int>);
+
+  gt xmin() const;
+  gt ymin() const;
+  gt zmin() const;
+  gt xmax() const;
+  gt ymax() const;
+  gt zmax() const;
+  gt xsize() const;
+  gt ysize() const;
+  gt zsize() const;
+};
+ 
+// Creates an empty grid
+template<class T> inline grid<T>::grid()
+{
+  cells=0;
+  parent=0;
+
+  min=matrix<gt>(3,0,0.);
+  max=matrix<gt>(3,0,0.);
+  n=matrix<int>(3,0,0);
+  ntot=0;
+}
+
+template<class T> inline grid<T>::~grid()
+{
+  if(cells)
+    delete[] cells;
+}
+
+// Returns the array index for the (x,y,z)'th cell
+template<class T> inline grid<T>::index(int x, int y, int z) const
+{
+  return x+ n[0]*y + n[0]*n[1]*z;
+}
+
+// Returns the array index for the (x,y,z)'th cell
+template<class T> inline grid<T>::index(const matrix<int>& i) const
+{
+  return i[0]+ n[0]*i[1] + n[0]*n[1]*i[2];
+}
+
+// Constructir creates a celldata along with the grid cell itself
+template<class T> inline grid_cell<T>::grid_cell()
+{
+  sub=0;
+  host=0;
+  celldata = new T;
+}
+
+template<class T> inline grid_cell<T>::~grid_cell()
+{
+  if(sub)
+    delete sub;
+
+  if(celldata)
+    delete celldata;
+}
+
+template<class T> inline T* grid_cell<T>::data()
+{
+  return celldata;
+}
+
+template<class T> inline matrix<int> grid_cell<T>::indices()
+{
+  matrix<int> ind(3,1,0,0,0);
+  int n = int(this - host->cells);
+  
+  int t1=host->n[0]*host->n[1];
+  ind[2]=n/t1;
+  int t2=n%t1;
+  ind[1]=t2/host->n[0];
+  ind[0]=t2%host->n[0];
+
+  return ind;
+}
+
+template<class T> inline gt grid_cell<T>::xmin() const
+{
+  return min[0];
+}
+
+template<class T> inline gt grid_cell<T>::ymin() const
+{
+  return min[1];
+}
+
+template<class T> inline gt grid_cell<T>::zmin() const
+{
+  return min[2];
+}
+
+template<class T> inline gt grid_cell<T>::xmax() const
+{
+  return min[0]+xsize();
+}
+
+template<class T> inline gt grid_cell<T>::ymax() const
+{
+  return min[1]+ysize();
+}
+
+template<class T> inline gt grid_cell<T>::zmax() const
+{
+  return min[2]+zsize();
+}
+
+template<class T> inline gt grid_cell<T>::xsize() const
+{
+  return host->size[0];
+}
+
+template<class T> inline gt grid_cell<T>::ysize() const
+{
+  return host->size[1];
+}
+
+template<class T> inline gt grid_cell<T>::zsize() const
+{
+  return host->size[2];
+}
+
+// Defines a grid of specified positions and number of cells
+template<class T> grid<T>::grid(gt x1, gt x2, int nx, 
+				gt y1, gt y2, int ny, 
+				gt z1, gt z2, int nz)
+{
+  parent=0;
+  min=matrix<gt>(x1,y1,z1);
+  max=matrix<gt>(x2,y2,z2);
+  n=matrix<int>(3,1,nx,ny,nz);
+  ntot=nx*ny*nz;
+  size=matrix<gt>( (x2-x1)/nx, (y2-y1)/ny, (z2-z1)/nz);
+
+  int ntot=nx*ny*nz;
+  cells = new grid_cell<T>[ntot];
+  for(int x=0;x<nx;x++)
+    for(int y=0;y<ny;y++)
+      for(int z=0;z<nz;z++) {
+	int i=index(x,y,z);
+	cells[i].host = this;
+	cells[i].min = matrix<gt>(min[0]+x*size[0],
+				  min[1]+y*size[1],
+				  min[2]+z*size[2]);
+      }
+}
+
+// Find the grid cell that contains position p
+// If this grid is not in a top-level grid, we accept positions outside
+// of the grid itself, and just use the closest coordinates. This is used
+// by neighbour.
+template<class T> grid_cell<T>* grid<T>::locate(matrix<gt> p)
+{
+  int a=int(floor((p[0]-min[0])/size[0]));
+  int b=int(floor((p[1]-min[1])/size[1]));
+  int c=int(floor((p[2]-min[2])/size[2]));
+  matrix<int> ind(3,1,a,b,c);
+
+  // Check for out of bonds
+  int err=0;
+  for(int i=0;i<3;i++) {
+    if(ind[i]<0) {
+      ind[i]=0;
+      err=1;
+    }
+    if(ind[i]>=n[i]) {
+      ind[i]=n[i]-1;
+      err=1;
+    }
+  }
+
+  // If we are in a top-level grid, return null for out-of-bonds
+  if(!parent && err)
+      return 0;
+
+  int i=index(ind);
+	       
+  if(cells[i].sub)
+    // That cell is subdivided, recurse in the locate call
+    return cells[i].sub->locate(p);
+  else
+    return &(cells[i]);
+}
+
+
+// Finds cell neighbour by projecting position p in dimension dim, direction
+// dir (dir<0, lower, dir>=0 upper)
+// if p is not within this cell weird things may or may not happen, but the
+// call doesn't make sense then anyway. (For non-refined grids p doesn't
+// matter.)
+template<class T> grid_cell<T>* grid_cell<T>::neighbour(matrix<gt> p,
+							int dim, int dir)
+{
+  int d=(dir>=0?1:-1);
+
+  // Get indices of this grid cell
+  matrix<int> ind=indices();
+
+  // Get the neighbouring cell
+  ind[dim]+=d;
+
+  // now we have 3 possibilities: normal grid cell, refined grid cell or
+  // outside the grid (in which case it may mean outside grid volume or
+  // in another higher-level cell
+  if( (ind[dim]<0) || (ind[dim]>=host->n[dim]) )
+    // outside grid
+    if(!host->parent)
+      // The grid we're in is the top-level grid, so the point must be
+      // outside the entire grid volume - return null
+      return 0;
+    else {
+      // there is a higher-level grid - find the neighbouring cell in that grid
+      matrix<int> ind2=host->parent->indices();
+      ind2[dim]+=d;
+      grid_cell<T>* upone = 
+	&(host->parent->host->cells[ host->parent->host->index(ind2) ]);
+
+      // Now, if that cell is not refined it's easy, just return it
+      if(!upone->sub)
+	return upone;
+      else
+	// It IS refined - call locate on that subgrid to find the cell
+	// (this relies on a part of locate that should not be used by the
+	// public, but they can't use it since they can never get pointers to
+	// the refined grids anyway) (i hope...)
+	return upone->sub->locate(p);
+    }
+  else {
+    // So it is a valid grid cell in this grid
+    // Get a pointer to it
+    grid_cell<T>* nb = &(host->cells[host->index(ind)]);
+    
+    // If it's not refined it's easy, just return it.
+    if(!nb->sub)
+      return nb;
+    else
+      // It IS refined, so just like above we call locate in that subgrid
+      return nb->sub->locate(p);
+
+  }
+
+}  
+
+// Refines this grid_cell to contain a sub grid of subn cells
+template<class T> void grid_cell<T>::refine(matrix<int> subn)
+{
+  if(!sub) {
+    sub = new grid<T>(xmin(),xmax(),subn[0],ymin(),ymax(),subn[1],
+		   zmin(),zmax(),subn[2]);
+    sub->parent =this;
+  }
+  else
+    cerr << "This cell is already refined, pucko!\n";
+
+  // Delete the celldata
+  if(celldata)
+    delete celldata;
+  celldata=0;
+}
+
+// Recursively reset the counters in the whole grid
+template<class T> void grid<T>::init()
+{
+  counter=0;
+
+  for(int i=0;i<ntot;i++)
+    if(cells[i].sub)
+      cells[i].sub->init();
+}
+
+// successively (and recursively) returns a pointer to all 
+// the grid cells in the grid
+template<class T> grid_cell<T>* grid<T>::enumerate()
+{
+  if(counter>=ntot)
+    return 0;
+
+  grid_cell<T>* p;
+
+  if(cells[counter].sub) {
+    // The current cell is has a subgrid, recurse
+    p=cells[counter].sub->enumerate();
+
+    // If the call returned zero, there are no more in that subgrid and we can
+    // advance to the next in this grid
+    if(!p) {
+      counter++;
+      p=&cells[counter];
+      // advance for next invocation
+      counter++;
+    }
+  }
+  else {
+    // the current cell is a valid cell to return
+    p=&cells[counter];
+    // advance for next invocation
+    counter++;
+  }
+
+  return p;
+}
+
+#endif

File src/optical.cc

+// optical.cc
+// classes for optical data, grid cell data and dust
+// $Id$
+// $Log$
+
+#include "optical.h"
+#include <stdlib.h>
+
+const double pi=atan(1)*4;
+
+dust_grain::dust_grain()
+{
+  kappa_a=0;
+  kappa_s=1;
+}
+
+dust_grain::~dust_grain()
+{}
+
+void dust_grain::scatter(ray* r) const
+{
+  // For now we only have isotropic scattering, so we simply redraw the
+  // random direction
+
+  // Get the emission direction from an isotropic distribution
+  rt phi=drand48()*2*pi;
+  rt theta=acos(2*drand48()-1);
+  r->dir[0]=sin(theta)*cos(phi);
+  r->dir[1]=sin(theta)*sin(phi);
+  r->dir[2]=cos(theta);
+
+  // zero the travelled optical depth
+  r->tau_s =0;
+  r->tau_a =0;
+}
+
+optical_data::optical_data()
+{
+  density=0;
+  dust=0;
+  emissivity=0;
+  em=0;
+  edep=0;
+}
+
+optical_data::~optical_data()
+{}
+

File src/optical.h

+// optical.h
+// classes for optical data, grid cell data and dust
+// $Id$
+// $Log$
+
+#ifndef __optical__
+#define __optical__
+
+typedef double ot;
+
+#include "ray.h"
+
+// *** class dust_grain ***
+// dust optical characteristics
+
+class dust_grain {
+private:
+  // these will be replaced with a wavelength-dependent thing later
+  ot kappa_a; // absorption opacity 
+  ot kappa_s; // scattering opacity
+
+  // we need to define some kind of phase function for scattering as well
+
+public:
+  dust_grain();
+  ~dust_grain();
+
+  ot op_scat(ot lambda) const {return kappa_s;};
+  ot op_abs(ot lambda) const {return kappa_a;};
+
+  void scatter(ray*) const; // scatters a light ray
+};
+
+
+
+// *** class emitter ***
+// emitter spectral characteristics
+
+class emitter {
+private:
+  // here we will put Bruzual-Charlot models or something
+  // but for now never mind, it just returns one.
+
+public:
+  emitter();
+  ~emitter();
+
+  ot emissivity(ot lambda) {return 1;};
+};
+
+
+// *** class optical_data ***
+
+// This class contains the grid cell optical characteristics
+// For now only one species of dust and emissivity is allowed, later we
+// can add lists and then the simple functions must sum over all contributions
+class optical_data {
+private:
+public:
+  // dust density and characteristics
+  // in the future we can augment this to several dust species
+  ot density;
+  dust_grain* dust;
+
+  // emissivity and characteristics
+  // here we will need several emitters as well (possibly one per grid cell
+  // if we pre-run the stellar synthesis models)
+  ot emissivity;
+  emitter* em;
+
+  // energy deposited by absorbed rays
+  ot edep;
+
+public:
+  optical_data();
+  ~optical_data();
+
+  // Returns opacities for this grid cell
+  ot op_abs(ot l) const {return density*dust->op_abs(l);};
+  ot op_scat(ot l) const {return density*dust->op_scat(l);};;
+
+  // Deposits energy in this grid cell
+  void deposit(ot e) {edep+=e;};
+
+  // Scatters a ray against the dust
+  void scatter(ray* r) {dust->scatter(r);};
+};
+
+#endif
+// ray.h
+// A light ray class
+// $Id$
+// $Log$
+
+#ifndef __ray__
+#define __ray__
+
+#include "matrix.h"
+
+typedef double rt;
+
+class ray {
+private:
+
+public:
+  matrix<rt> dir;
+  matrix<rt> pos;
+  rt lambda;
+  rt intensity;
+  // Optical depth travelled by the ray
+  rt tau_s; 
+  rt tau_a; 
+
+  ray();
+  ray(matrix<rt>&,matrix<rt>&,rt,rt);
+  ~ray();
+
+  rt energy(); // Returns enerhy of the ray
+};
+
+inline ray::ray()
+{
+  dir=matrix<rt>(0.,0.,1.);
+  pos=matrix<rt>(0.,0.,0.);
+  lambda=1;
+  intensity=1;
+}
+
+inline ray::ray(matrix<rt>& d, matrix<rt>& o, rt l, rt i)
+{
+  dir=d;
+  pos=o;
+  lambda=l;
+  intensity=i;
+}
+
+inline ray::~ray()
+{}
+
+// The energy of the ray is its intensity divided by its wavelength
+// (in some kind of units)
+inline rt ray::energy()
+{
+  return intensity/lambda;
+}
+
+#endif
+#include "grid.h"
+#include "xfer.h"
+#include "optical.h"
+#include "ray.h"
+#include <iostream.h>
+#include "counter.h"
+#include <stdlib.h>
+
+const double pi=atan(1)*4;
+
+main()
+{
+  srand48((long)998723871);
+  cout << drand48() << '\t' << drand48() << '\n';
+  grid<optical_data> g(-5,5,10,-5,5,10,-5,5,10);
+  grid_cell<optical_data>* c;
+  dust_grain d;
+
+  // Set a dust grain for all cells
+  g.init();
+  while(c=g.enumerate()) {
+    c->data()->dust=&d;
+    c->data()->density=0.5;
+  }
+  c=g.locate(matrix<double>(0.5,3.5,0.5));
+  c->data()->density=1;
+  
+  xfer x(&g);
+  counter cc(1000);
+  cout << "Shooting\n";
+
+  for(int i=0;i<10000;i++) {
+       x.shoot();
+       cc.incr();
+  }
+  cout << "Projecting\n";
+  x.out->project(pi/2,0,15*pi/180,100,100);
+}
+// xfer.cc
+// Radiative transfer stuff
+// $Id$
+// $Log$
+
+#include "xfer.h"
+#include <stdlib.h>
+#include <math.h>
+
+const double pi=atan(1)*4;
+
+xfer::xfer(grid<optical_data>* grid)
+{
+  g=grid;
+  out = new emergence;
+  r=new ray;
+  nray=0;
+}
+
+xfer::~xfer()
+{
+  if(r)
+    delete r;
+  if(out)
+    delete out;
+}
+
+// Propagates the ray from its position to the end of the current grid cell
+void xfer::propagate()
+{
+  const optical_data* curdata= cur->data();
+
+  // Compute the line element to the grid boundary
+  rt dlx,dly,dlz;
+
+  // we need to handle divide by zero here
+  // x-boundary 
+  if(r->dir[0] <0)
+    // heading in negative x-dir
+    dlx=(cur->xmin() - r->pos[0])/r->dir[0];
+  else
+    // heading in positive x-dir
+    dlx=(cur->xmax() - r->pos[0])/r->dir[0];
+
+  // y-boundary 
+  if(r->dir[1] <0)
+    // heading in negative y-dir
+    dly=(cur->ymin() - r->pos[1])/r->dir[1];
+  else
+    // heading in positive y-dir
+    dly=(cur->ymax() - r->pos[1])/r->dir[1];
+
+  // z-boundary 
+  if(r->dir[2] <0)
+    // heading in negative x-dir
+    dlz=(cur->zmin() - r->pos[2])/r->dir[2];
+  else
+    // heading in positive x-dir
+    dlz=(cur->zmax() - r->pos[2])/r->dir[2];
+
+  // Find the *closest* boundary
+  // b contains the dim# of the closest boundary, dl the distance
+  int b;
+  rt dl;
+  if(dlx<dly) {
+    b=0;
+    dl=dlx;
+  }
+  else {
+    b=1;
+    dl=dly;
+  }
+
+  if(dlz<dl) {
+    b=2;
+    dl=dlz;
+  }
+
+  // Check if we have reached optical depth for absorption/scattering
+  rt dtau_a = dl*curdata->op_abs(r->lambda);
+  rt dtau_s = dl*curdata->op_scat(r->lambda);
+
+  // Which will take place first? (if any)
+  // Calculate fractional optical distance to absorption/scattering events
+  rt dist_a = (tau_a - r->tau_a)/dtau_a;
+  rt dist_s = (tau_s - r->tau_s)/dtau_s;
+
+  // 3 things can happen: absorption, scattering or just moving to the edge
+  // of the grid cell. Any of these end the execution of propagate()
+  if( (dist_a < dist_s) && (dist_a<1) )
+    // We have absorption - kill it and deposit energy!
+    absorb();    
+
+  else if( (dist_s<dist_a) && (dist_s<1) ) {
+    // We have scattering
+
+    // Update ray position to pos of scattering event
+    r->pos += dist_s*dl*r->dir;
+
+    // and scatter it
+    scatter();
+  }
+  else {
+    // We just propagate to the edge of the cell
+
+    // Update position
+    r->pos += dl*r->dir;
+
+    // Update cur to point to the next grid cell
+    int dir=int(2/r->dir[b]);
+    cur = cur->neighbour(r->pos,b,dir);
+
+    // Update optical depths
+    r->tau_a += dtau_a;
+    r->tau_s += dtau_s;
+
+    // Check if the ray exited the grid completely
+    if(!cur)
+      // It did, bin it as emergent
+      emerge();
+  }
+}
+
+// Deposit the ray energy into the cell and delete the ray
+void xfer::absorb()
+{
+  cur->data()->deposit( r->energy() );
+
+  // Set cur=0 to mark ray is gone
+  cur=0;
+}
+
+// Scatter the ray at the current position
+void xfer::scatter()
+{
+  // Scatter the ray
+  cur->data()->scatter(r);
+
+  // Generate new interaction optical depths
+  tau_a = -log(drand48());
+  tau_s = -log(drand48());
+}
+
+// Put the ray into the histogram of emerging rays
+void xfer::emerge()
+{
+  out->add(r);
+
+  // Set cur=0 to mark ray has left the building
+  cur=0;
+}
+
+
+// Generates a ray and propagates it to its demise
+void xfer::shoot()
+{
+
+  // Generate a ray
+  emit();
+  //cout << r->pos[0] << '\t' << r->pos[1] << '\t' << r->pos[2] << '\n';
+
+  // Propagate ray until it dies
+  while(cur) {
+    propagate();
+    //cout << r->pos[0] << '\t' << r->pos[1] << '\t' << r->pos[2] << '\n';
+  }
+  nray++;
+  /*
+  cout << r->pos << '\t' << r->dir << '\n';
+  double a=dot(r->pos,r->dir);
+  double b=dot(r->pos,r->pos)-15*15;
+  double l=-a+sqrt(a*a-b);
+  matrix<double> p=r->pos+l*r->dir;
+  cout << l << '\t' << p << '\t' << p.mag() << '\n';
+  */
+}
+
+// Emits a ray according to the emissivity distribution
+void xfer::emit()
+{
+  // Finding a correct distribution for emission position and wavelength
+  // is somewhat tricky.
+  // CHEAT!
+  
+  r->pos[0]=drand48()-0.5;
+  r->pos[1]=drand48()-0.5;
+  r->pos[2]=drand48()-0.5;
+  
+  
+  // Get the emission direction from an isotropic distribution
+  rt phi=drand48()*2*pi;
+  rt theta=acos(2*drand48()-1);
+  r->dir[0]=sin(theta)*cos(phi);
+  r->dir[1]=sin(theta)*sin(phi);
+  r->dir[2]=cos(theta);
+
+  // CHEAT
+  r->lambda=1;
+  r->intensity=1;
+
+  // zero the travelled optical depth
+  r->tau_s =0;
+  r->tau_a =0;
+
+  // Generate new interaction optical depths
+  tau_a = -log(drand48());
+  tau_s = -log(drand48());
+
+  // Get the current grid cell
+  cur = g->locate( r->pos );
+}
+// xfer.h
+// The radiative transfer class
+// Takes care of the operations of ray propagation and scattering
+
+#ifndef __xfer__
+#define __xfer__
+
+#include "grid.h"
+#include "ray.h"
+#include "optical.h"
+#include "emergence.h"
+
+class xfer {
+friend main();
+private:
+  // The Grid
+  grid<optical_data>* g;
+
+  // Ray tracing variables
+  ray* r;
+  grid_cell<optical_data>* cur; // Current grid_cell of the ray
+  rt tau_a; // Optical depth of absorption
+  rt tau_s; // Optical depth of scattering
+  int nray;
+
+  // The emerging rays
+  emergence* out;
+
+  void emit();      // Emits a ray
+  void propagate(); // Propagates the ray to the boundary of the grid cell
+  void absorb();    // Deposits ray energy into the grid cell
+  void scatter();   // Scatters ray
+  void emerge();    // Puts the ray into the histogram of emerging rays
+
+
+public:
+  xfer(grid<optical_data>*);
+  ~xfer();
+
+  // Transfers one ray until it exits the volume or is absorbed
+  void shoot();
+};
+
+#endif