code: public static RandomVariableInterface bachelierOptionValue(
  			RandomVariableInterface forward,
  			RandomVariableInterface volatility,
  			double optionMaturity,
  			double optionStrike,
  			RandomVariableInterface payoffUnit)
  	{
  		if(optionMaturity < 0) {
  			return forward.mult(0.0);
  		}
  		else
  		{
  			RandomVariableInterface integratedVolatility = volatility.mult(Math.sqrt(optionMaturity));
  			RandomVariableInterface dPlus	= forward.sub(optionStrike).div(integratedVolatility);
  
  			UnivariateFunction cummulativeNormal = new UnivariateFunction() {
  				public double value(double x) { return NormalDistribution.cumulativeDistribution(x); }
  			};
  			UnivariateFunction densityNormal	= new UnivariateFunction() {
  				public double value(double x) { return NormalDistribution.density(x); }
  			};
  			RandomVariableInterface valueAnalytic = dPlus.apply(cummulativeNormal).mult(forward.sub(optionStrike))
  					.add(dPlus.apply(densityNormal).mult(integratedVolatility)).mult(payoffUnit);
  
  			return valueAnalytic;
  		}
  	}
docs: Calculates the option value of a call, i.e., the payoff max(S(T)-K,0) P, where S follows a
  normal process with constant volatility, i.e., a Bachelier model
  \[
  \mathrm{d} S(t) = r S(t) \mathrm{d} t + \sigma \mathrm{d}W(t)
  \]
  
  @param forward The forward of the underlying \( F = S(T) \exp(r T) \).
  @param volatility The Bachelier volatility \( \sigma \).
  @param optionMaturity The option maturity T.
  @param optionStrike The option strike.
  @param payoffUnit The payoff unit (e.g., the discount factor)
  @return Returns the value of a European call option under the Bachelier model.