added law 11 into fission pop

main.cpp

@@ -310,7 +310,7 @@ int main(int argc, char* argv[]){
 		std::vector<std::string> topes (n_topes);
 		std::vector<float>    fracs_fuel  (n_topes);
 		std::vector<float>    fracs_water (n_topes);
-		topes[0] = "94239.03c";
+		topes[0] = "92233.03c";
 		topes[1] = "1002.03c";
 		fracs_fuel[0] = 1;
 		fracs_fuel[1] = 1;  
@@ -358,7 +358,7 @@ int main(int argc, char* argv[]){
 
 	whistory hist ( N , geom );
 	hist.set_print_level(4);
-	hist.set_dump_level(3);
+	hist.set_dump_level(1);
 	hist.set_device(0);
 	hist.init();
 	hist.print_xs_data();

pop_fission.cu

@@ -300,6 +300,59 @@ __global__ void pop_fission_kernel(unsigned N, cross_section_data* d_xsdata, par
 			mu  = 2.0*get_rand(&rn)-1.0;
 			phi = 2.0*pi*get_rand(&rn);
 
+		}
+		else if( this_law == 11 ){  // energy-dependent maxwellian
+
+			// get tabulated parameters
+			float a0 = edist_lower.var[0];
+			float a1 = edist_upper.var[0];
+			float b0 = edist_lower.cdf[0];
+			float b1 = edist_upper.cdf[0];
+			float U  = edist_lower.pdf[0];
+			float e0 = edist_lower.erg;
+			float e1 = edist_upper.erg;
+			float  a = 0.0;
+			float  b = 0.0;
+			float  g = 0.0;
+			float  c = 0.0;
+			sampled_E = 99999.0;
+
+			// interpolate T
+			if (e1==e0 | edist_lower.intt==1){  // in top bin, both values are the same
+				a = a0;
+				b = b0;
+			}
+			else if (edist_lower.intt==2){// lin-lin interpolation
+				a  = (a1 - a0)/(e1 - e0) * (this_E - e0) + a0;
+				b  = (b1 - b0)/(e1 - e0) * (this_E - e0) + b0;
+				c  = 1.0 + a*b/8.0;
+				g  = sqrtf( c*c - 1.0 ) + c;
+			}
+			else{
+				printf("dont know what to do!\n");
+			}
+
+			// restriction
+			while (sampled_E > this_E - U){
+
+				// rejection sample
+				rn1 = get_rand(&rn);
+				rn2 = get_rand(&rn);
+				sampled_E = -a*g*logf(rn1);
+				c = (1.0-g)*(1.0-logf(rn1)) - logf(rn2);
+				while ( c*c > b*sampled_E ) {
+					rn1 = get_rand(&rn);
+					rn2 = get_rand(&rn);
+					sampled_E = -a*g*logf(rn1);
+					c = (1.0-g)*(1.0-logf(rn1)) - logf(rn2);
+				}
+
+			}
+
+			// isotropic mu/phi
+			mu  = 2.0*get_rand(&rn)-1.0;
+			phi = 2.0*pi*get_rand(&rn);
+
 		}
 		else{
 			printf("LAW %u NOT HANDLED IN FISSION POP!\n",this_law);

unionize.py

@@ -152,9 +152,13 @@ def _unionize(self):
 				rxn = table.reactions[MT]
 				if hasattr(rxn,"ang_energy_in"):
 					self.MT_E_grid=numpy.union1d(self.MT_E_grid,rxn.ang_energy_in)
-				if hasattr(rxn,"energy_dist") and rxn.energy_dist.law!=3 and rxn.energy_dist.law!=66:
-					#print table.name, MT, "law",rxn.energy_dist.law
-					self.MT_E_grid=numpy.union1d(self.MT_E_grid,rxn.energy_dist.energy_in)
+				if hasattr(rxn,"energy_dist"):# and rxn.energy_dist.law!=3 and rxn.energy_dist.law!=66:
+					if hasattr(rxn.energy_dist,"energy_in"):
+						self.MT_E_grid=numpy.union1d(self.MT_E_grid,rxn.energy_dist.energy_in)
+					if hasattr(rxn.energy_dist,"energya_in"):
+						self.MT_E_grid=numpy.union1d(self.MT_E_grid,rxn.energy_dist.energya_in)
+					if hasattr(rxn.energy_dist,"energyb_in"):
+						self.MT_E_grid=numpy.union1d(self.MT_E_grid,rxn.energy_dist.energyb_in)
 
 		self.num_main_E   = self.MT_E_grid.__len__()
 
@@ -809,8 +813,14 @@ def _get_energy_data(self,row,col):
 		#print MTnum
 
 		# do the cases
-		if hasattr(rxn,"energy_dist") and hasattr(rxn.energy_dist,"energy_in"):
-			# there is no higher lavel angular table, everything is in energy_dist
+		if hasattr(rxn,"energy_dist") and ( hasattr(rxn.energy_dist,"energy_in") or hasattr(rxn.energy_dist,"energya_in")):
+			# unionize a/b for law 11 and set it as energy_in, interpolate a/b values to new grid
+			# just in case they have differrent grids...
+			if hasattr(rxn.energy_dist,"energya_in"):
+				rxn.energy_dist.energy_in	= numpy.union1d(rxn.energy_dist.energya_in,rxn.energy_dist.energyb_in)
+				rxn.energy_dist.a			= numpy.interp( rxn.energy_dist.energy_in, rxn.energy_dist.energya_in, rxn.energy_dist.a )
+				rxn.energy_dist.b			= numpy.interp( rxn.energy_dist.energy_in, rxn.energy_dist.energyb_in, rxn.energy_dist.b )
+			# there is no higher level table, everything is in energy_dist
 			# find where this energy lies on this grid
 			upper_index = next((i for i, x in enumerate(this_E < rxn.energy_dist.energy_in) if x), len(rxn.energy_dist.energy_in))
 			lower_index = upper_index - 1
@@ -878,6 +888,13 @@ def _get_energy_data(self,row,col):
 					upper_cdf = numpy.array([rxn.energy_dist.U])
 					lower_pdf = numpy.array([0])
 					upper_pdf = numpy.array([0])
+				elif hasattr(rxn.energy_dist,"a"):  # e dep maxwellian 
+					lower_var = numpy.array([rxn.energy_dist.a[lower_index]])
+					upper_var = numpy.array([rxn.energy_dist.a[upper_index]])
+					lower_cdf = numpy.array([rxn.energy_dist.b[lower_index]])
+					upper_cdf = numpy.array([rxn.energy_dist.b[upper_index]])
+					lower_pdf = numpy.array([rxn.energy_dist.U])
+					upper_pdf = numpy.array([rxn.energy_dist.U])
 				else:
 					print "UNHANDLED ENERGY DIST CONTENTS"