program assignment4_1b;

uses Graph, crt, math, sysutils;

var sigma, I0, maxx, maxy:real;
    j, grDriver, grMode, ErrCode, decx, decy:integer;
    n:array[0..10] of integer = (2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96);
    I, x, y:array[0..10] of real;



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

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

end;



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

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

begin

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

  for j:=1 to (n-1) do begin

    if odd(j) then	yodd := yodd+func(j*(upperlim-lowerlim)/n+lowerlim)
    else		yeven:=yeven+func(j*(upperlim-lowerlim)/n+lowerlim);

  end;

  I:=(I+2*yeven+4*yodd)*(upperlim-lowerlim)/(3*n);

end;



procedure length(n:real; var digits:integer);

begin

  digits:=0;

  while n>=1 do begin

    n:=n/10;
    digits:=digits+1;

  end;

end;



procedure scale(n:real; var out:real);

var digits, i, mul:integer;

begin

  if (n<0) then mul:=-1
  else mul:=1;

  n:=n*mul;

  length(n, digits);

  while n>=100 do n:=n/10;

  n:=n/5;
  n:=ceil(n)*5;

  for i:=3 to digits do n:=n*10;

  out:=n*mul;

end;



procedure putonx(n:real; x, y, dec:integer);

var digits:integer;
    output:ShortString;

begin

  length(n,digits);

  SetLength(output,FloatToText(@output[1], n, ffFixed, digits, dec));

  MoveTo((x-ceil(TextWidth(output)/2)), y);
  OutText(output);

end;



procedure putony(n:real; x, y, dec:integer);

var digits:integer;
    output:ShortString;

begin

  length(n,digits);

  SetLength(output,FloatToText(@output[1], n, ffFixed, digits, dec));

  x:=x-ceil(TextWidth(output));
  y:=y-ceil(TextHeight(output)/2);

  if x<0 then x:=0;
  if y<0 then y:=0;

  MoveTo(x, y);
  OutText(output);

end;



procedure drawcross(x, y, maxx, maxy:real; size:integer);

var xpixel, ypixel:integer;

begin

  xpixel:=round(x*1220/maxx)+59;
  ypixel:=round(y*900/maxy)+59;

  Line(xpixel-size, ypixel-size, xpixel+size, ypixel+size);
  Line(xpixel+size, ypixel-size, xpixel-size, ypixel+size);

end;



begin

  I0:=0.68268949213708589717;

  writeln('What should the standard deviation (sigma) be?');
  readln(sigma);

  writeln('  n |      I     |      x     |       y      ');
  writeln('----+------------+------------+--------------');

  for j:=0 to 10 do begin
    integrate(@gauss, -sigma, sigma, n[j], I[j]);

    x[j]:=ln(n[j]);
    y[j]:=ln(I[j]-I0+1E-6);

    writeln(' ',n[j]:2,' | ',I[j]:10:8,' | ',x[j]:10:8,' | ',y[j]:12:8);

  end;




  grDriver:=Detect;
  InitGraph(grDriver, grMode, '');
  ErrCode:=GraphResult;

  if ErrCode <> grOk then begin

    writeln('Graphics error:', GraphErrorMsg(ErrCode));
    readln; halt;

  end else begin

    Line(59, 0, 59, 959);
    Line(0, 59, 1279, 59);

    {Calculate the maximum values for the axes.}

    scale(x[10], maxx);
    scale(y[10], maxy);

    {Add axes tick marks and scales.}

    MoveTo(646, 15);
    OutText('log(n)');

    MoveTo(10, 464);
    OutText('y');

    if maxx/10=round(maxx/10) then decx:=0 else decx:=1;
    if maxy/10=round(maxy/10) then decy:=0 else decy:=1;

    for j:=1 to 10 do begin
      putonx(j*maxx/10, (j*122+59), 34, decx);
      putony(j*maxy/10, 44, (j*90+54), decy);

      Line((j*122+59), 49, (j*122+59), 59);
      Line(49, (j*90+54), 59, (j*90+54));
    end;

    for j:=0 to 10 do begin
      drawcross(x[j], y[j], maxx, maxy, 3);
    end;

  end;




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

end.