#!/usr/bin/env python
# (c) 2015 US EPA
# Originally Created - 20160812; post Draft TAR of 2016
#


#Set up
import pandas as pd

scriptversion = '20160924.1'
#py_version = pd.DataFrame(pd.Series(['Script Version']))
#py_version['Value'] = scriptversion

#Set the Inputs
inputFile = pd.ExcelFile('OMEGA_ICBT_Inputs_PD2016.xlsx')
outputFile = 'OMEGA_ICBT_Outputs_20161118_icm_aeoR.xlsx'

#Read Inputs
Basic = inputFile.parse('Basic', index_col=0)
General = inputFile.parse('General', index_col=0)
Fuel_Price = inputFile.parse('FuelPrices_Inputs')
Energy = inputFile.parse('Energy', index_col=0)
CO2Analysis = inputFile.parse('CO2_Analysis', index_col=0)
Sales = inputFile.parse('Sales_Analysis', index_col=0)
pop_in = inputFile.parse('population', index_col=0)
MAR_in = inputFile.parse('MAR', index_col=0)
survival = inputFile.parse('survival', index_col=0)
EconYear = inputFile.parse('EconYear')
VehicleEF = inputFile.parse('Vehicle_EFs', index_col=0)
ElecEF = inputFile.parse('Elec_EFs', index_col=0)
RefineryEF = inputFile.parse('Refinery_EFs', index_col=0)
AC = inputFile.parse('AC_Outputs', index_col=0)
Mass = inputFile.parse('Mass', index_col=0)
Insurance = inputFile.parse('Insurance')
print("Inputs read")

maxage = int(Basic.ix['maxage']['Value'])
startMY = int(Basic.ix['startMY']['Value'])
endMY = int(Basic.ix['endMY']['Value'])
startCY = int(Basic.ix['startCY']['Value'])
endCY = int(Basic.ix['endCY']['Value'])
vmt_GrowthRate = Basic.ix['VMT_GrowthRate']['Value']
basepopyear = int(Basic.ix['basepopyear']['Value'])
DR3 = General.ix['DR3']['Value']
DR7 = General.ix['DR7']['Value']
DR5 = General.ix['DR5']['Value']
DR2_5 = General.ix['DR2_5']['Value']
sDR3 = str(DR3)
sDR7 = str(DR7)
sDR5 = str(DR5)
sDR2_5 = str(DR2_5)
discountYr=General.ix['Discount values to year']['Value']
case_ref = str(General.ix['Reference case standards']['Value'])
case_control = str(General.ix['Control case standards']['Value'])
scenario = General.ix['Scenario for 2017-2025 standards']['Value']
baseline = General.ix['Baseline_fleetMY']['Value']
version = str(Basic.ix['Version']['Value'])
USD_Basis = str(General.ix['USD_Basis']['Value'])
taxes = float(General.ix['Sales taxes on increased vehicle costs']['Value'])
foreignOilFrac=float(General.ix['Foreign Oil Fraction']['Value'])
loanrate4=float(General.ix['LoanRate_4yr']['Value'])
loanrate5=float(General.ix['LoanRate_5yr']['Value'])
loanrate6=float(General.ix['LoanRate_6yr']['Value'])

million = 1000000
metric_to_uston = 1.10231131
print("Parameters set")

#>>>>>>>>>>>>>>>>>>> Create Calcs Matrix >>>>>>>>>>>>>>>>

#range of MYs

MYs = list(range(startMY,endMY,1))


#add the ages
calcs=pd.DataFrame(pd.Series(sorted(MYs*maxage),name='MY')) #creates a column of data showing startMY maxage times through endMY maxage times

#now join that column with a column showing ages 0 through maxage for each MY
calcs=calcs.join(
    pd.Series(list(range(0,maxage))*(endMY-startMY+1),name='age'),
                  how='inner')

#add the CYs of those ages - this simply adds the age to the MY to get the CY, so MY2010 plus age 3 gives CY 2013 (i.e., CY2013 has some 3 year old MY2010 vehicles)
calcs=calcs.join(
    pd.Series(calcs.MY+calcs.age,name='CY'),
                how='inner')

#add sales - this joins everything denoted as Sales (see inputs above) to the growing calcs DataFrame
calcs=calcs.join(Sales, on='MY')

#add CO2 & electicity emission factors - this joins everything denoted as CO2Analysis to the growing calcs DataFrame
calcs=calcs.join(CO2Analysis, on='MY')

#Create VMT matrix
vmt_Matrix=pd.DataFrame(pd.Series(range(startCY, endCY + maxage), name='CYs'))

#Create the factors which are used for VMT growth (**=exp)
#Note that startCY is set to 2009 since the baseVMT values are from 2009, therefore VMT growth is relative to 2009
VMTfactors= ([(1 + vmt_GrowthRate)] * len(vmt_Matrix.CYs)) ** (vmt_Matrix.CYs - startCY) #this calculates 1+vmt_GrowthRate a total of vmt_Matrix.CYs times (i.e., startCY through endCY+maxage times) and raises each to the exponent shown)

#List with the range of ages
age=range(0,maxage)

#Define a function that grows the VMT and multiplies by survival
def generate_VMT(baseVMT,vsurvival):
    vmt_Matrix.Base= baseVMT * vsurvival
    #print(vmt_Matrix.Base)
    VMTsched=[]
    for i in range(0,(len(VMTfactors))):
        VMTa=pd.Series(VMTfactors[i] * vmt_Matrix.Base, name=vmt_Matrix.CYs[i], index=age)
        if i==0:
            VMTsched=pd.DataFrame(VMTa, index=age)
            VMTsched.insert(0,'age',range(0,maxage))
        else:
            VMTsched=VMTsched.join(VMTa)
    return(VMTsched)

#Apply the function, and send the survival weighted VMT to a CSV
carVMT=generate_VMT(MAR_in.MAR_car,survival.surv_car)
truckVMT=generate_VMT(MAR_in.MAR_truck,survival.surv_truck)

def reshape_VMT(vtype,outVMT):
    #Reshape the VMT to match the form of the calculation matrix
    outVMTr=pd.core.reshape.melt(outVMT, id_vars=['age'])
    #print("This is reshaped ref case {0} VMT:\n".format(vtype), outVMTr)
    outVMTr.columns=['age', 'CY', 'miPerYear_'+vtype+'_ref']
    #create a joining id so that the calcs matrix can join the VMT
    outVMTr['MY']=outVMTr.CY-outVMTr.age
    outVMTr['id']=outVMTr.MY.map(str)+"-"+outVMTr.age.map(str)
    outVMTr=outVMTr.set_index('id')
    outVMTr=outVMTr.drop(['age','MY','CY'], axis = 1)
    calcs['id']=calcs.MY.map(str)+"-"+calcs.age.map(str)
    calcs2=calcs.set_index('id')
    ###Join VMT
    calcs2=calcs2.join(outVMTr)
    return(calcs2)

#Apply the Function
calcs=reshape_VMT('car',carVMT)
calcs=reshape_VMT('truck',truckVMT)

#create population matrix - one matrix for cars and one for trucks
# this population code is incomplete and it's not clear where Ari was headed with this
'''
populationC=pd.DataFrame(index=range(0,maxage))
populationT=pd.DataFrame(index=range(0,maxage))
populationC['age']=range(0,maxage)
populationT['age']=range(0,maxage)
populationC.insert(1,basepopyear,pop_in.pop_Car/survival.surv_Car)
populationT.insert(1,basepopyear,pop_in.pop_Truck/survival.surv_Truck)
cyears=range(basepopyear+1,endCY+1)

#A loop that generates the sales matrix.  This does not include survival. Survival is handled in the VMTs.
Sales['MY']=Sales.index
for i in range(0,len(cyears)):
        popC=pd.Series(name=cyears[i],index=age)
        popT=pd.Series(name=cyears[i],index=age)
        for r in age:
            if r==0:
                popC[r]=Sales.sales_Car_classic[[l for l,x in enumerate(Sales['MY']) if x == cyears[i]]]
                popT[r]=Sales.sales_Truck_classic[[l for l,x in enumerate(Sales['MY']) if x == cyears[i]]]
            else:
                popC[r]=populationC[cyears[i]-1][r-1]
                popT[r]=populationT[cyears[i]-1][r-1]
        populationC=populationC.join(popC)
        populationT=populationT.join(popT)
'''
#>>>>>>>>>> Cost per mile section >>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#join fuel Prices
fuelprice = pd.DataFrame(Fuel_Price, columns=['CY', '$_perGal_retail_diesel_'+case_ref, '$_perGal_pretax_diesel_'+case_ref,
                                              '$_perGal_retail_gasoline_'+case_ref, '$_perGal_pretax_gasoline_'+case_ref,
                                              '$_perkWh_retail_elec', '$_perkWh_pretax_elec'])
calcs=pd.merge(left=calcs, right=fuelprice, on='CY')
calcs=calcs.set_index(calcs.MY.map(str)+"-"+calcs.age.map(str))

#Convert to MPG for cost per mile - Should this be retail, even though the CO2 is E0?  I don't think so - thats mixing and matching.
#Setting to E0 CO2 density, so that the energy matches.....?
# This corrects to 5-cycle era MPG since the gap captures, in part, the energy content of 5-cycle era fuel
# and the CO2_gpm_road values used here include the gap
# NOTE-20160814: double check this - equations use road CO2
indoleneCO2density=Energy.ix['CO2_gpG_Indolene']['Value']
calcs['mpg_retail_car_ref']= indoleneCO2density / calcs['gpm_CO2_road_car_ref']
calcs['mpg_retail_truck_ref']=indoleneCO2density/calcs['gpm_CO2_road_truck_ref']
calcs['mpg_retail_car_control']=indoleneCO2density/calcs['gpm_CO2_road_car_control']
calcs['mpg_retail_truck_control']=indoleneCO2density/calcs['gpm_CO2_road_truck_control']

#Calculate Cost per mile assuming all liquid is gasoline
calcs['$_perMi_retail_gasoline_car_ref']=1/calcs['mpg_retail_car_ref'] * calcs['$_perGal_retail_gasoline_'+case_ref]
calcs['$_perMi_retail_gasoline_truck_ref']=1/calcs['mpg_retail_truck_ref']*calcs['$_perGal_retail_gasoline_'+case_ref]
calcs['$_perMi_retail_gasoline_car_control']=1/calcs['mpg_retail_car_control']*calcs['$_perGal_retail_gasoline_'+case_ref]
calcs['$_perMi_retail_gasoline_truck_control']=1/calcs['mpg_retail_truck_control']*calcs['$_perGal_retail_gasoline_'+case_ref]

#Calculate Cost per mile for electric
calcs['$_perMi_retail_elec_car_ref']=calcs['kWhPerMi_road_car_ref']*calcs['$_perkWh_retail_elec']
calcs['$_perMi_retail_elec_truck_ref']=calcs['kWhPerMi_road_truck_ref']*calcs['$_perkWh_retail_elec']
calcs['$_perMi_retail_elec_car_control']=calcs['kWhPerMi_road_car_control']*calcs['$_perkWh_retail_elec']
calcs['$_perMi_retail_elec_truck_control']=calcs['kWhPerMi_road_truck_control']*calcs['$_perkWh_retail_elec']

#Calculate Combined Cost per Mile
calcs['$_perMi_retail_car_ref']=calcs['$_perMi_retail_gasoline_car_ref']+calcs['$_perMi_retail_elec_car_ref']
calcs['$_perMi_retail_truck_ref']=calcs['$_perMi_retail_gasoline_truck_ref']+calcs['$_perMi_retail_elec_truck_ref']
calcs['$_perMi_retail_car_control']=calcs['$_perMi_retail_gasoline_car_control']+calcs['$_perMi_retail_elec_car_control']
calcs['$_perMi_retail_truck_control']=calcs['$_perMi_retail_gasoline_truck_control']+calcs['$_perMi_retail_elec_truck_control']
print ("Cost per mile complete")

#>>>>>>>>>>>>> Rebound VMT impact and VMT schedules including rebound (policy VMT schedule) >>>>>>>>>>>>>>>>>>>>>>>

#Calculate Rebound Impact
calcs['vmtincreasecar']=(calcs['$_perMi_retail_car_control']-calcs['$_perMi_retail_car_ref'])/calcs['$_perMi_retail_car_ref']*General.ix['Rebound Effect']['Value']
calcs['vmtincreasetruck']=(calcs['$_perMi_retail_truck_control']-calcs['$_perMi_retail_truck_ref'])/calcs['$_perMi_retail_truck_ref']*General.ix['Rebound Effect']['Value']

#Calculate Policy VMT Schedule
calcs['miPerYear_car_control']=calcs['vmtincreasecar']*calcs['miPerYear_car_ref']+calcs['miPerYear_car_ref']
calcs['miPerYear_truck_control']=calcs['vmtincreasetruck']*calcs['miPerYear_truck_ref']+calcs['miPerYear_truck_ref']
calcs['miPerYear_trucktocar_control']=calcs['vmtincreasecar']*calcs['miPerYear_truck_ref']+calcs['miPerYear_truck_ref']

#Multiply Sales by Miles for Total VMT - Note that cars in the policy case have a specific special case.
calcs['vmt_car_ref']=calcs['sales_car_MOVES']*calcs['miPerYear_car_ref']+calcs['sales_trucktocar']*calcs['miPerYear_truck_ref']
calcs['vmt_truck_ref']=calcs['sales_truck_CAFE']*calcs['miPerYear_truck_ref']
calcs['vmt_car_control']=calcs['sales_car_MOVES']*calcs['miPerYear_car_control']+calcs['sales_trucktocar']*calcs['miPerYear_trucktocar_control']
calcs['vmt_truck_control']=calcs['sales_truck_CAFE']*calcs['miPerYear_truck_control']
calcs['vmt_car_rebound']=calcs['vmt_car_control']-calcs['vmt_car_ref']
calcs['vmt_truck_rebound']=calcs['vmt_truck_control']-calcs['vmt_truck_ref']
calcs['vmt_fleet_rebound']=calcs['vmt_car_rebound']+calcs['vmt_truck_rebound']

#>>>>>>>>>>>>>>>>>> Fuel consumption section >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#Calculate Total Fuel Consumption (indolene)
indoleneCO2density=Energy.ix['CO2_gpG_Indolene']['Value']
calcs['gal_E0_car_ref']=calcs['vmt_car_ref']/ (indoleneCO2density / calcs['gpm_CO2_road_car_ref'])
calcs['gal_E0_car_control']=calcs['vmt_car_control']/ (indoleneCO2density / calcs['gpm_CO2_road_car_control'])
calcs['gal_E0_truck_ref']=calcs['vmt_truck_ref']/ (indoleneCO2density / calcs['gpm_CO2_road_truck_ref'])
calcs['gal_E0_truck_control']=calcs['vmt_truck_control']/ (indoleneCO2density / calcs['gpm_CO2_road_truck_control'])
calcs['gal_E0_car_delta']=calcs['gal_E0_car_control']-calcs['gal_E0_car_ref']
calcs['gal_E0_truck_delta']=calcs['gal_E0_truck_control']-calcs['gal_E0_truck_ref']
calcs['gal_E0_fleet_delta']=calcs['gal_E0_truck_delta']+calcs['gal_E0_car_delta']

#Calculate Total Fuel Consumption (Retail Gallons)
ratioenergy=Energy.ix['BTUpG_Indolene']['Value']/Energy.ix['BTUpG_GasolineRetail2025']['Value']
calcs['gal_retail_car_ref']=ratioenergy*calcs['gal_E0_car_ref']
calcs['gal_retail_car_control']=ratioenergy*calcs['gal_E0_car_control']
calcs['gal_retail_truck_ref']=ratioenergy*calcs['gal_E0_truck_ref']
calcs['gal_retail_truck_control']=ratioenergy*calcs['gal_E0_truck_control']
calcs['gal_retail_car_delta']=calcs['gal_retail_car_control']-calcs['gal_retail_car_ref']
calcs['gal_retail_truck_delta']=calcs['gal_retail_truck_control']-calcs['gal_retail_truck_ref']
calcs['gal_retail_fleet_delta']=calcs['gal_retail_truck_delta']+calcs['gal_retail_car_delta']

#>>>>>>>>>>>>> Monetized fuel section >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#Calculate Total Retail Fuel Savings (Retail Dollars)

calcs['$_retail_petrol_car_delta']=calcs['gal_retail_car_delta']*calcs['$_perGal_retail_gasoline_'+case_ref]
calcs['$_retail_petrol_truck_delta']=calcs['gal_retail_truck_delta']*calcs['$_perGal_retail_gasoline_'+case_ref]
calcs['$_retail_petrol_fleet_delta']=calcs['gal_retail_fleet_delta']*calcs['$_perGal_retail_gasoline_'+case_ref]

#Calculate Total Pre-tax Fuel Savings (Pre-tax Dollars)
calcs['$_pretax_petrol_car_delta']=calcs['gal_retail_car_delta']*calcs['$_perGal_pretax_gasoline_'+case_ref]
calcs['$_pretax_petrol_truck_delta']=calcs['gal_retail_truck_delta']*calcs['$_perGal_pretax_gasoline_'+case_ref]
calcs['$_pretax_petrol_fleet_delta']=calcs['gal_retail_fleet_delta']*calcs['$_perGal_pretax_gasoline_'+case_ref]

#>>>>>>>>>>>>> Rebound miles driven (to get rebound fuel consumption and increased fuel costs) >>>>>>>>>>>>>>>>>

#Calculate Rebound VMT
#calcs['carreboundVMT']=calcs['vmt_car_control']-calcs['vmt_car_ref']
#calcs['truckreboundVMT']=calcs['vmt_truck_control']-calcs['vmt_truck_ref']
#calcs['vmt_fleet_rebound']=calcs['carreboundVMT']+calcs['truckreboundVMT']

#Calculate Rebound Fuel Consumption (Retail gallons) and Costs (Retail dollars)

calcs['gal_retail_fleet_rebound']= (calcs['vmt_car_rebound'] / (indoleneCO2density / calcs['gpm_CO2_road_car_control'])) + (calcs['vmt_truck_rebound'] / (indoleneCO2density / calcs['gpm_CO2_road_truck_control']))
calcs['$_retail_petrol_fleet_rebound']=ratioenergy*calcs['gal_retail_fleet_rebound']*calcs['$_perGal_retail_gasoline_'+case_ref]

#>>>>>>>>>>>> CO2 inventories (Metric tons) >>>>>>>>>>>>>>>>>>>>>>>>>>>>

#Calculate Tailpipe CO2 inventories in Tons
calcs['ton_CO2_road_car_ref']=(calcs['vmt_car_ref']*calcs['gpm_CO2_road_car_ref'])/million
calcs['ton_CO2_road_car_control']=(calcs['vmt_car_control']* (calcs['gpm_CO2_road_car_control']))/million
calcs['ton_CO2_road_truck_ref']=(calcs['vmt_truck_ref']*(calcs['gpm_CO2_road_truck_ref']))/million
calcs['ton_CO2_road_truck_control']=(calcs['vmt_truck_control']* (calcs['gpm_CO2_road_truck_control']))/million
calcs['ton_CO2_road_car_delta']=calcs['ton_CO2_road_car_control']-calcs['ton_CO2_road_car_ref']
calcs['ton_CO2_road_truck_delta']=calcs['ton_CO2_road_truck_control']-calcs['ton_CO2_road_truck_ref']
calcs['ton_CO2_road_fleet_delta']=calcs['ton_CO2_road_truck_delta']+calcs['ton_CO2_road_car_delta']

#Calculate Indirect A/C inventory - Remember, the general inventory is included in the onroad gap, and the delta should be captured in the CO2 numbers that run through the OMEGA model.
#Thus, this is just reference.
ACbyIndex=pd.DataFrame(calcs['MY']).join(AC, on='MY')
calcs['ton_CO2_road_ACefficiency_car_ref']=(calcs['vmt_car_ref']*ACbyIndex['acER_eff_car_ref'])/million
calcs['ton_CO2_road_ACefficiency_car_control']=(calcs['vmt_car_control']* (ACbyIndex['acER_eff_car_control']))/million
calcs['ton_CO2_road_ACefficiency_truck_ref']=(calcs['vmt_truck_ref']*(ACbyIndex['acER_eff_truck_ref']))/million
calcs['ton_CO2_road_ACefficiency_truck_control']=(calcs['vmt_truck_control']* (ACbyIndex['acER_eff_truck_control']))/million
calcs['ton_CO2_road_ACefficiency_car_delta']=calcs['ton_CO2_road_ACefficiency_car_control']-calcs['ton_CO2_road_ACefficiency_car_ref']
calcs['ton_CO2_road_ACefficiency_truck_delta']=calcs['ton_CO2_road_ACefficiency_truck_control']-calcs['ton_CO2_road_ACefficiency_truck_ref']
calcs['ton_CO2_road_ACefficiency_fleet_delta']=calcs['ton_CO2_road_ACefficiency_truck_delta']+calcs['ton_CO2_road_ACefficiency_car_delta']

#Calculate Indirect A/C EO gasoline impact - remember, this is part of the total already calculated above
#and is just shown for example purposes; and remember that indoleneCO2density=Energy['CO2_gpG_Indolene']
calcs['gal_E0_ACefficiency_car_ref']= calcs['ton_CO2_road_ACefficiency_car_ref'] / indoleneCO2density
calcs['gal_E0_ACefficiency_car_control']= calcs['ton_CO2_road_ACefficiency_car_control'] / indoleneCO2density
calcs['gal_E0_ACefficiency_truck_ref']= calcs['ton_CO2_road_ACefficiency_truck_ref'] / indoleneCO2density
calcs['gal_E0_ACefficiency_truck_control']= calcs['ton_CO2_road_ACefficiency_truck_control'] / indoleneCO2density
calcs['gal_E0_ACefficiency_car_delta']= calcs['ton_CO2_road_ACefficiency_car_delta'] / indoleneCO2density
calcs['gal_E0_ACefficiency_truck_delta']= calcs['ton_CO2_road_ACefficiency_truck_delta'] / indoleneCO2density
calcs['gal_E0_ACefficiency_fleet_delta']= calcs['ton_CO2_road_ACefficiency_fleet_delta'] / indoleneCO2density

#>>>>>>>>>>>>> Direct AC impacts >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#Calculate Direct A/C impact in Tons
#HFC is based on reference VMT only - so that the same leakage rate produces the same totals in both cases.
#Not clear that driving more will cause people to refill refrigerant more (not necessarily increasing durability?)
calcs['ton_CO2e_road_ACleakage_car_ref']=calcs['vmt_car_ref']*ACbyIndex['acER_leakage_car_ref']/million
calcs['ton_CO2e_road_ACleakage_car_control']=calcs['vmt_car_ref']* (ACbyIndex['acER_leakage_car_control'])/million
calcs['ton_CO2e_road_ACleakage_truck_ref']=calcs['vmt_truck_ref']*(ACbyIndex['acER_leakage_truck_ref'])/million
calcs['ton_CO2e_road_ACleakage_truck_control']=calcs['vmt_truck_ref']* (ACbyIndex['acER_leakage_truck_control'])/million
calcs['ton_CO2e_road_ACleakage_car_delta']=calcs['ton_CO2e_road_ACleakage_car_control']-calcs['ton_CO2e_road_ACleakage_car_ref']
calcs['ton_CO2e_road_ACleakage_truck_delta']=calcs['ton_CO2e_road_ACleakage_truck_control']-calcs['ton_CO2e_road_ACleakage_truck_ref']
calcs['ton_CO2e_road_ACleakage_fleet_delta']=calcs['ton_CO2e_road_ACleakage_truck_delta']+calcs['ton_CO2e_road_ACleakage_car_delta']

#>>>>>>>>>>>>>>>>>>>> Electricity consumption and costs >>>>>>>>>>>>>>>>>>

#Calculate Electricity Consumed at vehicle
calcs['kWh_retail_car_ref']=calcs['vmt_car_ref']*calcs['kWhPerMi_road_car_ref']
calcs['kWh_retail_car_control']=calcs['vmt_car_control']*calcs['kWhPerMi_road_car_control']
calcs['kWh_retail_truck_ref']=calcs['vmt_truck_ref']*calcs['kWhPerMi_road_truck_ref']
calcs['kWh_retail_truck_control']=calcs['vmt_truck_control']*calcs['kWhPerMi_road_truck_control']
calcs['kWh_retail_car_delta']=calcs['kWh_retail_car_control']-calcs['kWh_retail_car_ref']
calcs['kWh_retail_truck_delta']=calcs['kWh_retail_truck_control']-calcs['kWh_retail_truck_ref']
calcs['kWh_retail_fleet_delta']=calcs['kWh_retail_truck_delta']+calcs['kWh_retail_car_delta']

#Calculate Electricity Consumed at Plant
tranloss=General.ix['Transmission Loss']['Value']
calcs['kWh_plant_car_ref']=calcs['kWh_retail_car_ref']/tranloss
calcs['kWh_plant_car_control']=calcs['kWh_retail_car_control'] / tranloss
calcs['kWh_plant_truck_ref']=calcs['kWh_retail_truck_ref'] / tranloss
calcs['kWh_plant_truck_control']=calcs['kWh_retail_truck_control']/ tranloss
calcs['kWh_plant_car_delta']=calcs['kWh_retail_car_delta']/tranloss
calcs['kWh_plant_truck_delta']=calcs['kWh_retail_truck_delta']/ tranloss
calcs['kWh_plant_fleet_delta']=calcs['kWh_retail_fleet_delta']/ tranloss

#Calculate Cost of Electricity Emissions
calcs['$_retail_elec_car_delta']=calcs['kWh_plant_car_delta']*calcs['$_perkWh_retail_elec']
calcs['$_retail_elec_truck_delta']=calcs['kWh_plant_truck_delta']*calcs['$_perkWh_retail_elec']
calcs['$_retail_elec_fleet_delta']=calcs['kWh_plant_fleet_delta']*calcs['$_perkWh_retail_elec']

#>>>>>>>>>>>>>>>>>>> Electricity consumption impacts on criteria pollutants >>>>>>>>>>>>>>>>>>>>>>>>>

#Calculate Electricity Emissions
ElecEFbyIndex=pd.DataFrame(calcs['CY']).join(ElecEF, on='CY')

#define a function for electricity emissions that we want in short tons
def elecemissions(pol, elec_plant_ef):
    b='uston_up_' + pol
    z1=pd.DataFrame((pd.Series(calcs['kWh_plant_car_ref'] * elec_plant_ef / million * metric_to_uston, name= b + '_elec_car_ref')))
    z=z1.join([
    pd.DataFrame(pd.Series(calcs['kWh_plant_car_control'] * elec_plant_ef / million * metric_to_uston, name= b + '_elec_car_control')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_truck_ref'] * elec_plant_ef / million * metric_to_uston, name= b + '_elec_truck_ref')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_truck_control'] * elec_plant_ef / million * metric_to_uston, name= b + '_elec_truck_control')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_car_delta'] * elec_plant_ef / million * metric_to_uston, name= b + '_elec_car_delta')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_truck_delta'] * elec_plant_ef / million * metric_to_uston, name= b + '_elec_truck_delta')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_fleet_delta'] * elec_plant_ef / million * metric_to_uston, name= b + '_elec_fleet_delta')),
    ])
    return(z)

#pd.DataFrame()

#Run Electricity Function for all pollutants, combine together, also add emissions for later rates.
calcs=calcs.join([
elecemissions('VOC', ElecEFbyIndex['gpkWh_VOC']),
elecemissions('CO', ElecEFbyIndex['gpkWh_CO']),
elecemissions('NOx', ElecEFbyIndex['gpkWh_NOx']),
elecemissions('PM25', ElecEFbyIndex['gpkWh_PM25']),
elecemissions('SOx', ElecEFbyIndex['gpkWh_SOx']),
#elecemissions('CO2',ElecEFbyIndex['gpkWh_CO2'])/metric_to_uston,
#elecemissions('N2O',ElecEFbyIndex['gpkWh_N2O'])/metric_to_uston,
#elecemissions('CH4',ElecEFbyIndex['gpkWh_CH4'])/metric_to_uston,
elecemissions('Benzene', ElecEFbyIndex['gpkWh_Benzene']),
elecemissions('13Butadiene', ElecEFbyIndex['gpkWh_13Butadiene']),
elecemissions('Formaldehyde', ElecEFbyIndex['gpkWh_Formaldehyde']),
elecemissions('Acetaldehyde', ElecEFbyIndex['gpkWh_Acetaldehyde']),
elecemissions('Acrolein', ElecEFbyIndex['gpkWh_Acrolein'])
])

#define a function for electricity emissions that we want in metric tons
def elecemissions2(pol, elec_plant_ef):
    b = 'ton_up_' + pol
    z1=pd.DataFrame((pd.Series(calcs['kWh_plant_car_ref'] * elec_plant_ef / million, name= b + '_elec_car_ref')))
    z=z1.join([
    pd.DataFrame(pd.Series(calcs['kWh_plant_car_control'] * elec_plant_ef / million, name= b + '_elec_car_control')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_truck_ref'] * elec_plant_ef / million, name= b + '_elec_truck_ref')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_truck_control'] * elec_plant_ef / million, name= b + '_elec_truck_control')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_car_delta'] * elec_plant_ef / million, name= b + '_elec_car_delta')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_truck_delta'] * elec_plant_ef / million, name= b + '_elec_truck_delta')),
    pd.DataFrame(pd.Series(calcs['kWh_plant_fleet_delta'] * elec_plant_ef / million, name= b + '_elec_fleet_delta')),
    ])
    return(z)

calcs=calcs.join([
elecemissions2('CO2', ElecEFbyIndex['gpkWh_CO2']),
elecemissions2('N2O', ElecEFbyIndex['gpkWh_N2O']),
elecemissions2('CH4', ElecEFbyIndex['gpkWh_CH4'])
])

#>>>>>>>>>>>>>>>>>>>>> Critieria emissions from rebound miles (assumes all are driven using gasoline) >>>>>>>>>>>>>>>>>>>>>>>>>

#Calculate Emissions due to Rebound -the rebound VMT is by new definition, and the rates are by old.
calcs['vmt_MOVES_car_rebound']=calcs['sales_car_MOVES']*calcs['miPerYear_car_control']-calcs['sales_car_MOVES']*calcs['miPerYear_car_ref']
calcs['vmt_MOVES_truck_rebound']=calcs['sales_truck_CAFE']*calcs['miPerYear_truck_control']-calcs['sales_truck_CAFE']*calcs['miPerYear_truck_ref']
calcs['vmt_trucktocar_rebound']=calcs['sales_trucktocar']*calcs['miPerYear_trucktocar_control']-calcs['sales_trucktocar']*calcs['miPerYear_truck_ref']

rebframe=pd.DataFrame([

#Cars
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['VOC_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['VOC_gpm_Car']) / million * metric_to_uston, name='uston_down_VOC_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['CO_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['CO_gpm_Car']) / million * metric_to_uston, name='uston_down_CO_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['NOx_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['NOx_gpm_Car']) / million * metric_to_uston, name='uston_down_NOx_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['PM2.5_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['PM2.5_gpm_Car']) / million * metric_to_uston, name='uston_down_PM25_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['Benzene_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['Benzene_gpm_Car']) / million * metric_to_uston, name='uston_down_Benzene_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['1,3 Butadiene_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['1,3 Butadiene_gpm_Car']) / million * metric_to_uston, name='uston_down_13Butadiene_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['Formaldehyde_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['Formaldehyde_gpm_Car']) / million * metric_to_uston, name='uston_down_Formaldehyde_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['Acetaldehyde_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['Acetaldehyde_gpm_Car']) / million * metric_to_uston, name='uston_down_Acetaldehyde_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['Acrolein_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['Acrolein_gpm_Car']) / million * metric_to_uston, name='uston_down_Acrolein_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['CH4_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['CH4_gpm_Car'])/million,name='uston_down_CH4_car_rebound'),
pd.Series((calcs['vmt_trucktocar_rebound']*VehicleEF['N20_gpm_Truck']+calcs['vmt_MOVES_car_rebound']*VehicleEF['N20_gpm_Car'])/million,name='uston_down_N2O_car_rebound'),

#Different Algorithm for SO2
pd.Series(-calcs['gal_E0_car_delta'] * VehicleEF['SO2_gpG_Car'] / million * metric_to_uston, name='uston_down_SOx_car_rebound'),

#Trucks
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['VOC_gpm_Truck'] / million * metric_to_uston, name='uston_down_VOC_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['CO_gpm_Truck'] / million * metric_to_uston, name='uston_down_CO_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['NOx_gpm_Truck'] / million * metric_to_uston, name='uston_down_NOx_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['PM2.5_gpm_Truck'] / million * metric_to_uston, name='uston_down_PM25_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['Benzene_gpm_Truck'] / million * metric_to_uston, name='uston_down_Benzene_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['1,3 Butadiene_gpm_Truck'] / million * metric_to_uston, name='uston_down_13Butadiene_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['Formaldehyde_gpm_Truck'] / million * metric_to_uston, name='uston_down_Formaldehyde_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['Acetaldehyde_gpm_Truck'] / million * metric_to_uston, name='uston_down_Acetaldehyde_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound'] * VehicleEF['Acrolein_gpm_Truck'] / million * metric_to_uston, name='uston_down_Acrolein_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound']*VehicleEF['CH4_gpm_Truck']/million,name='uston_down_CH4_truck_rebound'),
pd.Series(calcs['vmt_MOVES_truck_rebound']*VehicleEF['N20_gpm_Truck']/million,name='uston_down_N2O_truck_rebound'),

#Different Algorithm for SO2
pd.Series(-calcs['gal_E0_truck_delta'] * VehicleEF['SO2_gpG_Truck'] / million * metric_to_uston, name='uston_down_SOx_truck_rebound')
])

#Transpose Output and join to calcs matrix
rebframeT=rebframe.T
calcs=calcs.join(rebframeT)

#>>>>>>>>>>>>>>>>>>>>>> Calculate Upstream Fuel Production Emissions >>>>>>>>>>>>>>>>>>>>>>>>>

#Get Fuel in Energy
calcs['millionBTU_refinery_delta']=(calcs['gal_E0_car_delta']+calcs['gal_E0_truck_delta'])*Energy.ix['BTUpG_Indolene']['Value']/million

#Create a dataframe that makes the EFs into by index rather than by CY
RefineryEFbyIndex=pd.DataFrame(calcs['CY']).join(RefineryEF, on='CY')

#Create a DF for the inventories
refineryemissions=pd.DataFrame(calcs['CY'])

#A loop that goes through the emission rates and generates inventories (more elegant than approach above)
for i in (range(1, len(RefineryEFbyIndex.columns))):
    namer='uston_up_' + RefineryEFbyIndex.columns[i] + '_refinery_delta'
    refineryemissions[namer]=pd.Series(calcs['millionBTU_refinery_delta'] * RefineryEFbyIndex[RefineryEFbyIndex.columns[i]] / million * metric_to_uston, name=(namer))

#Join it back to the main calculation matrix.
refineryemissions=refineryemissions.drop('CY', axis=1)
calcs=calcs.join(refineryemissions)

#Make the CO2 pollutants into Metric tons
calcs['ton_up_CH4_refinery_delta']= calcs['uston_up_CH4_refinery_delta'] / metric_to_uston
calcs['ton_up_N2O_refinery_delta']= calcs['uston_up_N2O_refinery_delta'] / metric_to_uston
calcs['ton_up_CO2_refinery_delta']= calcs['uston_up_CO2_refinery_delta'] / metric_to_uston

#Monetized Pollutants Upstream, these are dtype='object' so convert to float64 when used
EconYearbyIndex=pd.DataFrame(calcs['CY']).join(EconYear, on='CY')

def nonGHG_upBenefits(pol, rate, range):
    b='$_'+rate+'_'+range+'_up_'+pol
    c='uston_up_'+pol+'_refinery_delta'
    d='uston_up_'+pol+'_elec_fleet_delta'
    calcs[b]=pd.Series((calcs[c]+calcs[d])*EconYearbyIndex[pol+'_up_'+rate+'_'+range], dtype='float64')
    return(calcs)

def nonGHG_downBenefits(pol, rate, range):
    b='$_'+rate+'_'+range+'_down_'+pol
    c='uston_down_'+pol+'_car_rebound'
    d='uston_down_'+pol+'_truck_rebound'
    calcs[b]=pd.Series((calcs[c]+calcs[d])*EconYearbyIndex[pol+'_down_'+rate+'_'+range], dtype='float64')
    return(calcs)
#upstream, 3%
nonGHG_upBenefits('PM25', '0.03', 'low')
nonGHG_upBenefits('PM25', '0.03', 'mid')
nonGHG_upBenefits('PM25', '0.03', 'high')
nonGHG_upBenefits('NOx', '0.03', 'low')
nonGHG_upBenefits('NOx', '0.03', 'mid')
nonGHG_upBenefits('NOx', '0.03', 'high')
nonGHG_upBenefits('SOx', '0.03', 'low')
nonGHG_upBenefits('SOx', '0.03', 'mid')
nonGHG_upBenefits('SOx', '0.03', 'high')
#upstream, 7%
nonGHG_upBenefits('PM25', '0.07', 'low')
nonGHG_upBenefits('PM25', '0.07', 'mid')
nonGHG_upBenefits('PM25', '0.07', 'high')
nonGHG_upBenefits('NOx', '0.07', 'low')
nonGHG_upBenefits('NOx', '0.07', 'mid')
nonGHG_upBenefits('NOx', '0.07', 'high')
nonGHG_upBenefits('SOx', '0.07', 'low')
nonGHG_upBenefits('SOx', '0.07', 'mid')
nonGHG_upBenefits('SOx', '0.07', 'high')
#downstream, 3%
nonGHG_downBenefits('PM25', '0.03', 'low')
nonGHG_downBenefits('PM25', '0.03', 'mid')
nonGHG_downBenefits('PM25', '0.03', 'high')
nonGHG_downBenefits('NOx', '0.03', 'low')
nonGHG_downBenefits('NOx', '0.03', 'mid')
nonGHG_downBenefits('NOx', '0.03', 'high')
nonGHG_downBenefits('SOx', '0.03', 'low')
nonGHG_downBenefits('SOx', '0.03', 'mid')
nonGHG_downBenefits('SOx', '0.03', 'high')
#downstream, 7%
nonGHG_downBenefits('PM25', '0.07', 'low')
nonGHG_downBenefits('PM25', '0.07', 'mid')
nonGHG_downBenefits('PM25', '0.07', 'high')
nonGHG_downBenefits('NOx', '0.07', 'low')
nonGHG_downBenefits('NOx', '0.07', 'mid')
nonGHG_downBenefits('NOx', '0.07', 'high')
nonGHG_downBenefits('SOx', '0.07', 'low')
nonGHG_downBenefits('SOx', '0.07', 'mid')
nonGHG_downBenefits('SOx', '0.07', 'high')

calcs['$_0.03_low_PM25']=calcs[['$_0.03_low_up_PM25', '$_0.03_low_down_PM25']].sum(axis=1)
calcs['$_0.03_low_NOx']=calcs[['$_0.03_low_up_NOx', '$_0.03_low_down_NOx']].sum(axis=1)
calcs['$_0.03_low_SOx']=calcs[['$_0.03_low_up_SOx', '$_0.03_low_down_SOx']].sum(axis=1)
calcs['$_0.03_mid_PM25']=calcs[['$_0.03_mid_up_PM25', '$_0.03_mid_down_PM25']].sum(axis=1)
calcs['$_0.03_mid_NOx']=calcs[['$_0.03_mid_up_NOx', '$_0.03_mid_down_NOx']].sum(axis=1)
calcs['$_0.03_mid_SOx']=calcs[['$_0.03_mid_up_SOx', '$_0.03_mid_down_SOx']].sum(axis=1)
calcs['$_0.03_high_PM25']=calcs[['$_0.03_high_up_PM25', '$_0.03_high_down_PM25']].sum(axis=1)
calcs['$_0.03_high_NOx']=calcs[['$_0.03_high_up_NOx', '$_0.03_high_down_NOx']].sum(axis=1)
calcs['$_0.03_high_SOx']=calcs[['$_0.03_high_up_SOx', '$_0.03_high_down_SOx']].sum(axis=1)
calcs['$_0.07_low_PM25']=calcs[['$_0.07_low_up_PM25', '$_0.07_low_down_PM25']].sum(axis=1)
calcs['$_0.07_low_NOx']=calcs[['$_0.07_low_up_NOx', '$_0.07_low_down_NOx']].sum(axis=1)
calcs['$_0.07_low_SOx']=calcs[['$_0.07_low_up_SOx', '$_0.07_low_down_SOx']].sum(axis=1)
calcs['$_0.07_mid_PM25']=calcs[['$_0.07_mid_up_PM25', '$_0.07_mid_down_PM25']].sum(axis=1)
calcs['$_0.07_mid_NOx']=calcs[['$_0.07_mid_up_NOx', '$_0.07_mid_down_NOx']].sum(axis=1)
calcs['$_0.07_mid_SOx']=calcs[['$_0.07_mid_up_SOx', '$_0.07_mid_down_SOx']].sum(axis=1)
calcs['$_0.07_high_PM25']=calcs[['$_0.07_high_up_PM25', '$_0.07_high_down_PM25']].sum(axis=1)
calcs['$_0.07_high_NOx']=calcs[['$_0.07_high_up_NOx', '$_0.07_high_down_NOx']].sum(axis=1)
calcs['$_0.07_high_SOx']=calcs[['$_0.07_high_up_SOx', '$_0.07_high_down_SOx']].sum(axis=1)

#Monetize SC_CO2, SC_CH4 and SCC N2O
CH4GWP=General.ix['CH4-GWP']['Value']
N2OGWP=General.ix['N2O-GWP']['Value']
HFCGWP=General.ix['HFC-GWP']['Value']

#Add up TP CO2e reductions
calcs['mmt_CO2e_road_fleet_delta']=(calcs['ton_CO2_road_fleet_delta']+(calcs['uston_down_CH4_car_rebound']+calcs[
    'uston_down_CH4_truck_rebound'])*CH4GWP+
                               (calcs['uston_down_N2O_car_rebound']+calcs['uston_down_N2O_truck_rebound'])*N2OGWP+calcs['ton_CO2e_road_ACleakage_fleet_delta']*HFCGWP)/million

#Add up all the CH4, N2O, CO2, HFC, do NOT use GWP to convert to CO2e but remember that HFC is in CO2e - convert Tons to MMT
calcs['mmt_CH4_total_delta']=(calcs['ton_up_CH4_refinery_delta']+calcs['ton_up_CH4_elec_fleet_delta']+calcs['uston_down_CH4_car_rebound']+ calcs['uston_down_CH4_truck_rebound'])/million
calcs['mmt_N2O_total_delta']=(calcs['ton_up_N2O_refinery_delta']+calcs['ton_up_N2O_elec_fleet_delta']+calcs['uston_down_N2O_car_rebound']+calcs['uston_down_N2O_truck_rebound'])/million
calcs['mmt_CO2_total_delta']=(calcs['ton_up_CO2_refinery_delta']+calcs['ton_up_CO2_elec_fleet_delta']+calcs['ton_CO2_road_fleet_delta'])/million
calcs['mmt_HFC_total_delta']=calcs['ton_CO2e_road_ACleakage_fleet_delta']*HFCGWP/million
calcs['mmt_CO2e_total_delta']=calcs['mmt_HFC_total_delta']/HFCGWP+calcs['mmt_CO2_total_delta']+calcs['mmt_N2O_total_delta']*N2OGWP+calcs['mmt_CH4_total_delta']*CH4GWP
#Apply each of 4 SCC streams to each of 3 GHGs
#calcs['totdeltaCO2e_MMT']=calcs['totdeltaHFCCO2e_MMT']+calcs['totdeltaCO2_MMT']+calcs['totdeltaN2OCO2e_MMT']+calcs['totdeltaCH4CO2e_MMT']

def scc(ghg, rate):
    b='$_SC' + ghg + '_' + rate
    calcs[b]=pd.Series(calcs['mmt_'+ghg+'_total_delta']*EconYearbyIndex['SC'+ghg+'_'+rate]*million, dtype='float64')
    return(calcs)

scc('CO2', '5')
scc('CO2', '3')
scc('CO2', '25')
scc('CO2', '395')
scc('CH4', '5')
scc('CH4', '3')
scc('CH4', '25')
scc('CH4', '395')
scc('N2O', '5')
scc('N2O', '3')
scc('N2O', '25')
scc('N2O', '395')

calcs['$_SCGHG_5']=calcs[['$_SCCO2_5', '$_SCCH4_5', '$_SCN2O_5']].sum(axis=1)
calcs['$_SCGHG_3']=calcs[['$_SCCO2_3', '$_SCCH4_3', '$_SCN2O_3']].sum(axis=1)
calcs['$_SCGHG_25']=calcs[['$_SCCO2_25', '$_SCCH4_25', '$_SCN2O_25']].sum(axis=1)
calcs['$_SCGHG_395']=calcs[['$_SCCO2_395', '$_SCCH4_395', '$_SCN2O_395']].sum(axis=1)

def TravelBenefits(veh):
    b='$_TravelBenefit_'+veh
    c=calcs['$_perMi_retail_'+veh+'_control']
    d=calcs['$_perMi_retail_'+veh+'_ref']
    e=calcs['vmt_'+veh+'_rebound']
    calcs[b]=-(c*e) + 0.5*e*(c-d)
    return(calcs)

TravelBenefits('car')
TravelBenefits('truck')
calcs['$_TravelBenefit_fleet']=calcs[['$_TravelBenefit_car', '$_TravelBenefit_truck']].sum(axis=1)

#Monetize Congestion, Accidents, Noise - static, doesn't change by CY
calcs['$_congestion_car']=calcs['vmt_car_rebound']*General.ix['Congestion']['Value']
calcs['$_congestion_truck']=calcs['vmt_truck_rebound']*General.ix['Congestion']['Value']
calcs['$_crashes_car']=calcs['vmt_car_rebound']*General.ix['Accidents']['Value']
calcs['$_crashes_truck']=calcs['vmt_truck_rebound']*General.ix['Accidents']['Value']
calcs['$_noise_car']=calcs['vmt_car_rebound']*General.ix['Noise']['Value']
calcs['$_noise_truck']=calcs['vmt_truck_rebound']*General.ix['Noise']['Value']
calcs['$_NCCP_car']=calcs[['$_congestion_car','$_crashes_car','$_noise_car']].sum(axis=1)
calcs['$_NCCP_truck']=calcs[['$_congestion_truck','$_crashes_truck','$_noise_truck']].sum(axis=1)
calcs['$_NCCP_fleet']=calcs[['$_NCCP_car','$_NCCP_truck']].sum(axis=1)

#Monetize Energy Security impacts
calcs['$_EnergySecurity']=pd.Series(EconYearbyIndex['Energy Security']*calcs['gal_E0_fleet_delta']*foreignOilFrac, dtype='float64')

#Monetize Refueling impacts
def refueling(veh):
    b='$_Refueling_'+veh
    c=General.ix['Fuel Tank volume']['Value']
    d=General.ix['Average tank refill']['Value']
    e=General.ix['Time per refueling']['Value']/60
    f='gal_retail_'+veh+'_delta'
    g=General.ix['Benefits discounted']['Value']
    calcs[b]=pd.Series((calcs[f]/(c*d))*e*(1-g)*EconYearbyIndex['wage_rate'], dtype='float64')
    return(calcs)

refueling('car')
refueling('truck')
calcs['$_Refueling_fleet']=calcs[['$_Refueling_car', '$_Refueling_truck']].sum(axis=1)

#Mass related fatalities
Mass=Mass.ix[:, 0:5]
calcsTemp=calcs.copy()
calcsTemp=calcsTemp.merge(Mass, on='MY', how='left')
calcsTemp=calcsTemp.set_index(calcsTemp.MY.map(str)+'-'+calcsTemp.age.map(str))
calcs['fatalities_car_ref']=calcs['vmt_car_ref']*calcsTemp['car_ref']
calcs['fatalities_truck_ref']=calcs['vmt_truck_ref']*calcsTemp['truck_ref']
calcs['fatalities_car_control']=calcs['vmt_car_control']*calcsTemp['car_control']
calcs['fatalities_truck_control']=calcs['vmt_truck_control']*calcsTemp['truck_control']
calcs['fatalities_car_delta']=calcs['fatalities_car_control']-calcs['fatalities_car_ref']
calcs['fatalities_truck_delta']=calcs['fatalities_truck_control']-calcs['fatalities_truck_ref']
calcs['fatalities_fleet_delta']=calcs['fatalities_car_delta']+calcs['fatalities_truck_delta']

#Add Legacy Population for CY analysis (this is Ari's old message)
'''
#These population calcs provide populations of MY2014 and later vehicles for use in maintenance calcs
for i in range(startMY, endMY):
    yr=str(i)
    popC_temp=pd.DataFrame(populationC['age'])
    popT_temp=pd.DataFrame(populationT['age'])
    popC_temp['MY']=''
    popT_temp['MY']=''
    popC_temp['CY']=i
    popT_temp['CY']=i
    popC_temp['population_Car']=populationC[i]*survival['surv_car']
    popT_temp['population_Truck']=populationT[i]*survival['surv_truck']
    if i==startMY:
        populationC_reshaped=pd.DataFrame(popC_temp)
        populationT_reshaped=pd.DataFrame(popT_temp)
    else:
        populationC_reshaped=populationC_reshaped.append(popC_temp)
        populationT_reshaped=populationT_reshaped.append(popT_temp)

populationC_reshaped=populationC_reshaped.set_index(populationC_reshaped.MY.map(str)+'-'+populationC_reshaped.age.map(str))
populationT_reshaped=populationT_reshaped.set_index(populationT_reshaped.MY.map(str)+'-'+populationT_reshaped.age.map(str))
population=populationC_reshaped.join(populationT_reshaped['population_Truck'])
population['CY']=population[['MY','age']].sum(axis=1)
population=population[['MY', 'age', 'CY', 'population_Car', 'population_Truck']]
population.fillna(0, inplace=True)
population['population_fleet']=population[['population_Car','population_Truck']].sum(axis=1)
calcs=calcs.join(population['population_Car'])
calcs=calcs.join(population['population_Truck'])
calcs=calcs.join(population['population_fleet'])
'''

#Add FFV analysis Portion

#>>>>>>>>>>>>>>>>  Vehicle program costs >>>>>>>>>>>>>>>
#first, tech costs
vehcosts=inputFile.parse('VehCost', index_col=0)
vehcosts=vehcosts.ix[:,1:5]
calcs=calcs.join(vehcosts)
calcs.fillna(0, inplace=True)
calcs['$_perVeh_Tech_car_delta']=calcs['$_perVeh_Tech_car_control']-calcs['$_perVeh_Tech_car_ref']
calcs['$_perVeh_Tech_truck_delta']=calcs['$_perVeh_Tech_truck_control']-calcs['$_perVeh_Tech_truck_ref']
calcs['$_perVeh_Tech_fleet_delta']=(calcs[['$_perVeh_Tech_car_delta', 'sales_car_CAFE']].product(axis=1)+calcs[['$_perVeh_Tech_truck_delta','sales_truck_CAFE']].product(axis=1))/calcs[['sales_car_CAFE','sales_truck_CAFE']].sum(axis=1)
calcs['$_Tech_car_ref']=calcs[['sales_car_CAFE', '$_perVeh_Tech_car_ref']].product(axis=1)
calcs['$_Tech_car_control']=calcs[['sales_car_CAFE', '$_perVeh_Tech_car_control']].product(axis=1)
calcs['$_Tech_truck_ref']=calcs[['sales_truck_CAFE', '$_perVeh_Tech_truck_ref']].product(axis=1)
calcs['$_Tech_truck_control']=calcs[['sales_truck_CAFE', '$_perVeh_Tech_truck_control']].product(axis=1)
calcs['$_Tech_car_delta']=calcs['$_Tech_car_control']-calcs['$_Tech_car_ref']
calcs['$_Tech_truck_delta']=calcs['$_Tech_truck_control']-calcs['$_Tech_truck_ref']
calcs['$_Tech_fleet_delta']=calcs[['$_Tech_car_delta', '$_Tech_truck_delta']].sum(axis=1)

#second, maintenance costs
maintTechPens=inputFile.parse('MaintTechPens_Analysis', skiprows=1)
calcsTemp=pd.DataFrame(calcs, columns=['MY', 'age', 'vmt_car_ref', 'vmt_truck_ref', 'vmt_car_control', 'vmt_truck_control'])
maintTechPensbyIndex=calcsTemp.merge(maintTechPens, on='MY', how='left')
maintTechPensbyIndex=maintTechPensbyIndex.set_index(maintTechPensbyIndex.MY.map(str)+"-"+maintTechPensbyIndex.age.map(str))
maintPE=inputFile.parse('MaintMetricsPE', index_col='Tech')
# cpe=cost per event; mpe=miles per event
cpe_LRRT1=maintPE.ix['LRRT1']['CostPE']
cpe_LRRT2=maintPE.ix['LRRT2']['CostPE']
cpe_Diesel=maintPE.ix['Diesel']['CostPE']
cpe_EVoilChange=maintPE.ix['EV_OilChange']['CostPE']
cpe_EVairFilter=maintPE.ix['EV_AirFilter']['CostPE']
cpe_EVengineCoolant=maintPE.ix['EV_EngineCoolant']['CostPE']
cpe_EVsparkPlugs=maintPE.ix['EV_SparkPlugs']['CostPE']
cpe_EVbatteryCoolant=maintPE.ix['EV_BatteryCoolant']['CostPE']
cpe_EVbatteryCheck=maintPE.ix['EV_BatteryCheck']['CostPE']
cpe_PHEVbatteryCoolant=maintPE.ix['PHEV_BatteryCoolant']['CostPE']
cpe_PHEVbatteryCheck=maintPE.ix['EV_BatteryCheck']['CostPE']
mpe_LRRT1=maintPE.ix['LRRT1']['MilesPE']
mpe_LRRT2=maintPE.ix['LRRT2']['MilesPE']
mpe_Diesel=maintPE.ix['Diesel']['MilesPE']
mpe_EVoilChange=maintPE.ix['EV_OilChange']['MilesPE']
mpe_EVairFilter=maintPE.ix['EV_AirFilter']['MilesPE']
mpe_EVengineCoolant=maintPE.ix['EV_EngineCoolant']['MilesPE']
mpe_EVsparkPlugs=maintPE.ix['EV_SparkPlugs']['MilesPE']
mpe_EVbatteryCoolant=maintPE.ix['EV_BatteryCoolant']['MilesPE']
mpe_EVbatteryCheck=maintPE.ix['EV_BatteryCheck']['MilesPE']
mpe_PHEVbatteryCoolant=maintPE.ix['PHEV_BatteryCoolant']['MilesPE']
mpe_PHEVbatteryCheck=maintPE.ix['EV_BatteryCheck']['MilesPE']

maintenance=pd.DataFrame(calcsTemp['MY'])

def calcMaint(veh, case, system, cpe, mpe):
    b='$_Maintenance_'+system+'_'+veh+'_'+case
    maintenance[b]=(calcsTemp['vmt_'+veh+'_'+case]/mpe)*cpe*maintTechPensbyIndex[system+'_'+veh+'_'+case]
    return(maintenance)

#pass Car, ref case
calcMaint('car', 'ref', 'LRRT1', cpe_LRRT1, mpe_LRRT1)
calcMaint('car', 'ref', 'LRRT2', cpe_LRRT2, mpe_LRRT2)
calcMaint('car', 'ref', 'Diesel', cpe_Diesel, mpe_Diesel)
calcMaint('car', 'ref', 'EV_OilChange', cpe_EVoilChange, mpe_EVoilChange)
calcMaint('car', 'ref', 'EV_AirFilter', cpe_EVairFilter, mpe_EVairFilter)
calcMaint('car', 'ref', 'EV_EngineCoolant', cpe_EVengineCoolant, mpe_EVengineCoolant)
calcMaint('car', 'ref', 'EV_SparkPlugs', cpe_EVsparkPlugs, mpe_EVsparkPlugs)
calcMaint('car', 'ref', 'EV_BatteryCoolant', cpe_EVbatteryCoolant, mpe_EVbatteryCoolant)
calcMaint('car', 'ref', 'EV_BatteryCheck', cpe_EVbatteryCheck, mpe_EVbatteryCheck)
calcMaint('car', 'ref', 'PHEV_BatteryCoolant', cpe_PHEVbatteryCoolant, mpe_PHEVbatteryCoolant)
calcMaint('car', 'ref', 'PHEV_BatteryCheck', cpe_PHEVbatteryCheck, mpe_PHEVbatteryCheck)
#pass truck, ref case
calcMaint('truck', 'ref', 'LRRT1', cpe_LRRT1, mpe_LRRT1)
calcMaint('truck', 'ref', 'LRRT2', cpe_LRRT2, mpe_LRRT2)
calcMaint('truck', 'ref', 'Diesel', cpe_Diesel, mpe_Diesel)
calcMaint('truck', 'ref', 'EV_OilChange', cpe_EVoilChange, mpe_EVoilChange)
calcMaint('truck', 'ref', 'EV_AirFilter', cpe_EVairFilter, mpe_EVairFilter)
calcMaint('truck', 'ref', 'EV_EngineCoolant', cpe_EVengineCoolant, mpe_EVengineCoolant)
calcMaint('truck', 'ref', 'EV_SparkPlugs', cpe_EVsparkPlugs, mpe_EVsparkPlugs)
calcMaint('truck', 'ref', 'EV_BatteryCoolant', cpe_EVbatteryCoolant, mpe_EVbatteryCoolant)
calcMaint('truck', 'ref', 'EV_BatteryCheck', cpe_EVbatteryCheck, mpe_EVbatteryCheck)
calcMaint('truck', 'ref', 'PHEV_BatteryCoolant', cpe_PHEVbatteryCoolant, mpe_PHEVbatteryCoolant)
calcMaint('truck', 'ref', 'PHEV_BatteryCheck', cpe_PHEVbatteryCheck, mpe_PHEVbatteryCheck)
#pass Car, control case
calcMaint('car','control', 'LRRT1', cpe_LRRT1, mpe_LRRT1)
calcMaint('car','control', 'LRRT2', cpe_LRRT2, mpe_LRRT2)
calcMaint('car','control', 'Diesel', cpe_Diesel, mpe_Diesel)
calcMaint('car','control', 'EV_OilChange', cpe_EVoilChange, mpe_EVoilChange)
calcMaint('car','control', 'EV_AirFilter', cpe_EVairFilter, mpe_EVairFilter)
calcMaint('car','control', 'EV_EngineCoolant', cpe_EVengineCoolant, mpe_EVengineCoolant)
calcMaint('car','control', 'EV_SparkPlugs', cpe_EVsparkPlugs, mpe_EVsparkPlugs)
calcMaint('car','control', 'EV_BatteryCoolant', cpe_EVbatteryCoolant, mpe_EVbatteryCoolant)
calcMaint('car','control', 'EV_BatteryCheck', cpe_EVbatteryCheck, mpe_EVbatteryCheck)
calcMaint('car','control', 'PHEV_BatteryCoolant', cpe_PHEVbatteryCoolant, mpe_PHEVbatteryCoolant)
calcMaint('car','control', 'PHEV_BatteryCheck', cpe_PHEVbatteryCheck, mpe_PHEVbatteryCheck)
#pass truck, control case
calcMaint('truck','control', 'LRRT1', cpe_LRRT1, mpe_LRRT1)
calcMaint('truck','control', 'LRRT2', cpe_LRRT2, mpe_LRRT2)
calcMaint('truck','control', 'Diesel', cpe_Diesel, mpe_Diesel)
calcMaint('truck','control', 'EV_OilChange', cpe_EVoilChange, mpe_EVoilChange)
calcMaint('truck','control', 'EV_AirFilter', cpe_EVairFilter, mpe_EVairFilter)
calcMaint('truck','control', 'EV_EngineCoolant', cpe_EVengineCoolant, mpe_EVengineCoolant)
calcMaint('truck','control', 'EV_SparkPlugs', cpe_EVsparkPlugs, mpe_EVsparkPlugs)
calcMaint('truck','control', 'EV_BatteryCoolant', cpe_EVbatteryCoolant, mpe_EVbatteryCoolant)
calcMaint('truck','control', 'EV_BatteryCheck', cpe_EVbatteryCheck, mpe_EVbatteryCheck)
calcMaint('truck','control', 'PHEV_BatteryCoolant', cpe_PHEVbatteryCoolant, mpe_PHEVbatteryCoolant)
calcMaint('truck','control', 'PHEV_BatteryCheck', cpe_PHEVbatteryCheck, mpe_PHEVbatteryCheck)
#maintenance.drop(labels=['MY'], axis=1, inplace=True)
maintenance.fillna(0, inplace=True)
maintenance['$_Maintenance_car_ref']=maintenance[['$_Maintenance_LRRT1_car_ref', '$_Maintenance_LRRT2_car_ref', '$_Maintenance_Diesel_car_ref',
                                                  '$_Maintenance_EV_OilChange_car_ref', '$_Maintenance_EV_AirFilter_car_ref', '$_Maintenance_EV_EngineCoolant_car_ref',
                                                  '$_Maintenance_EV_SparkPlugs_car_ref', '$_Maintenance_EV_BatteryCoolant_car_ref', '$_Maintenance_EV_BatteryCheck_car_ref',
                                                  '$_Maintenance_PHEV_BatteryCoolant_car_ref', '$_Maintenance_PHEV_BatteryCheck_car_ref']].sum(axis=1)
maintenance['$_Maintenance_truck_ref']=maintenance[['$_Maintenance_LRRT1_truck_ref', '$_Maintenance_LRRT2_truck_ref', '$_Maintenance_Diesel_truck_ref',
                                                  '$_Maintenance_EV_OilChange_truck_ref', '$_Maintenance_EV_AirFilter_truck_ref', '$_Maintenance_EV_EngineCoolant_truck_ref',
                                                  '$_Maintenance_EV_SparkPlugs_truck_ref', '$_Maintenance_EV_BatteryCoolant_truck_ref', '$_Maintenance_EV_BatteryCheck_truck_ref',
                                                  '$_Maintenance_PHEV_BatteryCoolant_truck_ref', '$_Maintenance_PHEV_BatteryCheck_truck_ref']].sum(axis=1)
maintenance['$_Maintenance_fleet_ref']=maintenance[['$_Maintenance_car_ref', '$_Maintenance_truck_ref']].sum(axis=1)
maintenance['$_Maintenance_car_control']=maintenance[['$_Maintenance_LRRT1_car_control', '$_Maintenance_LRRT2_car_control', '$_Maintenance_Diesel_car_control',
                                                  '$_Maintenance_EV_OilChange_car_control', '$_Maintenance_EV_AirFilter_car_control', '$_Maintenance_EV_EngineCoolant_car_control',
                                                  '$_Maintenance_EV_SparkPlugs_car_control', '$_Maintenance_EV_BatteryCoolant_car_control', '$_Maintenance_EV_BatteryCheck_car_control',
                                                  '$_Maintenance_PHEV_BatteryCoolant_car_control', '$_Maintenance_PHEV_BatteryCheck_car_control']].sum(axis=1)
maintenance['$_Maintenance_truck_control']=maintenance[['$_Maintenance_LRRT1_truck_control', '$_Maintenance_LRRT2_truck_control', '$_Maintenance_Diesel_truck_control',
                                                  '$_Maintenance_EV_OilChange_truck_control', '$_Maintenance_EV_AirFilter_truck_control', '$_Maintenance_EV_EngineCoolant_truck_control',
                                                  '$_Maintenance_EV_SparkPlugs_truck_control', '$_Maintenance_EV_BatteryCoolant_truck_control', '$_Maintenance_EV_BatteryCheck_truck_control',
                                                  '$_Maintenance_PHEV_BatteryCoolant_truck_control', '$_Maintenance_PHEV_BatteryCheck_truck_control']].sum(axis=1)
maintenance['$_Maintenance_fleet_control']=maintenance[['$_Maintenance_car_control', '$_Maintenance_truck_control']].sum(axis=1)
maintenance['$_Maintenance_car_delta']=maintenance['$_Maintenance_car_control']-maintenance['$_Maintenance_car_ref']
maintenance['$_Maintenance_truck_delta']=maintenance['$_Maintenance_truck_control']-maintenance['$_Maintenance_truck_ref']
maintenance['$_Maintenance_fleet_delta']=maintenance['$_Maintenance_fleet_control']-maintenance['$_Maintenance_fleet_ref']
maintenanceMerge=pd.DataFrame(maintenance, columns=['$_Maintenance_car_ref', '$_Maintenance_truck_ref', '$_Maintenance_fleet_ref',
                                                    '$_Maintenance_car_control', '$_Maintenance_truck_control', '$_Maintenance_fleet_control',
                                                    '$_Maintenance_car_delta', '$_Maintenance_truck_delta', '$_Maintenance_fleet_delta'])
calcs=calcs.join(maintenanceMerge, how='left')

# pull together monetized values for easier summaries
#nonGHG benefits
calcs['$_0.03_low_nonGHG']=calcs[['$_0.03_low_NOx', '$_0.03_low_PM25', '$_0.03_low_SOx']].sum(axis=1)
calcs['$_0.03_mid_nonGHG']=calcs[['$_0.03_mid_NOx', '$_0.03_mid_PM25', '$_0.03_mid_SOx']].sum(axis=1)
calcs['$_0.03_high_nonGHG']=calcs[['$_0.03_high_NOx', '$_0.03_high_PM25', '$_0.03_high_SOx']].sum(axis=1)
calcs['$_0.07_low_nonGHG']=calcs[['$_0.07_low_NOx', '$_0.07_low_PM25', '$_0.07_low_SOx']].sum(axis=1)
calcs['$_0.07_mid_nonGHG']=calcs[['$_0.07_mid_NOx', '$_0.07_mid_PM25', '$_0.07_mid_SOx']].sum(axis=1)
calcs['$_0.07_high_nonGHG']=calcs[['$_0.07_high_NOx', '$_0.07_high_PM25', '$_0.07_high_SOx']].sum(axis=1)

# benefits excluding fuel, and using the 3/7% mid nonGHG values; Note these cannot be passed thru the discount function because SCC is discounted differently than other metrics
calcs['$_Benefits_0.03_mid_SCGHG_5']=calcs[['$_0.03_mid_nonGHG', '$_SCGHG_5', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
calcs['$_Benefits_0.03_mid_SCGHG_3']=calcs[['$_0.03_mid_nonGHG', '$_SCGHG_3', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
calcs['$_Benefits_0.03_mid_SCGHG_25']=calcs[['$_0.03_mid_nonGHG', '$_SCGHG_25', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
calcs['$_Benefits_0.03_mid_SCGHG_395']=calcs[['$_0.03_mid_nonGHG', '$_SCGHG_395', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
calcs['$_Benefits_0.07_mid_SCGHG_5']=calcs[['$_0.07_mid_nonGHG', '$_SCGHG_5', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
calcs['$_Benefits_0.07_mid_SCGHG_3']=calcs[['$_0.07_mid_nonGHG', '$_SCGHG_3', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
calcs['$_Benefits_0.07_mid_SCGHG_25']=calcs[['$_0.07_mid_nonGHG', '$_SCGHG_25', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
calcs['$_Benefits_0.07_mid_SCGHG_395']=calcs[['$_0.07_mid_nonGHG', '$_SCGHG_395', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)

# net benefits, using the 3/7% mid nonGHG values; Note that these cannot be passed thru the discount function
temp_net=calcs[['$_pretax_petrol_fleet_delta', '$_Tech_fleet_delta', '$_Maintenance_fleet_delta']].sum(axis=1)
calcs['$_NetBenefits_0.03_mid_SCGHG_5']=calcs['$_Benefits_0.03_mid_SCGHG_5'] + temp_net
calcs['$_NetBenefits_0.03_mid_SCGHG_3']=calcs['$_Benefits_0.03_mid_SCGHG_3'] + temp_net
calcs['$_NetBenefits_0.03_mid_SCGHG_25']=calcs['$_Benefits_0.03_mid_SCGHG_25'] + temp_net
calcs['$_NetBenefits_0.03_mid_SCGHG_395']=calcs['$_Benefits_0.03_mid_SCGHG_395'] + temp_net
calcs['$_NetBenefits_0.07_mid_SCGHG_5']=calcs['$_Benefits_0.07_mid_SCGHG_5'] + temp_net
calcs['$_NetBenefits_0.07_mid_SCGHG_3']=calcs['$_Benefits_0.07_mid_SCGHG_3'] + temp_net
calcs['$_NetBenefits_0.07_mid_SCGHG_25']=calcs['$_Benefits_0.07_mid_SCGHG_25'] + temp_net
calcs['$_NetBenefits_0.07_mid_SCGHG_395']=calcs['$_Benefits_0.07_mid_SCGHG_395'] + temp_net

##########################################################
#######        now work on payback calcs      ############
##########################################################
survival_2 = inputFile.parse('survival')
payback=pd.DataFrame(calcs, columns=['MY', 'age', 'CY'])
payback=payback.merge(Insurance, on='age', how='left', sort=False)
payback=payback.merge(survival_2, on='age', how='left', sort=False)
payback=payback.set_index(payback.MY.map(str)+'-'+payback.age.map(str))

payback['population_car_CAFE']=calcs['sales_car_CAFE']*payback['surv_car']
payback['population_truck_CAFE']=calcs['sales_truck_CAFE']*payback['surv_truck']
payback['population_fleet']=payback[['population_car_CAFE', 'population_truck_CAFE']].sum(axis=1)
payback['population_car_MOVES']=calcs['sales_car_MOVES']*payback['surv_car']
payback['population_trucktocar']=calcs['sales_trucktocar']*payback['surv_truck']
payback['population_truck_MOVES']=calcs['sales_truck_MOVES']*payback['surv_truck']

# first, maintenance per vehicle per year
payback['$_perVeh_Maintenance_car_delta']=calcs['$_Maintenance_car_delta']/payback['population_car_CAFE']
payback['$_perVeh_Maintenance_truck_delta']=calcs['$_Maintenance_truck_delta']/payback['population_truck_CAFE']
payback['$_perVeh_Maintenance_fleet_delta']=calcs['$_Maintenance_fleet_delta']/payback['population_fleet']

# for fuel costs, need to use ref case VMT since payback does not consider rebound miles; so here
# we use $_perMi in the ref & control cases and VMT in the ref case
def calc_FuelPerVehPerYear(veh):
    b = '$_perVeh_Fuel_' + veh + '_delta'
    c = '$_perMi_retail_' + veh + '_ref'
    d = '$_perMi_retail_' + veh + '_control'
    e = 'miPerYear_' + veh + '_ref'
    payback[b] = (calcs[d]-calcs[c])*calcs[e]
    return(payback)

calc_FuelPerVehPerYear('car')
calc_FuelPerVehPerYear('truck')
payback['$_perVeh_Fuel_fleet_delta']=(payback['population_car_CAFE']*payback['$_perVeh_Fuel_car_delta']+
            payback['population_truck_CAFE']*payback['$_perVeh_Fuel_truck_delta'])/payback['population_fleet']

# pull in tech costs per veh deltas
payback['$_perVeh_Tech_car_delta']=calcs['$_perVeh_Tech_car_delta']
payback['$_perVeh_Tech_truck_delta']=calcs['$_perVeh_Tech_truck_delta']
payback['$_perVeh_Tech_fleet_delta']=calcs['$_perVeh_Tech_fleet_delta']

# taxes are only on the original purchase, not each year
def get_taxes(veh, tax):
    b='$_perVeh_Taxes_'+veh+'_delta'
    c='$_perVeh_Tech_'+veh+'_delta'
    payback[b]=payback[c]*tax
    return(payback)

get_taxes('car', taxes)
get_taxes('truck', taxes)
get_taxes('fleet', taxes)

# payback_temp provides tech costs every year that can be multiplied by insurance costs to get increased insurance costs each year
# this needs to be sorted by MY then age to do the ffill
payback_temp=pd.DataFrame(calcs, columns=['MY', 'age', 'CY', 'sales_car_CAFE', 'sales_truck_CAFE'])
payback_temp=payback_temp.join(vehcosts)
payback_temp=payback_temp.sort_values(by=['MY', 'age'], ascending=[True, True])
payback_temp['$_perVeh_Tech_car_delta']=payback_temp['$_perVeh_Tech_car_control']-payback_temp['$_perVeh_Tech_car_ref']
payback_temp['$_perVeh_Tech_truck_delta']=payback_temp['$_perVeh_Tech_truck_control']-payback_temp['$_perVeh_Tech_truck_ref']
payback_temp['$_perVeh_Tech_fleet_delta']=(payback_temp['$_perVeh_Tech_car_delta']*payback_temp['sales_car_CAFE']+
    payback_temp['$_perVeh_Tech_truck_delta']*payback_temp['sales_truck_CAFE'])/payback_temp[['sales_car_CAFE','sales_truck_CAFE']].sum(axis=1)
#payback_temp.fillna(method='ffill', inplace=True)
#payback_temp.fillna(0, inplace=True)

# now resort the payback dataframe to be consistent with payback_temp
payback=payback.sort_values(by=['MY', 'age'], ascending=[True, True])

# now that payback is structured like payback_temp, add taxes to payback_temp vehicle costs for use in loan amount calcs
payback_temp['$_LoanAmount_car'] = payback_temp['$_perVeh_Tech_car_delta'] + payback['$_perVeh_Taxes_car_delta']
payback_temp['$_LoanAmount_truck'] = payback_temp['$_perVeh_Tech_truck_delta'] + payback['$_perVeh_Taxes_truck_delta']
payback_temp['$_LoanAmount_fleet'] = payback_temp['$_perVeh_Tech_fleet_delta'] + payback['$_perVeh_Taxes_fleet_delta']
payback_temp.fillna(method='ffill', inplace=True)
payback_temp.fillna(0, inplace=True)

# increased insurance costs; new on 20160924 - include survival rates
def insurance(veh):
    b='$_perVeh_Insurance_'+veh+'_delta'
    c='$_perVeh_Tech_'+veh+'_delta'
    d='surv_' + veh
    payback[b]=payback_temp[c]*payback['C&C_asPctofBasePrice']*payback[d]
    return(payback)

insurance('car')
insurance('truck')
#insurance('fleet')
payback['$_perVeh_Insurance_fleet_delta'] = (payback['$_perVeh_Insurance_car_delta'] * payback['population_car_CAFE'] \
                                + payback['$_perVeh_Insurance_truck_delta'] * payback['population_truck_CAFE'])/payback['population_fleet']

# now generate the payments on loan purchases which would include both the increased tech and taxes
#payback_loan = pd.DataFrame(payback, columns=['MY', 'age', 'CY'])
#payback_loan['$_LoanAmount_car'] = payback[['$_perVeh_Tech_car_delta', '$_perVeh_Tech_Taxes_car_delta']].sum(axis=1)
#payback_loan['$_LoanAmount_truck'] = payback[['$_perVeh_Tech_truck_delta', '$_perVeh_Tech_Taxes_truck_delta']].sum(axis=1)
#payback_loan['$_LoanAmount_fleet'] = payback[['$_perVeh_Tech_fleet_delta', '$_perVeh_Tech_Taxes_fleet_delta']].sum(axis=1)

#payment = pd.DataFrame(payback_temp, columns=['age', '$_LoanAmount_car', '$_LoanAmount_truck', '$_LoanAmount_fleet'])
payment = pd.DataFrame(payback_temp, columns=['age', '$_LoanAmount_car', '$_LoanAmount_truck'])

def loan(rate, period):
    # the payment formula used is PMT=PV*(rate/(1-((1+rate)^-period))); the PV changes by veh but the other portion is constant for each rate & period
    # this also uses the effective interest rate rather than than the simple annual rate
    # this also includes survival to place payment terms on same basis as fuel savings, maintenance, etc.
    constant = (rate/12)/(1-((1+rate/12)**(-period*12))) * 12
    payment_temp = payment.copy()
    payment_temp.loc[(payment_temp['age'] < period), 'car'] = payment_temp['$_LoanAmount_car'] * constant
    payment_temp.loc[(payment_temp['age'] < period), 'truck'] = payment_temp['$_LoanAmount_truck'] * constant
    #payment_temp.loc[(payment_temp['age'] < period), 'fleet'] = payment_temp['$_LoanAmount_fleet'] * constant
    #payment_temp.drop(payment.columns, axis=1, inplace=True)
    payment_temp.fillna(0, inplace=True)
    payback['$_LoanPayment_car_' + str(period) + 'yr'] = payment_temp['car'] * payback['surv_car']
    payback['$_LoanPayment_truck_' + str(period) + 'yr'] = payment_temp['truck'] * payback['surv_truck']
    payback['$_LoanPayment_fleet_' + str(period) + 'yr'] = (payback['$_LoanPayment_car_' + str(period) + 'yr'] * payback['population_car_CAFE'] \
                        + payback['$_LoanPayment_truck_' + str(period) + 'yr'] * payback['population_truck_CAFE']) / payback['population_fleet']
    #payback['$_LoanPayment_fleet_' + str(period) + 'yr'] = payment_temp['fleet']
    return(payback)

loan(loanrate4, 4)
loan(loanrate5, 5)
loan(loanrate6, 6)

# now add tech costs, insurance & taxes to get the delta purchase costs for a cash purchase
payback['$_perVeh_Cash_car_delta']=payback[['$_perVeh_Tech_car_delta','$_perVeh_Insurance_car_delta','$_perVeh_Taxes_car_delta']].sum(axis=1)
payback['$_perVeh_Cash_truck_delta']=payback[['$_perVeh_Tech_truck_delta','$_perVeh_Insurance_truck_delta','$_perVeh_Taxes_truck_delta']].sum(axis=1)
payback['$_perVeh_Cash_fleet_delta']=payback[['$_perVeh_Tech_fleet_delta','$_perVeh_Insurance_fleet_delta','$_perVeh_Taxes_fleet_delta']].sum(axis=1)

# now add loan payment (which includes taxes) & insurance to get the delta purchase costs for a 4 year loan on a new vehicle purchase
payback['$_perVeh_Loan4yr_car_delta'] = payback[['$_LoanPayment_car_4yr', '$_perVeh_Insurance_car_delta']].sum(axis=1)
payback['$_perVeh_Loan4yr_truck_delta'] = payback[['$_LoanPayment_truck_4yr', '$_perVeh_Insurance_truck_delta']].sum(axis=1)
payback['$_perVeh_Loan4yr_fleet_delta'] = payback[['$_LoanPayment_fleet_4yr', '$_perVeh_Insurance_fleet_delta']].sum(axis=1)

# now add loan payment (which includes taxes) & insurance to get the delta purchase costs for a 5 year loan on a new vehicle purchase
payback['$_perVeh_Loan5yr_car_delta'] = payback[['$_LoanPayment_car_5yr', '$_perVeh_Insurance_car_delta']].sum(axis=1)
payback['$_perVeh_Loan5yr_truck_delta'] = payback[['$_LoanPayment_truck_5yr', '$_perVeh_Insurance_truck_delta']].sum(axis=1)
payback['$_perVeh_Loan5yr_fleet_delta'] = payback[['$_LoanPayment_fleet_5yr', '$_perVeh_Insurance_fleet_delta']].sum(axis=1)

# now add loan payment (which includes taxes) & insurance to get the delta purchase costs for a 6 year loan on a new vehicle purchase
payback['$_perVeh_Loan6yr_car_delta'] = payback[['$_LoanPayment_car_6yr', '$_perVeh_Insurance_car_delta']].sum(axis=1)
payback['$_perVeh_Loan6yr_truck_delta'] = payback[['$_LoanPayment_truck_6yr', '$_perVeh_Insurance_truck_delta']].sum(axis=1)
payback['$_perVeh_Loan6yr_fleet_delta'] = payback[['$_LoanPayment_fleet_6yr', '$_perVeh_Insurance_fleet_delta']].sum(axis=1)


###################################################
#OUTPUTS
###################################################
#Identify those variables that are static and can't be summed.
nonadditivevars=['MY', 'age', 'CY', 'sales_car_CAFE', 'sales_truck_CAFE', 'sales_fleet',
                 'sales_car_MOVES/sales_car_CAFE', 'sales_car_MOVES', 'sales_truck_MOVES', 'sales_trucktocar', 'gpm_CO2_road_fleet_ref', 'gpm_CO2_road_car_ref',
                 'gpm_CO2_road_truck_ref', 'gpm_CO2_road_fleet_control', 'gpm_CO2_road_car_control', 'gpm_CO2_road_truck_control',
                 'kWhPerMi_road_fleet_ref', 'kWhPerMi_road_car_ref', 'kWhPerMi_road_truck_ref',
                 'kWhPerMi_road_fleet_control', 'kWhPerMi_road_car_control', 'kWhPerMi_road_truck_control',
                 '$_perGal_retail_diesel_'+case_ref, '$_perGal_pretax_diesel_'+case_ref, '$_perGal_retail_gasoline_'+case_ref,'$_perGal_pretax_gasoline_'+case_ref,
#                 '$_perGal_retail_diesel_aeoL','$_perGal_pretax_diesel_aeoL','$_perGal_retail_gasoline_aeoL','$_perGal_pretax_gasoline_aeoL',
#                 '$_perGal_retail_diesel_aeoH','$_perGal_pretax_diesel_aeoH','$_perGal_retail_gasoline_aeoH','$_perGal_pretax_gasoline_aeoH',
#                 'price_CNG_retail_ref','price_CNG_untaxed_ref','price_CNG_retail_low','price_CNG_untaxed_low','price_CNG_retail_high','price_CNG_untaxed_high',
                 '$_perkWh_retail_elec','$_perkWh_pretax_elec','mpg_retail_car_ref', 'mpg_retail_truck_ref', 'mpg_retail_car_control',
                 'mpg_retail_truck_control','$_perMi_retail_gasoline_car_ref','$_perMi_retail_gasoline_truck_ref','$_perMi_retail_gasoline_car_control',
                 '$_perMi_retail_gasoline_truck_control','$_perMi_retail_elec_car_ref','$_perMi_retail_elec_truck_ref','$_perMi_retail_elec_car_control','$_perMi_retail_elec_truck_control',
                 '$_perMi_retail_car_ref','$_perMi_retail_truck_ref','$_perMi_retail_car_control','$_perMi_retail_truck_control']

#Identify those variables that are monetized, with the exception of SCC variable which are handled below (because of the unique discounting of monetized SCC)
monetized_vars_3=['$_retail_petrol_car_delta','$_retail_petrol_truck_delta','$_retail_petrol_fleet_delta','$_pretax_petrol_car_delta','$_pretax_petrol_truck_delta','$_pretax_petrol_fleet_delta',
                '$_retail_elec_car_delta', '$_retail_elec_truck_delta', '$_retail_elec_fleet_delta',
                '$_EnergySecurity','$_Refueling_car','$_Refueling_truck','$_Refueling_fleet','$_retail_petrol_fleet_rebound','$_TravelBenefit_car','$_TravelBenefit_truck','$_TravelBenefit_fleet',
                '$_congestion_car','$_congestion_truck','$_crashes_car','$_crashes_truck','$_noise_car','$_noise_truck',
                '$_0.03_low_up_PM25','$_0.03_mid_up_PM25','$_0.03_high_up_PM25',
                '$_0.03_low_up_NOx','$_0.03_mid_up_NOx','$_0.03_high_up_NOx',
                '$_0.03_low_up_SOx','$_0.03_mid_up_SOx','$_0.03_high_up_SOx',
                '$_0.03_low_down_PM25','$_0.03_mid_down_PM25','$_0.03_high_down_PM25',
                '$_0.03_low_down_NOx','$_0.03_mid_down_NOx','$_0.03_high_down_NOx',
                '$_0.03_low_down_SOx','$_0.03_mid_down_SOx','$_0.03_high_down_SOx',
                '$_0.03_low_PM25','$_0.03_mid_PM25','$_0.03_high_PM25',
                '$_0.03_low_NOx','$_0.03_mid_NOx','$_0.03_high_NOx',
                '$_0.03_low_SOx','$_0.03_mid_SOx','$_0.03_high_SOx',
                '$_NCCP_car','$_NCCP_truck','$_NCCP_fleet',
                '$_Tech_car_ref','$_Tech_car_control','$_Tech_truck_ref','$_Tech_truck_control',
                '$_Tech_car_delta','$_Tech_truck_delta','$_Tech_fleet_delta',
                '$_Maintenance_car_ref','$_Maintenance_car_control','$_Maintenance_truck_ref','$_Maintenance_truck_control',
                '$_Maintenance_car_delta','$_Maintenance_truck_delta','$_Maintenance_fleet_delta',
                '$_0.03_low_nonGHG','$_0.03_mid_nonGHG','$_0.03_high_nonGHG']

monetized_vars_7=['$_retail_petrol_car_delta','$_retail_petrol_truck_delta','$_retail_petrol_fleet_delta','$_pretax_petrol_car_delta','$_pretax_petrol_truck_delta','$_pretax_petrol_fleet_delta',
                '$_retail_elec_car_delta', '$_retail_elec_truck_delta', '$_retail_elec_fleet_delta',
                '$_EnergySecurity','$_Refueling_car','$_Refueling_truck','$_Refueling_fleet','$_retail_petrol_fleet_rebound','$_TravelBenefit_car','$_TravelBenefit_truck','$_TravelBenefit_fleet',
                '$_congestion_car','$_congestion_truck','$_crashes_car','$_crashes_truck','$_noise_car','$_noise_truck',
                '$_0.07_low_up_PM25','$_0.07_mid_up_PM25','$_0.07_high_up_PM25',
                '$_0.07_low_up_NOx','$_0.07_mid_up_NOx','$_0.07_high_up_NOx',
                '$_0.07_low_up_SOx','$_0.07_mid_up_SOx','$_0.07_high_up_SOx',
                '$_0.07_low_down_PM25','$_0.07_mid_down_PM25','$_0.07_high_down_PM25',
                '$_0.07_low_down_NOx','$_0.07_mid_down_NOx','$_0.07_high_down_NOx',
                '$_0.07_low_down_SOx','$_0.07_mid_down_SOx','$_0.07_high_down_SOx',
                '$_0.07_low_PM25','$_0.07_mid_PM25','$_0.07_high_PM25',
                '$_0.07_low_NOx','$_0.07_mid_NOx','$_0.07_high_NOx',
                '$_0.07_low_SOx','$_0.07_mid_SOx','$_0.07_high_SOx',
                '$_NCCP_car','$_NCCP_truck','$_NCCP_fleet',
                '$_Tech_car_ref','$_Tech_car_control','$_Tech_truck_ref','$_Tech_truck_control',
                '$_Tech_car_delta','$_Tech_truck_delta','$_Tech_fleet_delta',
                '$_Maintenance_car_ref','$_Maintenance_car_control','$_Maintenance_truck_ref','$_Maintenance_truck_control',
                '$_Maintenance_car_delta','$_Maintenance_truck_delta','$_Maintenance_fleet_delta',
                '$_0.07_low_nonGHG','$_0.07_mid_nonGHG','$_0.07_high_nonGHG']

#Identify SCC monetized variables
monetized_SCC5=['$_SCCO2_5','$_SCCH4_5','$_SCN2O_5', '$_SCGHG_5']
monetized_SCC3=['$_SCCO2_3','$_SCCH4_3','$_SCN2O_3','$_SCGHG_3']
monetized_SCC25=['$_SCCO2_25','$_SCCH4_25','$_SCN2O_25', '$_SCGHG_25']
monetized_SCC395=['$_SCCO2_395','$_SCCH4_395','$_SCN2O_395', '$_SCGHG_395']

#identify payback variables
#note that Cumumlative operating savings cannot be discounted, we first need to discount all streams, then do a
#cumulative sum of the discounted streams
payback_vars=['$_perVeh_Maintenance_car_delta','$_perVeh_Tech_car_delta','$_perVeh_Insurance_car_delta','$_perVeh_Taxes_car_delta',
                '$_perVeh_Cash_car_delta','$_perVeh_Loan4yr_car_delta','$_perVeh_Loan5yr_car_delta','$_perVeh_Loan6yr_car_delta',
                '$_perVeh_Fuel_car_delta','$_perVeh_Maintenance_truck_delta','$_perVeh_Tech_truck_delta','$_perVeh_Insurance_truck_delta',
                '$_perVeh_Taxes_truck_delta','$_perVeh_Cash_truck_delta','$_perVeh_Loan4yr_truck_delta','$_perVeh_Loan5yr_truck_delta',
                '$_perVeh_Loan6yr_truck_delta','$_perVeh_Fuel_truck_delta','$_perVeh_Maintenance_fleet_delta','$_perVeh_Tech_fleet_delta',
                '$_perVeh_Insurance_fleet_delta','$_perVeh_Taxes_fleet_delta','$_perVeh_Cash_fleet_delta','$_perVeh_Loan4yr_fleet_delta',
                '$_perVeh_Loan5yr_fleet_delta','$_perVeh_Loan6yr_fleet_delta','$_perVeh_Fuel_fleet_delta']

summarycalcs=calcs.drop(nonadditivevars, axis=1)

#Create Summary by MY
MYSums=summarycalcs.groupby(pd.Series(calcs['MY'],name='MY_Key')).sum()
MYMeans=calcs[nonadditivevars].groupby(pd.Series(calcs['MY'],name='MY_Key')).mean()
MYSummary=MYMeans.join(MYSums)
MYSummary.drop(labels=['age', 'CY'], axis=1, inplace=True)

#Create Summary by CY
CYSums=summarycalcs.groupby(pd.Series(calcs['CY'],name='CY_Key')).sum()
CYMeans=calcs[nonadditivevars].groupby(pd.Series(calcs['CY'],name='CY_Key')).mean()
CYSummary=CYMeans.join(CYSums)
CYSummary.drop(labels=['MY', 'age'], axis=1, inplace=True)

#Create summary inputs
Input_sum=pd.DataFrame(pd.Series(['PythonScript version','InputSheet version','case_ref','case_control','scenario','USD_Basis'], name='Inputs'))
Input_sum=Input_sum.join(pd.DataFrame(pd.Series([scriptversion, version, case_ref, case_control, scenario, USD_Basis], name='Selection')))

def discount(col_source, rate, source, destination):
    for i in range(0,len(col_source)):
        col=col_source[i]
        destination[col]=source[col]/((1+rate)**(destination['CY']-discountYr+.5))
    return(destination)

# first, discount MY values using calcs as the source
MY_3pct_all=pd.DataFrame(calcs, columns=['MY', 'age', 'CY'])
MY_7pct_all=pd.DataFrame(calcs, columns=['MY', 'age', 'CY'])
MY_5pct_SCC_all=pd.DataFrame(calcs, columns=['MY', 'age', 'CY'])
MY_3pct_SCC_all=pd.DataFrame(calcs, columns=['MY', 'age', 'CY'])
MY_2_5pct_SCC_all=pd.DataFrame(calcs, columns=['MY', 'age', 'CY'])
MY_395pct_SCC_all=pd.DataFrame(calcs, columns=['MY', 'age', 'CY'])
discount(monetized_vars_3, 0.03, calcs, MY_3pct_all)
discount(monetized_vars_7, 0.07, calcs, MY_7pct_all)
discount(monetized_SCC5, 0.05, calcs, MY_5pct_SCC_all)
discount(monetized_SCC3, 0.03, calcs, MY_3pct_SCC_all)
discount(monetized_SCC25, 0.025, calcs, MY_2_5pct_SCC_all)
discount(monetized_SCC395, 0.03, calcs, MY_395pct_SCC_all)

MY_3pct=MY_3pct_all.groupby(pd.Series(MY_3pct_all['MY'], name='MY_Key')).sum()
MY_7pct=MY_7pct_all.groupby(pd.Series(MY_7pct_all['MY'], name='MY_Key')).sum()
MY_5pct_SCC=MY_5pct_SCC_all.groupby(pd.Series(MY_5pct_SCC_all['MY'], name='MY_Key')).sum()
MY_3pct_SCC=MY_3pct_SCC_all.groupby(pd.Series(MY_3pct_SCC_all['MY'], name='MY_Key')).sum()
MY_2_5pct_SCC=MY_2_5pct_SCC_all.groupby(pd.Series(MY_2_5pct_SCC_all['MY'], name='MY_Key')).sum()
MY_395pct_SCC=MY_395pct_SCC_all.groupby(pd.Series(MY_395pct_SCC_all['MY'], name='MY_Key')).sum()
MY_3pct.drop(labels=['MY', 'age', 'CY'], axis=1, inplace=True)
MY_7pct.drop(labels=['MY', 'age', 'CY'], axis=1, inplace=True)
MY_5pct_SCC.drop(labels=['MY', 'age', 'CY'], axis=1, inplace=True)
MY_3pct_SCC.drop(labels=['MY', 'age', 'CY'], axis=1, inplace=True)
MY_2_5pct_SCC.drop(labels=['MY', 'age', 'CY'], axis=1, inplace=True)
MY_395pct_SCC.drop(labels=['MY', 'age', 'CY'], axis=1, inplace=True)
MY_SCC=MY_5pct_SCC.join(MY_3pct_SCC)
MY_SCC=MY_SCC.join(MY_2_5pct_SCC)
MY_SCC=MY_SCC.join(MY_395pct_SCC)

# second, discount CY values using calcs as the source
CY_3pct_all=pd.DataFrame(calcs, columns=['CY'])
CY_7pct_all=pd.DataFrame(calcs, columns=['CY'])
CY_5pct_SCC_all=pd.DataFrame(calcs, columns=['CY'])
CY_3pct_SCC_all=pd.DataFrame(calcs, columns=['CY'])
CY_2_5pct_SCC_all=pd.DataFrame(calcs, columns=['CY'])
CY_395pct_SCC_all=pd.DataFrame(calcs, columns=['CY'])
discount(monetized_vars_3, 0.03, calcs, CY_3pct_all)
discount(monetized_vars_7, 0.07, calcs, CY_7pct_all)
discount(monetized_SCC5, 0.05, calcs, CY_5pct_SCC_all)
discount(monetized_SCC3, 0.03, calcs, CY_3pct_SCC_all)
discount(monetized_SCC25, 0.025, calcs, CY_2_5pct_SCC_all)
discount(monetized_SCC395, 0.03, calcs, CY_395pct_SCC_all)

CY_3pct=CY_3pct_all.groupby(pd.Series(CY_3pct_all['CY'], name='CY_Key')).sum()
CY_7pct=CY_7pct_all.groupby(pd.Series(CY_7pct_all['CY'], name='CY_Key')).sum()
CY_5pct_SCC=CY_5pct_SCC_all.groupby(pd.Series(CY_5pct_SCC_all['CY'], name='CY_Key')).sum()
CY_3pct_SCC=CY_3pct_SCC_all.groupby(pd.Series(CY_3pct_SCC_all['CY'], name='CY_Key')).sum()
CY_2_5pct_SCC=CY_2_5pct_SCC_all.groupby(pd.Series(CY_2_5pct_SCC_all['CY'], name='CY_Key')).sum()
CY_395pct_SCC=CY_395pct_SCC_all.groupby(pd.Series(CY_395pct_SCC_all['CY'], name='CY_Key')).sum()
CY_3pct.drop(labels=['CY'], axis=1, inplace=True)
CY_7pct.drop(labels=['CY'], axis=1, inplace=True)
CY_5pct_SCC.drop(labels=['CY'], axis=1, inplace=True)
CY_3pct_SCC.drop(labels=['CY'], axis=1, inplace=True)
CY_2_5pct_SCC.drop(labels=['CY'], axis=1, inplace=True)
CY_395pct_SCC.drop(labels=['CY'], axis=1, inplace=True)
CY_SCC=CY_5pct_SCC.join(CY_3pct_SCC)
CY_SCC=CY_SCC.join(CY_2_5pct_SCC)
CY_SCC=CY_SCC.join(CY_395pct_SCC)

# pull together monetized values for easier summaries
# First for MY analysis; benefits excluding fuel, and using the 3/7% mid nonGHG values
MY_3pct_temp=MY_3pct[['$_0.03_mid_nonGHG', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
MY_3pct['$_Benefits_0.03_mid_SCGHG_5']=MY_SCC['$_SCGHG_5'] + MY_3pct_temp
MY_3pct['$_Benefits_0.03_mid_SCGHG_3']=MY_SCC['$_SCGHG_3'] + MY_3pct_temp
MY_3pct['$_Benefits_0.03_mid_SCGHG_25']=MY_SCC['$_SCGHG_25'] + MY_3pct_temp
MY_3pct['$_Benefits_0.03_mid_SCGHG_395']=MY_SCC['$_SCGHG_395'] + MY_3pct_temp
MY_7pct_temp=MY_7pct[['$_0.07_mid_nonGHG', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
MY_7pct['$_Benefits_0.07_mid_SCGHG_5']=MY_SCC['$_SCGHG_5'] + MY_7pct_temp
MY_7pct['$_Benefits_0.07_mid_SCGHG_3']=MY_SCC['$_SCGHG_3'] + MY_7pct_temp
MY_7pct['$_Benefits_0.07_mid_SCGHG_25']=MY_SCC['$_SCGHG_25'] + MY_7pct_temp
MY_7pct['$_Benefits_0.07_mid_SCGHG_395']=MY_SCC['$_SCGHG_395'] + MY_7pct_temp

# MY analysis net benefits, using the 3/7% mid nonGHG values
# NOTE 20160816: should these net benefits include electricity expenditure deltas in addition to petrol?
MY_3pct_temp=MY_3pct[['$_pretax_petrol_fleet_delta', '$_Tech_fleet_delta', '$_Maintenance_fleet_delta']].sum(axis=1)
MY_3pct['$_NetBenefits_0.03_mid_SCGHG_5']=MY_3pct['$_Benefits_0.03_mid_SCGHG_5'] + MY_3pct_temp
MY_3pct['$_NetBenefits_0.03_mid_SCGHG_3']=MY_3pct['$_Benefits_0.03_mid_SCGHG_3'] + MY_3pct_temp
MY_3pct['$_NetBenefits_0.03_mid_SCGHG_25']=MY_3pct['$_Benefits_0.03_mid_SCGHG_25'] + MY_3pct_temp
MY_3pct['$_NetBenefits_0.03_mid_SCGHG_395']=MY_3pct['$_Benefits_0.03_mid_SCGHG_395'] + MY_3pct_temp
MY_7pct_temp=MY_7pct[['$_pretax_petrol_fleet_delta', '$_Tech_fleet_delta', '$_Maintenance_fleet_delta']].sum(axis=1)
MY_7pct['$_NetBenefits_0.07_mid_SCGHG_5']=MY_7pct['$_Benefits_0.07_mid_SCGHG_5'] + MY_7pct_temp
MY_7pct['$_NetBenefits_0.07_mid_SCGHG_3']=MY_7pct['$_Benefits_0.07_mid_SCGHG_3'] + MY_7pct_temp
MY_7pct['$_NetBenefits_0.07_mid_SCGHG_25']=MY_7pct['$_Benefits_0.07_mid_SCGHG_25'] + MY_7pct_temp
MY_7pct['$_NetBenefits_0.07_mid_SCGHG_395']=MY_7pct['$_Benefits_0.07_mid_SCGHG_395'] + MY_7pct_temp

# CY analysis; first benefits
CY_3pct_temp=CY_3pct[['$_0.03_mid_nonGHG', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
CY_3pct['$_Benefits_0.03_mid_SCGHG_5']=CY_SCC['$_SCGHG_5'] + CY_3pct_temp
CY_3pct['$_Benefits_0.03_mid_SCGHG_3']=CY_SCC['$_SCGHG_3'] + CY_3pct_temp
CY_3pct['$_Benefits_0.03_mid_SCGHG_25']=CY_SCC['$_SCGHG_25'] + CY_3pct_temp
CY_3pct['$_Benefits_0.03_mid_SCGHG_395']=CY_SCC['$_SCGHG_395'] + CY_3pct_temp
CY_7pct_temp=CY_7pct[['$_0.07_mid_nonGHG', '$_NCCP_fleet', '$_EnergySecurity', '$_TravelBenefit_fleet', '$_Refueling_fleet']].sum(axis=1)
CY_7pct['$_Benefits_0.07_mid_SCGHG_5']=CY_SCC['$_SCGHG_5'] + CY_7pct_temp
CY_7pct['$_Benefits_0.07_mid_SCGHG_3']=CY_SCC['$_SCGHG_3'] + CY_7pct_temp
CY_7pct['$_Benefits_0.07_mid_SCGHG_25']=CY_SCC['$_SCGHG_25'] + CY_7pct_temp
CY_7pct['$_Benefits_0.07_mid_SCGHG_395']=CY_SCC['$_SCGHG_395'] + CY_7pct_temp

# CY analysis net benefits
CY_3pct_temp=CY_3pct[['$_pretax_petrol_fleet_delta', '$_Tech_fleet_delta', '$_Maintenance_fleet_delta']].sum(axis=1)
CY_3pct['$_NetBenefits_0.03_mid_SCGHG_5']=CY_3pct['$_Benefits_0.03_mid_SCGHG_5'] + CY_3pct_temp
CY_3pct['$_NetBenefits_0.03_mid_SCGHG_3']=CY_3pct['$_Benefits_0.03_mid_SCGHG_3'] + CY_3pct_temp
CY_3pct['$_NetBenefits_0.03_mid_SCGHG_25']=CY_3pct['$_Benefits_0.03_mid_SCGHG_25'] + CY_3pct_temp
CY_3pct['$_NetBenefits_0.03_mid_SCGHG_395']=CY_3pct['$_Benefits_0.03_mid_SCGHG_395'] + CY_3pct_temp
CY_7pct_temp=CY_7pct[['$_pretax_petrol_fleet_delta', '$_Tech_fleet_delta', '$_Maintenance_fleet_delta']].sum(axis=1)
CY_7pct['$_NetBenefits_0.07_mid_SCGHG_5']=CY_7pct['$_Benefits_0.07_mid_SCGHG_5'] + CY_7pct_temp
CY_7pct['$_NetBenefits_0.07_mid_SCGHG_3']=CY_7pct['$_Benefits_0.07_mid_SCGHG_3'] + CY_7pct_temp
CY_7pct['$_NetBenefits_0.07_mid_SCGHG_25']=CY_7pct['$_Benefits_0.07_mid_SCGHG_25'] + CY_7pct_temp
CY_7pct['$_NetBenefits_0.07_mid_SCGHG_395']=CY_7pct['$_Benefits_0.07_mid_SCGHG_395'] + CY_7pct_temp

# third, discount payback values using payback as the source
# do not want to group these by MY/CY since we need output year-over-year
# payback needs a new discount function since it needs to be discounted to age 0
def discount_payback(col_source, rate, source, destination):
    for i in range(0,len(col_source)):
        col=col_source[i]
        destination[col]=source[col]/((1+rate)**(destination['CY']-destination['MY']+.5))
    return(destination)

payback_3pct_all=pd.DataFrame(payback, columns=['MY', 'age', 'CY'])
payback_7pct_all=pd.DataFrame(payback, columns=['MY', 'age', 'CY'])
discount_payback(payback_vars, 0.03, payback, payback_3pct_all)
discount_payback(payback_vars, 0.07, payback, payback_7pct_all)

# groupby MY; this is to do the cumulative sum of the discounted streams
grouped_3pct=payback_3pct_all.groupby(['MY'])
payback_3pct_all['$_CashCumOperCost_car_delta']=grouped_3pct['$_perVeh_Cash_car_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_3pct_all['$_CashCumOperCost_truck_delta']=grouped_3pct['$_perVeh_Cash_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_3pct_all['$_CashCumOperCost_fleet_delta']=grouped_3pct['$_perVeh_Cash_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_fleet_delta'].cumsum()
payback_3pct_all['$_Loan4yrCumOperCost_car_delta']=grouped_3pct['$_perVeh_Loan4yr_car_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_3pct_all['$_Loan4yrCumOperCost_truck_delta']=grouped_3pct['$_perVeh_Loan4yr_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_3pct_all['$_Loan4yrCumOperCost_fleet_delta']=grouped_3pct['$_perVeh_Loan4yr_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_fleet_delta'].cumsum()
payback_3pct_all['$_Loan5yrCumOperCost_car_delta']=grouped_3pct['$_perVeh_Loan5yr_car_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_3pct_all['$_Loan5yrCumOperCost_truck_delta']=grouped_3pct['$_perVeh_Loan5yr_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_3pct_all['$_Loan5yrCumOperCost_fleet_delta']=grouped_3pct['$_perVeh_Loan5yr_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_fleet_delta'].cumsum()
payback_3pct_all['$_Loan6yrCumOperCost_car_delta']=grouped_3pct['$_perVeh_Loan6yr_car_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_3pct_all['$_Loan6yrCumOperCost_truck_delta']=grouped_3pct['$_perVeh_Loan6yr_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_3pct_all['$_Loan6yrCumOperCost_fleet_delta']=grouped_3pct['$_perVeh_Loan6yr_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_3pct['$_perVeh_Fuel_fleet_delta'].cumsum()

grouped_7pct=payback_7pct_all.groupby(['MY'])
payback_7pct_all['$_CashCumOperCost_car_delta']=grouped_7pct['$_perVeh_Cash_car_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_7pct_all['$_CashCumOperCost_truck_delta']=grouped_7pct['$_perVeh_Cash_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_7pct_all['$_CashCumOperCost_fleet_delta']=grouped_7pct['$_perVeh_Cash_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_fleet_delta'].cumsum()
payback_7pct_all['$_Loan4yrCumOperCost_car_delta']=grouped_7pct['$_perVeh_Loan4yr_car_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_7pct_all['$_Loan4yrCumOperCost_truck_delta']=grouped_7pct['$_perVeh_Loan4yr_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_7pct_all['$_Loan4yrCumOperCost_fleet_delta']=grouped_7pct['$_perVeh_Loan4yr_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_fleet_delta'].cumsum()
payback_7pct_all['$_Loan5yrCumOperCost_car_delta']=grouped_7pct['$_perVeh_Loan5yr_car_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_7pct_all['$_Loan5yrCumOperCost_truck_delta']=grouped_7pct['$_perVeh_Loan5yr_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_7pct_all['$_Loan5yrCumOperCost_fleet_delta']=grouped_7pct['$_perVeh_Loan5yr_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_fleet_delta'].cumsum()
payback_7pct_all['$_Loan6yrCumOperCost_car_delta']=grouped_7pct['$_perVeh_Loan6yr_car_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_car_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_car_delta'].cumsum()
payback_7pct_all['$_Loan6yrCumOperCost_truck_delta']=grouped_7pct['$_perVeh_Loan6yr_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_truck_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_truck_delta'].cumsum()
payback_7pct_all['$_Loan6yrCumOperCost_fleet_delta']=grouped_7pct['$_perVeh_Loan6yr_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Maintenance_fleet_delta'].cumsum() + grouped_7pct['$_perVeh_Fuel_fleet_delta'].cumsum()

# Create the Output
writer=pd.ExcelWriter(outputFile)
calcs.to_excel(writer, 'Main Matrix')
MYSummary.to_excel(writer, 'MY_summary')

MY_3pct.to_excel(writer, 'MY_3pct')
MY_7pct.to_excel(writer, 'MY_7pct')
MY_SCC.to_excel(writer, 'MY_SCC')

CYSummary.to_excel(writer, 'CY_summary')
CY_3pct.to_excel(writer, 'CY_3pct')
CY_7pct.to_excel(writer, 'CY_7pct')
CY_SCC.to_excel(writer, 'CY_SCC')

maintenance.to_excel(writer, 'Maintenance', index=False)
payback.to_excel(writer, 'Payback', index=False)
payback_3pct_all.to_excel(writer, 'Payback_3pct', index=False)
payback_7pct_all.to_excel(writer, 'Payback_7pct', index=False)
#population.to_excel(writer, 'Populations')
Input_sum.to_excel(writer,'Inputs', index=False)
# py_version.to_excel(writer, 'script', index=False)
writer.save()
print('Complete')
