{$R+}  {Range check}

(* CDISUNIF.PAS

Calculate the uniform contribution distribution coefficients that
give the target internal rate of return for a given lifespan *)

program CDistUnifProgram;

uses SMUCTPU;

const progDesc = 'Calculate uniform contribution distribution for target IRR';
      progName = 'CDISUNIF';
      progDate = 'Thu 5-Jun-97';
      progVers = '1.0';
      maxYears = 50;

type realVectorType = array [0..maxYears] of real;

var lifeSpan
     : integer;

var targetIrr
     : real;

var b               { the sought contribution coefficients }
     : realVectorType;

(* Get the parameters *)
procedure GETPARAM (var lifeSpan  : integer;
                    var targetIRR : real);
var s : string;
    k : integer;
    target : real;
begin
  if ParamCount < 2 then begin
    writeln ('Usage: CDISUNIF lifeSpan targetIRR%');
    halt;
  end;
  s := ParamStr(1);
  Val (s, lifeSpan, k);
  if k > 0 then begin writeln ('Error in value ', s); halt; end;
  s := ParamStr(2);
  Val (s, target, k);
  if k > 0 then begin writeln ('Error in value ', s); halt; end;
  targetIRR := target / 100.0;
end; (* getparam *)

(* Auxiliary function to calculate net present value *)
function NPVFN (b        : realVectorType;
                rate     : real;
                lifeSpan : integer) : real;
var y : real;
    i : integer;
begin
  y := -1.0;
  for i := 0 to lifeSpan do
    y := y + b[i] / GENPOWFN (1.0 + rate, i);
  npvfn := y;
end;  (* npvfn *)

(* Internal rate of return for the contribution distribution *)
function IRRFN (b : realVectorType; lifeSpan : integer) : real;
const maxIter = 30;
var x, x1, x2, y1, y2 : real;
    j : integer;
begin
  x1 := 0.0;
  x2 := 0.2;
  for j := 1 to maxIter do begin
    y1 := NPVFN (b, x1, lifeSpan);
    y2 := NPVFN (b, x2, lifeSpan);
    x := (x1*y2 - x2*y1) / (y2 - y1);
    x1 := x2;
    x2 := x;
    irrfn := x;
    if abs(x2-x1) < 0.000001 then exit;
  end; {for}
  writeln ('The IRRFN algorithm failed in ', maxIter, ' iterations');
  halt;
end;  (* irrfn *)

(* Calculate the contribution coeffiefficient, uniform distribution *)
procedure CALCDIST (n     : integer;
                    r     : real;
                    var b : realVectorType);
var i, j : integer;
    sum  : real;
begin
  for i := 1 to n do begin
    sum := 0;
    for j := 1 to n do begin
      sum := sum + GENPOWFN (1.0 + r, -j)
    end;
    b[i] := 1.0 / sum;
  end;
end; (* calcdist *)


(* Present the results *)
procedure SHOW (lifeSpan : integer;
                b        : realVectorType);
var irr : real;
    i   : integer;
begin
  writeln;
  irr := IRRFN (b, lifeSpan);
  writeln ('IRR = ', 100.0 * irr :0:8, '%');
  writeln;
  for i := 1 to lifespan do writeln (b[i]:10:8, '  #b(', i, ')');
end;  (* show *)

(* Main program *)
begin
  HEADER (progName, progDesc, progVers, progDate, '');
  writeln;
  GETPARAM (lifeSpan, targetIRR);
  CALCDIST (lifespan, targetIRR, b);
  SHOW (lifeSpan, b);
end. (* CDisUnifProgram *)

