program assignment4_1c_commented;

var I, sigma, converge:real;



function gauss(x:real):real; far;
begin

//  This is the formula for the normal, or Gaussian, distribution, assuming that the mean is zero.

  gauss:=exp(-sqr(x/sigma)/2)/(sigma*sqrt(2*pi()));

end;



type funcparam = function(x:real):real;

procedure integrate(func:funcparam; lowerlim, upperlim, converge:real; var I:real);
var j, n:integer;
    yodd, yeven, Iprev, Icurr, h:real;

begin

{
  This function can integrate a general function for 1 variable using Simpson's rule:

  I = h(y0 + 4*y1 + 2*y2 + 4*y3 + 2*y4 + ... + 2*y(n-2) + 4*y(n-1) + yn)/3

  Give it a function, some limits, a precision to converge to and it will give you back the integral, I.

  This function attains the required precision by repeatedly doubling the number of strips (and hence y values).
}

  n:=2;	// Start with the smallest number of strips.

  yodd:=0;
  yeven:=0;
  I:=func(lowerlim)+func(upperlim);

  yodd:=func((upperlim-lowerlim)/2+lowerlim);

  Icurr:=(I+4*yodd)*(upperlim-lowerlim)/6;	// First stab at integration, but we still need something to compare it to.

  // And that's what this loop does, generates a second value for I using more strips. It will keep doing this until we have enough precision.
  repeat begin

    Iprev:=Icurr;	// Iprev is set to Icurr so that Icurr can be set to the new value at the end of the loop.

    h:=(upperlim-lowerlim)/(2*n);	// The width of each strip.

    yeven:=yeven+yodd;	// All of the odd strips from the previous pass will become even strips in this pass, there is no need to recalculate them.
    yodd:=0;		// Re-initialise yodd for the following loop.

    // This loop only calculates odd values, since all the odd values from previous passes are now the even values this time through.
    for j:=1 to n do yodd:=yodd+func((2*j-1)*h+lowerlim);

    Icurr:=(I+2*yeven+4*yodd)*h/3;	// We now have a new value for Icurr to compare to Iprev.

    n:=n*2;

  end;
  until (abs(Icurr-Iprev)<abs(converge));	// This is the comparison, if it fails the loop runs again, if it passes, we have our answer.

  I:=Icurr;

end;



begin

{
  writeln('What should the standard deviation (sigma) be?');
  readln(sigma);
  // Can uncomment this section if custom input is desired but there's ...
}

  sigma:=12.84375;	// ... something weird here by default to avoid masking errors from "nice" values (ie: 1).
  converge:=1E-6;	// When successive applications of Simpson's rule differ by less than this value, it stops.

  integrate(@gauss, -sigma, sigma, converge, I);	// Go!

  writeln('Between -sigma and sigma, I = ',I:12:10);			// Our output is printed nicely for the user.
  writeln('         erfc(-1/sqrt(2))-1 = 0.68268949213708589717');	// This is the correct value (to many decimal places) for reference.

  writeln();
  write('Press any key to continue.');
  readln;

end.