diff --git a/src/simulation/Simulation.cpp b/src/simulation/Simulation.cpp
index 1010a200b5..b2739ce151 100644
--- a/src/simulation/Simulation.cpp
+++ b/src/simulation/Simulation.cpp
@@ -1766,27 +1766,42 @@ inline int Simulation::is_wire_off(int x, int y)
 
 int Simulation::get_wavelength_bin(int *wm)
 {
-	int i, w0=30, wM=0;
+	int i, w0, wM, r;
 
-	if (!*wm)
+	if (!(*wm & 0x3FFFFFFF))
 		return -1;
 
+#ifdef __GNUC__
+	w0 = __builtin_ctz(*wm | 0xC0000000);
+	wM = 31 - __builtin_clz(*wm & 0x3FFFFFFF);
+#else
 	for (i=0; i<30; i++)
-		if (*wm & (1<<i)) {
+		if (*wm & (1<<i))
+		{
 			if (i < w0)
 				w0 = i;
 			if (i > wM)
 				wM = i;
 		}
+#endif
 
-	if (wM-w0 < 5)
-		return (wM+w0)/2;
+	if (wM - w0 < 5)
+		return wM + w0;
 
-	i = rand() % (wM-w0-3);
+	r = rand();
+	i = (r >> 1) % (wM-w0-4);
 	i += w0;
 
-	*wm &= 0x1F << i;
-	return i + 2;
+	if (r & 1)
+	{
+		*wm &= 0x1F << i;
+		return (i + 2) * 2;
+	}
+	else
+	{
+		*wm &= 0xF << i;
+		return (i + 2) * 2 - 1;
+	}
 }
 
 void Simulation::set_emap(int x, int y)
@@ -2072,7 +2087,7 @@ void Simulation::init_can_move()
 		 || destinationType == PT_CLNE || destinationType == PT_PCLN || destinationType == PT_BCLN || destinationType == PT_PBCN
 		 || destinationType == PT_WATR || destinationType == PT_DSTW || destinationType == PT_SLTW || destinationType == PT_GLOW
 		 || destinationType == PT_ISOZ || destinationType == PT_ISZS || destinationType == PT_QRTZ || destinationType == PT_PQRT
-		 || destinationType == PT_H2)
+		 || destinationType == PT_H2   || destinationType == PT_BGLA || destinationType == PT_C5)
 			can_move[PT_PHOT][destinationType] = 2;
 		if (destinationType != PT_DMND && destinationType != PT_INSL && destinationType != PT_VOID && destinationType != PT_PVOD && destinationType != PT_VIBR && destinationType != PT_BVBR && destinationType != PT_PRTI && destinationType != PT_PRTO)
 		{
@@ -2259,6 +2274,34 @@ int Simulation::try_move(int i, int x, int y, int nx, int ny)
 			}
 			else if ((r&0xFF) == PT_FILT)
 				parts[i].ctype = Element_FILT::interactWavelengths(&parts[r>>8], parts[i].ctype);
+			else if ((r&0xFF) == PT_C5)
+			{
+				if (parts[r>>8].life > 0 && (parts[r>>8].ctype & parts[i].ctype & 0xFFFFFFC0))
+				{
+					float vx = ((parts[r>>8].tmp << 16) >> 16) / 255.0f;
+					float vy = (parts[r>>8].tmp >> 16) / 255.0f;
+					float vn = parts[i].vx * parts[i].vx + parts[i].vy * parts[i].vy;
+					parts[i].ctype = (parts[r>>8].ctype & parts[i].ctype) >> 6;
+					parts[r>>8].life = 0;
+					parts[r>>8].ctype = 0;
+					// add momentum of photons to each other
+					parts[i].vx += vx;
+					parts[i].vy += vy;
+					// normalize velocity to original value
+					vn /= parts[i].vx * parts[i].vx + parts[i].vy * parts[i].vy;
+					vn = sqrtf(vn);
+					parts[i].vx *= vn;
+					parts[i].vy *= vn;
+				}
+				else if(!parts[r>>8].ctype && parts[i].ctype & 0xFFFFFFC0)
+				{
+					parts[r>>8].life = 1;
+					parts[r>>8].ctype = parts[i].ctype;
+					parts[r>>8].tmp = (0xFFFF & (int)(parts[i].vx * 255.0f)) | (0xFFFF0000 & (int)(parts[i].vy * 16711680.0f));
+					parts[r>>8].tmp2 = (0xFFFF & (int)((parts[i].x - x) * 255.0f)) | (0xFFFF0000 & (int)((parts[i].y - y) * 16711680.0f));
+					kill_part(i);
+				}
+			}
 			else if ((r&0xFF) == PT_INVIS)
 			{
 				if (pv[ny/CELL][nx/CELL]<=4.0f && pv[ny/CELL][nx/CELL]>=-4.0f)
@@ -2292,7 +2335,7 @@ int Simulation::try_move(int i, int x, int y, int nx, int ny)
 		}
 		else if (parts[i].type == PT_NEUT)
 		{
-			if ((r&0xFF) == PT_GLAS)
+			if ((r&0xFF) == PT_GLAS || (r&0xFF) == PT_BGLA)
 				if (rand() < RAND_MAX/10)
 					create_cherenkov_photon(i);
 		}
@@ -2535,7 +2578,7 @@ int Simulation::is_blocking(int t, int x, int y)
 	if (t & REFRACT) {
 		if (x<0 || y<0 || x>=XRES || y>=YRES)
 			return 0;
-		if ((pmap[y][x] & 0xFF) == PT_GLAS)
+		if ((pmap[y][x] & 0xFF) == PT_GLAS || (pmap[y][x] & 0xFF) == PT_BGLA)
 			return 1;
 		return 0;
 	}
@@ -3344,7 +3387,7 @@ void Simulation::create_cherenkov_photon(int pp)//photons from NEUT going throug
 
 	nx = (int)(parts[pp].x + 0.5f);
 	ny = (int)(parts[pp].y + 0.5f);
-	if ((pmap[ny][nx] & 0xFF) != PT_GLAS)
+	if ((pmap[ny][nx] & 0xFF) != PT_GLAS && (pmap[ny][nx] & 0xFF) != PT_BGLA)
 		return;
 
 	if (hypotf(parts[pp].vx, parts[pp].vy) < 1.44f)
@@ -4229,7 +4272,9 @@ void Simulation::UpdateParticles(int start, int end)
 					{
 						int rt = pmap[fin_y][fin_x] & 0xFF;
 						int lt = pmap[y][x] & 0xFF;
-						if ((rt==PT_GLAS && lt!=PT_GLAS) || (rt!=PT_GLAS && lt==PT_GLAS))
+						int rt_glas = (rt == PT_GLAS) || (rt == PT_BGLA);
+						int lt_glas = (lt == PT_GLAS) || (lt == PT_BGLA);
+						if ((rt_glas && !lt_glas) || (lt_glas && !rt_glas))
 						{
 							if (!get_normal_interp(REFRACT|t, parts[i].x, parts[i].y, parts[i].vx, parts[i].vy, &nrx, &nry)) {
 								kill_part(i);
@@ -4242,11 +4287,11 @@ void Simulation::UpdateParticles(int start, int end)
 								kill_part(i);
 								continue;
 							}
-							nn = GLASS_IOR - GLASS_DISP*(r-15)/15.0f;
+							nn = GLASS_IOR - GLASS_DISP*(r-30)/30.0f;
 							nn *= nn;
 							nrx = -nrx;
 							nry = -nry;
-							if (rt==PT_GLAS && lt!=PT_GLAS)
+							if (rt_glas && !lt_glas)
 								nn = 1.0f/nn;
 							ct1 = parts[i].vx*nrx + parts[i].vy*nry;
 							ct2 = 1.0f - (nn*nn)*(1.0f-(ct1*ct1));
@@ -4319,17 +4364,35 @@ void Simulation::UpdateParticles(int start, int end)
 					else if ((r & 0xFF) == PT_URAN) parts[i].ctype &= 0x003FC000;
 					else if ((r & 0xFF) == PT_GOLD) parts[i].ctype &= 0x3C038100;
 
-					if (get_normal_interp(t, parts[i].x, parts[i].y, parts[i].vx, parts[i].vy, &nrx, &nry)) {
-						dp = nrx*parts[i].vx + nry*parts[i].vy;
-						parts[i].vx -= 2.0f*dp*nrx;
-						parts[i].vy -= 2.0f*dp*nry;
+					if (get_normal_interp(t, parts[i].x, parts[i].y, parts[i].vx, parts[i].vy, &nrx, &nry))
+					{
+						if ((r & 0xFF) == PT_CRMC)
+						{
+							float r = (rand() % 101 - 50) * 0.01f, rx, ry, anrx, anry;
+							r = r * r * r;
+							rx = cosf(r); ry = sinf(r);
+							anrx = rx * nrx + ry * nry;
+							anry = rx * nry - ry * nrx;
+							dp = anrx*parts[i].vx + anry*parts[i].vy;
+							parts[i].vx -= 2.0f*dp*anrx;
+							parts[i].vy -= 2.0f*dp*anry;
+						}
+						else
+						{
+							dp = nrx*parts[i].vx + nry*parts[i].vy;
+							parts[i].vx -= 2.0f*dp*nrx;
+							parts[i].vy -= 2.0f*dp*nry;
+						}
 						// leave the actual movement until next frame so that reflection of fast particles and refraction happen correctly
-					} else {
+					}
+					else
+					{
 						if (t!=PT_NEUT)
 							kill_part(i);
 						continue;
 					}
-					if (!(parts[i].ctype&0x3FFFFFFF) && t == PT_PHOT) {
+					if (!(parts[i].ctype&0x3FFFFFFF) && t == PT_PHOT)
+					{
 						kill_part(i);
 						continue;
 					}
diff --git a/src/simulation/elements/C5.cpp b/src/simulation/elements/C5.cpp
index 18cbfa345c..bff28e603e 100644
--- a/src/simulation/elements/C5.cpp
+++ b/src/simulation/elements/C5.cpp
@@ -30,7 +30,7 @@ Element_C5::Element_C5()
 	HeatConduct = 88;
 	Description = "Cold explosive, set off by anything cold.";
 
-	Properties = TYPE_SOLID | PROP_NEUTPENETRATE;
+	Properties = TYPE_SOLID | PROP_NEUTPENETRATE | PROP_LIFE_DEC;
 
 	LowPressure = IPL;
 	LowPressureTransition = NT;
@@ -42,6 +42,7 @@ Element_C5::Element_C5()
 	HighTemperatureTransition = NT;
 
 	Update = &Element_C5::update;
+	Graphics = &Element_C5::graphics;
 }
 
 //#TPT-Directive ElementHeader Element_C5 static int update(UPDATE_FUNC_ARGS)
@@ -66,8 +67,57 @@ int Element_C5::update(UPDATE_FUNC_ARGS)
 					}
 				}
 			}
+	if (parts[i].ctype && !parts[i].life)
+	{
+		float vx = ((parts[i].tmp << 16) >> 16) / 255.0f;
+		float vy = (parts[i].tmp >> 16) / 255.0f;
+		float dx = ((parts[i].tmp2 << 16) >> 16) / 255.0f;
+		float dy = (parts[i].tmp2 >> 16) / 255.0f;
+		r = sim->create_part(-3, x, y, PT_PHOT);
+		if (r != -1)
+		{
+			parts[r].ctype = parts[i].ctype;
+			parts[r].x += dx;
+			parts[r].y += dy;
+			parts[r].vx = vx;
+			parts[r].vy = vy;
+			parts[r].temp = parts[i].temp;
+		}
+		parts[i].ctype = 0;
+		parts[i].tmp = 0;
+		parts[i].tmp2 = 0;
+	}
 	return 0;
 }
 
+//#TPT-Directive ElementHeader Element_C5 static int graphics(GRAPHICS_FUNC_ARGS)
+int Element_C5::graphics(GRAPHICS_FUNC_ARGS)
+
+{
+	if(!cpart->ctype)
+		return 0;
+
+	int x = 0;
+	*colr = *colg = *colb = 0;
+	for (x=0; x<12; x++) {
+		*colr += (cpart->ctype >> (x+18)) & 1;
+		*colb += (cpart->ctype >>  x)     & 1;
+	}
+	for (x=0; x<12; x++)
+		*colg += (cpart->ctype >> (x+9))  & 1;
+	x = 624/(*colr+*colg+*colb+1);
+	*colr *= x;
+	*colg *= x;
+	*colb *= x;
+
+	*firea = 100;
+	*firer = *colr;
+	*fireg = *colg;
+	*fireb = *colb;
+
+	*pixel_mode &= ~PMODE_FLAT;
+	*pixel_mode |= FIRE_ADD | PMODE_ADD | NO_DECO;
+	return 0;
+}
 
 Element_C5::~Element_C5() {}
diff --git a/src/simulation/elements/PHOT.cpp b/src/simulation/elements/PHOT.cpp
index c2d5630b67..9ef33dd5f6 100644
--- a/src/simulation/elements/PHOT.cpp
+++ b/src/simulation/elements/PHOT.cpp
@@ -79,7 +79,7 @@ int Element_PHOT::update(UPDATE_FUNC_ARGS)
 						sim->pv[y/CELL][x/CELL] -= 15.0f * CFDS;
 					}
 				}
-				else if((r&0xFF) == PT_QRTZ && !ry && !rx)//if on QRTZ
+				else if(((r&0xFF) == PT_QRTZ || (r&0xFF) == PT_PQRT) && !ry && !rx)//if on QRTZ
 				{
 					float a = (rand()%360)*3.14159f/180.0f;
 					parts[i].vx = 3.0f*cosf(a);
@@ -88,6 +88,15 @@ int Element_PHOT::update(UPDATE_FUNC_ARGS)
 						parts[i].ctype = 0x1F<<(rand()%26);
 					parts[i].life++; //Delay death
 				}
+				else if((r&0xFF) == PT_BGLA && !ry && !rx)//if on BGLA
+				{
+					float a = (rand()%101 - 50) * 0.001f;
+					float rx = cosf(a), ry = sinf(a), vx, vy;
+					vx = rx * parts[i].vx + ry * parts[i].vy;
+					vy = rx * parts[i].vy - ry * parts[i].vx;
+					parts[i].vx = vx;
+					parts[i].vy = vy;
+				}
 				else if ((r&0xFF) == PT_FILT && parts[r>>8].tmp==9)
 				{
 					parts[i].vx += ((float)(rand()%1000-500))/1000.0f;