(*************** Cash Flow Calculations During Construction *****************)
usdrsconversion = 55;
overnightcost=4000; (* in dollars per megawatt *)
capacity = 1650 *10^3; (*In Kilowatts *)
inflation = .05;
initialcostestimate = usdrsconversion*overnightcost*capacity;
startyear = 0;(*tells us, which year is the $0$ year *)
(*Find debt and equity for each year *)
cutoffyear = 6; (* This is put by hand. It is an estimate of when the company stops putting in its own equity. You may need to change this number, depending on the S-curve and the inflation rate if you demand that the final ratios should be 70:30*)
debtinterestrate = .1321; (*This corresponds to a nominal interest rate. To get the "real interest" rate subtract inflation. *)
eqrate = .1683;
(*This is again a nominal rate. The real rate is given by substracing the inflation rate *)
(*d = (1 + 0.7*0.7*debtinterestrate + 0.3*eqrate) ; *)(*discount rate used to discount future payments and energy production. This is a weighted average of the debt interest rate and the return on equity. The debt interest rate is weighted by 0.7*0.7. The first 0.7 comes from the fact that the debt is 0.7 of the total capital. The other 0.7 comes from the fact that the corporate tax rate is 30%, and you don't have to pay tax on the debt interest payments. *)
eqratenpcil = 0.14; (*This is the rate used by the NPCIL.*)

(******************* Expenditure Curve using data from Kudankulam ************)

expenditurecurve ={2454,4168, 5476, 7634, 9571, 10592, 11947, 13002, 13811, 14811, 15824}
(* This data comes from the MOSPI reports on the progressive expenditure on Kudankulam *)
years = Length[expenditurecurve];
yearlyexpenditure = expenditurecurve - PadLeft[expenditurecurve[[1;;years-1]], years];
escalation = Table[(1 + inflation)^(-n), {n, 0, years-1}];
inflationadjustyearlyexpend = yearlyexpenditure*escalation;
dscurve = inflationadjustyearlyexpend/Total[inflationadjustyearlyexpend]
scurve = Table[Total[dscurve[[1;;n]]], {n, 1, years}]
(* S curve is the cumulative amount spent each year. ds curve is the amount spent each year. *)

base = initialcostestimate*dscurve;
escalatedbase = base*Table[(1 + inflation)^n, {n, startyear, years+startyear-1}];
intondebt = totalexpenditure = equity = debt= Table[0, {i, 1, years}];

equity[[1]] = escalatedbase[[1]];
totalexpenditure[[1]] = escalatedbase[[1]];
For[i = 2, i < cutoffyear, i++, intsofar = Sum[debt[[j]], {j, 1, i-1}]*debtinterestrate;debt[[i]] = (1/2) (escalatedbase[[i]] +intsofar)/(1 - debtinterestrate/4); 
intondebt[[i]] = intsofar + (debtinterestrate/2)*debt[[i]];
totalexpenditure[[i]] = escalatedbase[[i]] + intsofar + (debtinterestrate/2)*debt[[i]];
equity[[i]] = debt[[i]]]
(***** Cutoff Year ***********************)
equity[[cutoffyear]] = x; 
intsofar =  Sum[debt[[j]], {j, 1, cutoffyear-1}]*debtinterestrate;
debt[[cutoffyear]] = (escalatedbase[[cutoffyear]] + intsofar  - equity[[cutoffyear]])/(1 - debtinterestrate/2);
intondebt[[cutoffyear]] = intsofar + (debtinterestrate/2)*debt[[cutoffyear]];
totalexpenditure[[cutoffyear]] = equity[[cutoffyear]] + debt[[cutoffyear]];
(***** Beyond the Cutoff Year **********)
(**** Now there is no equity. There is no business of a penultimate year, when one reactor starts functioning *********)

For[i = cutoffyear+1, i ≤  years, i++, intsofar = Sum[debt[[j]], {j, 1, i-1}]*debtinterestrate;debt[[i]] =  (escalatedbase[[i]] +intsofar)/(1 - debtinterestrate/2); 
intondebt[[i]] = intsofar + (debtinterestrate/2)*debt[[i]];totalexpenditure[[i]] = escalatedbase[[i]] + intsofar + (debtinterestrate/2)*debt[[i]];
equity[[i]] = 0]

totaldebt = Total[debt];
totalexp = Total[totalexpenditure];
sol = Solve[totaldebt == 0.7 totalexp, x][[1]];
cashflowtable = {base,  escalatedbase , debt , intondebt , equity, totalexpenditure} /. sol;

TableForm[cashflowtable/10^7]

roundedcashflowtable = Round[cashflowtable/10^7];

TableForm[roundedcashflowtable]

totalcost = Total[totalexpenditure /. sol];
totalequity = Total[equity /. sol];
totaldebt = Total[debt /. sol];

Print["Total Cost is ", totalcost/10^7, " crores"]
Print["Ratio of Total cost to initial cost estimate is ", totalcost/initialcostestimate]
Print["Total Equity is ", totalequity/10^7, " crores"]
Print["Total Debt is ", totaldebt/10^7, " crores"]

(******************************** Effective Rate on Equity **************************************************)

effectiverate[naiverate_?NumericQ] := Module[{r}, soleffectrate = FindRoot[naiverate*totalequity/Sum[(equity[[i]] (1+r)^(years - i) /. sol), {i, 1, years}] == r, {r,  naiverate}];
 r/. soleffectrate]
(*This function solves for the effective rate of return on equity *)

(*The equity spent in initial years is sitting idle. So we can self consistently solve for the effective return on equity*)
Print["The effective rate on equity is ", effectiverate[eqrate]]
Print["The effective rate on equity with the NPCIL's assumptions is ", effectiverate[eqratenpcil]]

FindRoot[effectiverate[t] == debtinterestrate,  {t, 0.25} ]

{{0.1, 1.0}, {0.06, 1.0}}

effectiverateplot = Plot[effectiverate[x], {x, 0.1, 0.372},PlotRange -> Automatic]

Export["c://users//suvrat//home//politics//epreconomics//effectiverateplot.png",effectiverateplot]

xmin = 0.1;
xmax = 0.372;
maxc = 100;
effectiveratelist = Table[{xmin+ c*(xmax - xmin)/maxc, effectiverate[xmin + c*(xmax - xmin)/maxc]}, {c, 1, maxc}];
(* For exporting into MS word! *)

(******************************* The cost of fuel *************************************)
(**** Fuel Usage **********)
burnup = 50; (* in MWth-day/kg*)
electout = 1650; (* gross electrical output in MW *)
thermalout = 4500; (*thermal output in MW*)
efficiency = electout/thermalout;
loadfactor = 0.8;

(*check Areva's figure for total lifetime energy production with a higher load factor *)
4500*60*365*.93

fuelkgperunit= 1000/(24*(burnup*efficiency)) (* kg per million units i.e. per 1000 MW-hr *)
kgperyear = (thermalout/burnup)*365*loadfactor

(****** production cost for LEU *******
Note that conservation of mass tells us that for a total feed, F, tail T, and product P, we must have:
F = T + P; x_f F = x_t T + x_p P;
and so F/P = (x_p - x_t)/(x_f - x_t)*)
xp = 5*(1/100);
xt = 0.2 * (1/100);
xf = 0.7*(1/100);
F = (xp - xt)/(xf - xt)
(*amount of raw uranium required to produce one unit of leu *)
V[x_] := (2 x - 1) Log[x/(1 - x)]
(* separative work per unit is  *)
Ns = V[xp] + (F - 1) V[xt] - F V[xf]

Cu = 150;
Cc = 10;
Cs = 160;
Cf = 250;
fuelrate = (F (Cu + Cc) + Ns Cs + Cf)*1.1

(******************** Calculation of Tariff **********************)


pmt[P_, r_,nyears_] := r (1 + r)^(nyears) P/((1 + r)^(nyears) - 1)
(*Payment function defined for later usage. The idea here is as follows. You want to make a constant payment each year, so that at the end of n years, you have recovered the initial amount and the compound interest it would have earned at the rate r. However, when you make a payment in the first year, you can put that in the bank and earn compound interest at the same rate for (n-1) years ... the second year payment is earning for (n-2) years ... etc. This leads to the formula 
c (charge) = P (1 + r)^n/Sum[(1 + r)^j, {j, 0, n-1}].
An alternate way to understand this is that when you make a payment for the first year, the amount you need to repay reduces, and so now you need to pay interest only on that amount. *)

(************* Calculation of Year by Year Tariff *******************)
lifetimeofreactor = 40;
yearlydetailedtariffs = Range[lifetimeofreactor];
tariffs = Range[lifetimeofreactor];
totalunitsperyear = capacity* 24*365*loadfactor;
capacityfactor = .93;
omrate = .02;
depreciationrate = .0528;
depreciationyears = Round[0.9/depreciationrate];
(*In this revenue model, depreciation happens at a constant nominal rate. So, after 25 years, 90% of the value is gone, and then we assume that depreciation stops.*)
netunits = totalunitsperyear*capacityfactor;
repaymentyears = 8; (*This is the number of years over which you repay the full debt *)
(*Return on Equity *)
recoveryfactor = 1; 
fuelrecoveryint = debtinterestrate;
(*This is the amount of the initial fuel cost that is recovered. Presumably you want all of it to be recovered, and not just 2/3, which is what Thakur assumes. Thakur also assumes an interest rate of 11% for this, but there is no reason for this to be different from the rate of interest on debt. *)
fuelrecoveryyears = 15;
initialfuelinventory = 127*10^3; (*The NPCIL figure is 84 * 10^3 *)
(*fuelkgperunit = 3.15; Thakur's usage per million units. See above for the EPR usage *)
initialfuelcost = initialfuelinventory*fuelrate*usdrsconversion*(1+inflation)^(years + startyear - 1);
intonworkingcap = .12;
earlyfuelinflationrate = .05;
latefuelinflationrate = .05;
variableinflationrate = .05;
(** Parameters for Calculating Working Capital ***)
wctc = .02; (*fraction of total cost *)
wcom = 1/12; (*fraction of O&M costs *)
wcfuel = 1/2; (*fraction of fuel cost*)
wctariff = 1/6; (*fraction of total tariff collection for the year *)
dccharge = .02;(* Set by hand *)

debug = Range[lifetimeofreactor];

calctariffs[multiplier_] := Module[{tariffs = Range[lifetimeofreactor],yearlydetailedtariffs = Range[lifetimeofreactor], totalcostmult, totalequitymult, totaldebtmult},totalcostmult = totalcost*multiplier; totalequitymult = totalequity*multiplier; totaldebtmult = totaldebt*multiplier;
(*This function calculates tariffs as a function of a "multiplier" which allows one to scale the overnight cost*)
For[n  =  0, n < lifetimeofreactor, n++,
cost = {};
variablescalation = (1 + variableinflationrate)^n;
(****** Return on Equity *****)
eq = totalequitymult;
returnonequity = eq*eqrate/netunits;
cost = Append[cost, {"returnonequity", returnonequity}];
(*Interest on Market Borrowings ********)
intmarketborr  = totaldebtmult*((2 repaymentyears - 1 - 2 n)/(2 repaymentyears))*debtinterestrate/netunits;
If[ n ≥  repaymentyears,  intmarketborr = 0];cost = Append[cost, {"intonmarketborrowings", intmarketborr}];
(****************** depreciation *****************)
depreciation = (depreciationrate*totalcostmult)/netunits;
If[n ≥ depreciationyears, depreciation = 0];
cost = Append[cost, {"depreciation", depreciation}];
(*** OM Costs **********)
omcosts = variablescalation*omrate*totalcostmult;
omcoststariff = omcosts/netunits;
cost = Append[cost, {"omcosts tariff", omcoststariff}];
(*fuel costs *)
totalinflation = (1 + earlyfuelinflationrate)^(years+startyear-1) (1 + latefuelinflationrate)^n;
fuelcost = (totalunitsperyear/10^6 )*fuelkgperunit*fuelrate*usdrsconversion*totalinflation;
fueltariff = fuelcost/netunits;
cost = Append[cost, {"fueltariff", fueltariff}];
(****** fuel recovery costs ******)

fuelrecoverycost = (pmt[ (recoveryfactor)*initialfuelcost, fuelrecoveryint, fuelrecoveryyears]/netunits);
If[n ≥ fuelrecoveryyears, fuelrecoverycost = 0];cost = Append[cost, {"fuelrecoverycost", fuelrecoverycost}];
(** These costs stop after the end of the fuel recovery years ***)
(*****************working capital ***************)
(*The logic is that the working capital is 2% of completion cost + 1 month of OM costs (which are also 2% of completion costs initially, but then escalate because of inflation) + 6 months of fuel + 2 months of energy receivables. The last figure depends on the tariff $t$. In the formula below, recall that omcosts and fuelcosts are already escalated. So we explicitly escalate only the part that depends on the "total cost"*)
workingcapital = variablescalation*wctc*totalcostmult + wcom*omcosts+ wcfuel*fuelcost + t*netunits*wctariff;

intworkingcapital = workingcapital*intonworkingcap/netunits;
cost = Append[cost, {"intworkingcapital", intworkingcapital}];
(********** Decommissioning ****************)
cost = Append[cost, {"decommissioning", dccharge}];

tariffsol = Solve[t == Total[Transpose[cost][[2,All]]], t][[1]];
yearlydetailedtariffs[[n+1]] = (cost/.tariffsol);
debug[[n+1]] = {variablescalation*wctc*totalcost , wcom*omcosts, wcfuel*fuelcost,  t*netunits*wctariff} /. tariffsol;
tariffs[[n+1]] = Total[Transpose[yearlydetailedtariffs[[n+1]]][[2, All]]]]; {yearlydetailedtariffs, tariffs}]

yearlydetailedtariffs2 = Range[lifetimeofreactor];
tariffs2 = Range[lifetimeofreactor];
yearlydetailedtariffs3 = Range[lifetimeofreactor];
tariffs3 = Range[lifetimeofreactor];
yearlydetailedtariffs4 = Range[lifetimeofreactor];
tariffs4 = Range[lifetimeofreactor];
yearlydetailedtariffs5 = Range[lifetimeofreactor];
tariffs5 = Range[lifetimeofreactor];
yearlydetailedtariffsanalytic = Range[lifetimeofreactor];
tariffsanalytic = Range[lifetimeofreactor];

{yearlydetailedtariffs, tariffs} =calctariffs[1.0];
{yearlydetailedtariffs2, tariffs2} =calctariffs[0.75];
{yearlydetailedtariffs3, tariffs3} =calctariffs[0.5];
{yearlydetailedtariffs4, tariffs4} =calctariffs[1.25];
{yearlydetailedtariffs5, tariffs5} =calctariffs[1.5];
{yearlydetailedtariffsanalytic, tariffsanalytic} =calctariffs[ζ];

tariffs

tariffplot = ListPlot[tariffs, PlotRange -> {{0,40},{6,18}}, Joined -> True, PlotMarkers -> Automatic]
tariffplotsampleovernight = ListPlot[{tariffs, tariffs2, tariffs3, tariffs4}, PlotRange -> {{0,40},{4,18}}, Joined -> True, PlotMarkers -> Automatic]

Clear[tariffslist];
Clear[data];
nopoints =9;
data = Range[nopoints*lifetimeofreactor];
count = 1;
firstyeartariffs=Range[nopoints];
lowmult = 0.5;
highmult = 1.5;
For[m = 1, m ≤ nopoints, m++, multiplier =lowmult + (highmult - lowmult)*(m-1)/(nopoints-1);tariffslist[m]= calctariffs[multiplier][[2]]; firstyeartariffs[[m]] = {multiplier*overnightcost,tariffslist[m][[1]]};For[m2 = 1, m2 ≤ Length[tariffslist[m]], m2++, data[[count++]] = {m2, multiplier*overnightcost,  tariffslist[m][[m2]]}]]

alldetailsplot = ListPlot3D[data, ColorFunction -> "Rainbow"]

firstyearvsovernightplot = ListPlot[firstyeartariffs]

(************** The tariff varies linearly with the overnight cost of construction. The array tariffanalytic contains data about these coefficients. The formula is: t_n = f_n + v_n (O/4000) *)
tariffsanalytic = Simplify[tariffsanalytic]
fn = Range[lifetimeofreactor];
vn = Range[lifetimeofreactor];
fn = tariffsanalytic /. ζ -> 0
vn = Chop[Simplify[(tariffsanalytic -  fn )/ζ]]

Export["c://users//suvrat//home//politics//epreconomics//tariffplot.png", tariffplot]
Export["c://users//suvrat//home//politics//epreconomics//tariffplotsampleovernight.png", tariffplotsampleovernight]
Export["c://users//suvrat//home//politics//epreconomics//firstyearvsovernight.png", firstyearvsovernightplot]
Export["c://users//suvrat//home//politics//epreconomics//alldetails.png", alldetailsplot]

(************ First Year Tariff **********)
TableForm[yearlydetailedtariffs[[1]]]

(*********** Levelized Tariff Calculation *************)
d = 0.7*0.7*debtinterestrate + 0.3*eqrate; (*discount rate used to discount future payments and energy production obtained as a weighted sum of capital. The debt is weighted with 0.7*0.7 because it is 70% of the total, and also because the corporate tax rate is 30% and you don't have to pay tax on the interest on debt. *)
Sum[tariffs[[n]] /(1 + d)^(n-1), {n, 1, lifetimeofreactor}]/Sum[1/(1 + d)^(n-1),  {n, 1, lifetimeofreactor}]
(** The significance of this calculation is a little unclear since it discounts electricity production and nominal tariff at the same discount rate. It may make more sense to discount electricity at a "real discount rate" and tariffs at a "nominal discount rate." **)

(********************** Print Out Details of Each Years Tariff ************)
For[n = 1, n ≤ lifetimeofreactor, n++, Print["********* Year ", n-1, "**************"]; Print[TableForm[yearlydetailedtariffs[[n]]]]; Print["===================="];
Print["Total Tariff is ", tariffs[[n]]];
Print["======================"];Print["*** Working Capital is ", yearlydetailedtariffs[[n,7,2]]*netunits/intonworkingcap]]

(******** Export this information to Excel *************)
For[n = 1, n ≤ lifetimeofreactor, n++, debug[[n]] = Join[yearlydetailedtariffs[[n,1;;8,2]], {tariffs[[n]]}]]; 
Export["c://users//suvrat//home//politics//epreconomics//detailedtariff.xls", Transpose[debug], "XLS"]

For[n = 1, n ≤ lifetimeofreactor, n++, debug[[n]] = {n, tariffs[[n]], tariffs2[[n]], tariffs3[[n]], tariffs4[[n]], tariffs5[[n]]}];
Export["c://users//suvrat//home//politics//epreconomics//tariffcompare.xls", debug, "XLS"]

(*debug = Range[nopoints];
For[n = 1, n ≤ nopoints, n++, debug[[n]] = {(lowmult + (highmult - lowmult)*m/nopoints)*overnightcost, firstyeartariffs[[n]]}];*)
(* PUT Correct formula here*)
Export["c://users//suvrat//home//politics//epreconomics//firstyeartariffs.xls", firstyeartariffs, "XLS"] 
Export["c://users//suvrat//home//politics//epreconomics//tariffs3d.xls", data, "XLS"] 
Export["c://users//suvrat//home//politics//epreconomics//tariffcoefficients.xls", Transpose[{vn,fn}], "XLS"]

Transpose[{vn,fn}]

debug[[1]]

firstyeartariffs[[1]]