Easily Implement Stochastic Differential Equations, Bayesian Methods, and More

/*
 *	Simple Bethe-Bloch radiation transport implementation, exploiting UxHw.
 *	You can copy and paste directly into Signaloid Cloud Developer Platform
 *	code editor.
 */
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <uxhw.h>

const float		kInitialEnergy					= 800.0;
const float		kProtonMass						= 938.3;
const float		kElectronMass					= 0.511;
const float		kZconcrete						= 11;
const float		kIconcrete						= 0.000286;
const float		kLconcrete						= 0.1;
const int		kPartitionSize					= 100;
const float		kCollisionProbabilityDensityPerL= 50;
const float		kLowerBoundOnInelasticEnergyLoss= 0.7;
const float		kFudgeFactor1					= 0.15*kZconcrete;
const float		kFudgeFactor2					= 2*kElectronMass/kIconcrete;

float
betheBloch(float energy)
{
	float	beta = pow (1 - kProtonMass /(energy + kProtonMass), 2);
	float	dEdx = (kFudgeFactor1/beta)*log(kFudgeFactor2*beta);

	return dEdx;
}

int
main(void)
{
	/*
	 *	When running on UxHw-enabled platforms, all `float` and `double` values
	 *	have an associated distribution, in addition to their default semantics.
	 */
	float	energy = kInitialEnergy;
	float	dx = kLconcrete / kPartitionSize;
	float	inelasticCollisionProbability = dx*kCollisionProbabilityDensityPerL;

	for (int i=0; i<kPartitionSize; i++)
	{
		float	dEdx = betheBloch(energy);
		float	dE = fmaxf(dEdx*dx, 0);

		energy -= dE;
		energy = UxHwDoubleMixture(	energy,
						energy*UxHwFloatUniformDist(kLowerBoundOnInelasticEnergyLoss, 1),
						1 - inelasticCollisionProbability);
	}

	/*
	 *	When running on UxHw-enabled platforms, the variable `energy` at this point has
	 *	both its default floating-point semantics as well as an associated distribution.
	 *	You can use UxHw API functions to query any of the properties of this associated
	 *	distribution, or even display the full distribution.
	 */
	printf("Energy Absorbed = %f (mean value of distributional information is %f)\n",
			kInitialEnergy - energy,
            UxHwFloatNthMoment(energy,1));

	return 0