""" Python script to construct Carrying Capacity Diagrams for the Scaled-RRF

method for the 2015 Annual PM2.5 SIP update.

This works for old 2008 Annual Averge worksheets only.

Usage: python "this file".py

"""

import sys

import win32com.client

import os

import glob

import shutil

import csv

#

def scale(ws2, OC_scale, EC_scale, OTR_scale, NOx_scale, NH4_scale, SOx_scale):

index = [ 4, 5, 6, 7 ]

for k in index:

ws2.Range("I" + str(k)).Value = (1. - ws2.Range("B" + str(k)).Value) * NH4_scale

ws2.Range("J" + str(k)).Value = (1. - ws2.Range("C" + str(k)).Value) * NOx_scale

ws2.Range("K" + str(k)).Value = (1. - ws2.Range("D" + str(k)).Value) * SOx_scale

ws2.Range("L" + str(k)).Value = (1. - ws2.Range("E" + str(k)).Value) * EC_scale

ws2.Range("M" + str(k)).Value = (1. - ws2.Range("F" + str(k)).Value) * OC_scale

ws2.Range("N" + str(k)).Value = (1. - ws2.Range("G" + str(k)).Value) * OTR_scale

# Substact from unity and copy back.

ws2.Range("B" + str(k)).Value = 1. - ws2.Range("I" + str(k)).Value

ws2.Range("C" + str(k)).Value = 1. - ws2.Range("J" + str(k)).Value

ws2.Range("D" + str(k)).Value = 1. - ws2.Range("K" + str(k)).Value

ws2.Range("E" + str(k)).Value = 1. - ws2.Range("L" + str(k)).Value

ws2.Range("F" + str(k)).Value = 1. - ws2.Range("M" + str(k)).Value

ws2.Range("G" + str(k)).Value = 1. - ws2.Range("N" + str(k)).Value

#

def pbw(ws):

Row_List = [14, 15, 16, 17, 24, 25, 26, 27]

mystring_1 = "temp 1.0 1.0 1 0 rh 0. ammonium 0. sulphate nitrate 0. 0. 0. 0. 3 3 3 3 

3 \

30 1 3 2 3 3 3 4 3 5 3 6 3 7 3 8 3 9 3 10 3 \

11 3 12 3 13 3 14 3 15 3 16 3 17 3 18 3 19 3 20 3 \

21 3 22 3 23 3 24 3 25 3 26 3 27 3 28 3 29 3 30 3 "

# temp,rh are fixed.

temp,rh = 294.15, 0.35

mystring_2 = mystring_1.replace("temp", str(temp))

mystring_3 = mystring_2.replace("rh", str(rh))

#

Values = []

pbw = []

#

for i in Row_List:

cellnum = "%s" % (i)

nitrate = ws.Range("C" + cellnum).Value/62.0

sulphate = ws.Range("D" + cellnum).Value/96.0

ammonium = (nitrate + 2.0*sulphate)

mystring_4 = mystring_3.replace("nitrate", str(nitrate))

mystring_5 = mystring_4.replace("sulphate", str(sulphate))

mystring_6 = mystring_5.replace("ammonium", str(ammonium))

Values.append([mystring_6])

#

shutil.copy('eaim.dat.tmp', 'eaim.dat')

f = open("eaim.dat", 'a')

f.write("\n".join((item[0] for item in Values)))

f.write("\n")

f.close()

#

os.system("eaim.exe")

#
# Search the output for (PBW in grams) and save it to PBW.

fout = open("eaim.rs1", 'r')

for line in fout:

if line.find("H2O(aq)") > -1:

pbw.append(float(line.split()[2]))

fout.close()

#

k = 0

for j in Row_List:

cellnum = "%s" % (j)

ws.Range("H" + cellnum).Value = pbw[k]

k = k + 1

def main():

# Open each book, loop over the diagram(copy the appropriate sheet to a

# scratch sheet, scale RRFs, calculate PBW, copy the DV to another book,

# delete the sheet), and close the book. Make sure the book does not need any

# intervension when you open the book or copy sheets.

# Initialize Excel application.

xl = win32com.client.gencache.EnsureDispatch("Excel.Application")

# Visible make the workbook visible. This allows updates and saves later.

xl.Visible = 1

# win32.com looks for the file in "My Documents". The following few lines

# tells win32.com where the file is.

cwd = os.getcwd()

print "%s" % cwd

# 2012 & 2020 emissions - With external adjustments - Baseline.

# From Steve Zelinka on Feb 13, 2015 (V1.01) + rules/incentives from Patricia

# on January 15, 2015.

PM25_base = 65.956

PM25_future = 60.181

# PM is now split in to EC, OC, and OTR (=PM-EC-OC).

EC_base = 9.711

EC_future = 6.498

OC_base = 12.108

OC_future = 11.013

OTR_base = 44.137

OTR_future = 43.306

NOx_base = 332.191 # With Madera update

NOx_future = 206.432

SOx_base = 8.085

SOx_future = 7.849

# 2014-2005 Baseline from the old Plan.

PM25_2008_base = 86.0

PM25_2008 = 22.7

EC_2008_base = 11.4

EC_2008 = 3.9
OC_2008_base = 26.1

OC_2008 = 7.5

OTR_2008_base = 48.5

OTR_2008 = 11.2

NOx_2008_base = 575.4

NOx_2008 = 284.2

SOx_2008_base = 26.4

SOx_2008 = 1.8

# SOx is not changed in the Carrying Capacity Diagram. So, we can

# calculate that scaling factor here.

SOx_scale = ((SOx_base - SOx_future)/SOx_base)*(SOx_2008_base/SOx_2008)

print "SOx_scale %s" % SOx_scale

# Sheets in the Carrying Capacity Diagram.

Scale_Sheets = ['BAC-CTRL_2013', \

'BEP-CTRL_2013', \

'FSF-CTRL_2013', 'MAD-CTRL_2013', 'MAD-CTRL_2014', \

'FSH-CTRL_2013', 'TRA-CTRL_2013', \

'CLO-CTRL_2013', \

'M14-CTRL_2013', 'MAN-CTRL_2013', 'TUR-CTRL_2013', \

'MRM-CTRL_2013', \

'MRC-CTRL_2013', \

'SOH-CTRL_2013', \

'VCS-CTRL_2013', 'HAN-CTRL_2013' ]

# Design Values(DVs)are written into a csv file.

dv_file = "CC_Diagrams_percent_2020.csv"

dvWriter = csv.writer(open(dv_file, 'wb'), delimiter=',')

#

files = glob.glob('4K-SENS_2013\*.xlsm')

print "%s" % files

for f in files:

print "In the files loop"

inname = f

print "%s" % inname

# Open the workbooks and select the worksheets.

wb = xl.Workbooks.Open(cwd + "/" + inname)

number_sheets = wb.Sheets.Count

print "Number of sheets %s" % number_sheets

for ns in range(number_sheets):

print "Sheet number %s" % ns

wsname = wb.Sheets(ns+1).Name

if wsname in Scale_Sheets:

print "In the worksheet %s" % wsname

dvWriter.writerow(wsname)

ws = wb.Sheets(ns+1)

for pm in range(30, 101, 5):

# for pm in range(100, 101, 5):

# DV array is used to write DVs into the csv file.

DV = []

# We are scaling from the 2013 future value.

# OC_scale = EC_scale = OTR_scale = \

# ((PM25_base - PM25_future * 

pm/100.)/PM25_base)*(PM25_2008_base/PM25_2008)

OC_scale = \

((OC_base - OC_future * pm/100.)/OC_base)*(OC_2008_base/OC_2008)

EC_scale = \

((EC_base - EC_future * pm/100.)/EC_base)*(EC_2008_base/EC_2008)

OTR_scale = \

((OTR_base - OTR_future * 
pm/100.)/OTR_base)*(OTR_2008_base/OTR_2008)

print "OC_scale %s" % OC_scale

print "EC_scale %s" % EC_scale

print "OTR_scale %s" % OTR_scale

for nox in range(30, 101, 5):

# for nox in range(100, 101, 5):

NOx_scale = NH4_scale = \

((NOx_base - NOx_future * 

nox/100.)/NOx_base)*(NOx_2008_base/NOx_2008)

print "NOx_scale %s" % NOx_scale

wb.Sheets(ns+1).Copy(After=wb.Sheets(number_sheets))

ws2 = wb.Sheets(number_sheets+1)

ws2.Name = wsname + "_CC"

# Scale the RRFs. We could just pass two scaling factors, but we pass all

# just in case they were made different in the future.

scale(ws2, OC_scale, EC_scale, OTR_scale, NOx_scale, 

NH4_scale, \

SOx_scale)

# Calculate particle bound water (PBW).

pbw(ws2)

# Write the future DV onto a different sheet.

DV_value = ws2.Range("O28").Value

DV.append(DV_value)

xl.DisplayAlerts=False

ws2.Delete()

xl.DisplayAlerts=True

dvWriter.writerow(DV)

wb.Close(SaveChanges=1)

#

if __name__ == "__main__":

sys.exit(main())