#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon May 20 09:19:42 2024

@author: jdaniel
"""

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
plt.rcParams["font.family"] = "Arial"
#plt.rcParams["font.family"] = "HelveticaNeue"

def makefig1(kernel):
    #Figure 1
    fact = (1.-.1)/100 #account for fact that solvent use is removed from kernel
    fig, ax1 = plt.subplots(figsize=(5,2.5), dpi=300)
    plt.plot(kernel.cprodemiss[0,:]*100,color='red')
    plt.plot(kernel.cinstallemiss[0,:]/fact,color='orange')
    plt.plot(kernel.cactiveemiss[0,:]/fact,color='yellow')
    plt.plot(kernel.cinactiveemiss[0,:]/fact,color='grey')
    plt.plot(kernel.cdecomemiss[0,:]/fact,color='blue')
    plt.xlim(0,75); plt.ylim(0,20)
    plt.title('Cumulative Emission'); plt.ylabel('Gg'); plt.xlabel('Year After Production')
    ax1.tick_params(axis='both', direction='in', top=True, right=True)
    ax1.tick_params(axis='x', pad = 5)
    plt.savefig('figure1.png',bbox_inches='tight')
    plt.savefig('figure1.pdf',bbox_inches='tight')
    plt.show()
    
def getRange(ibot,itop,numyears,resultsmc,attrUse):
    tmparr = np.zeros([len(resultsmc),numyears])
    a, b = np.zeros(numyears), np.zeros(numyears)
    for i in range(0,len(resultsmc)):
        tmparr[i,:] = np.sum(getattr(resultsmc[i],attrUse),axis=0)
    for i in range(0,numyears):
        tmp = np.sort(tmparr[:,i])
        a[i] = tmp[ibot]
        b[i] = tmp[itop] 
    return(a,b)

def getOABank():
    file = '/Users/jdaniel/new/OA_banks_30Jun.xlsx'
    fileinfo = pd.read_excel(file,header=0)
    s = np.shape(fileinfo)
    bank = np.zeros(s)
    bank[:,:] = fileinfo.loc[:,:]
    return(bank)

def getLickley():
    file1 = '/Users/jdaniel/new/lickley_low.txt'
    file2 = '/Users/jdaniel/new/lickley_high.txt'
    fileinfo1 = pd.read_csv(file1, sep='\s+', header=None)
    fileinfo2 = pd.read_csv(file2, sep='\s+', header=None)
    s_low = np.shape(fileinfo1)
    s_high = np.shape(fileinfo2)
    low = np.zeros(s_low)
    high = np.zeros(s_high)
    low[:] = fileinfo1.iloc[:,:]
    high[:] = fileinfo2.iloc[:,:]
    highout = np.interp(low[:,0],high[:,0],high[:,1])
    return(low,highout)

def getWestern():
    file = 'HCFC-141b_Global_annual_emissions.xlsx'
    sheet1, sheet2 = 'AGAGE','NOAA'
    fileinfo1 = pd.read_excel(file,sheet_name=sheet1,header=12)
    fileinfo2 = pd.read_excel(file,sheet_name=sheet2,header=12)
    s_agage = np.shape(fileinfo1)
    s_noaa = np.shape(fileinfo2)
    print(s_agage)
    w_agage = np.zeros(s_agage)
    w_noaa = np.zeros(s_noaa)
    w_agage[:] = fileinfo1.iloc[:,:]
    w_noaa[:] = fileinfo2.iloc[:,:]
    print(np.shape(fileinfo1),np.shape(w_agage))
    return(w_agage,w_noaa)

def outputPresent(years,markets,regions,consump_bymarket,consump_bymarket_mc,results,resultsmc,weibull,weibullem):
    indexcheck = 20
    dim = results.solventemiss.shape
    emissions = results.solventemiss+results.installemiss+results.activeemiss+results.decomemiss+results.inactiveemiss
    cemissions = results.csolventemiss+results.cinstallemiss+results.cactiveemiss+results.cdecomemiss+results.cinactiveemiss
    emcheck = np.sum(emissions[:,0:indexcheck-1]) 
    cemcheck = np.sum(cemissions[:,indexcheck-1],axis=0)
    activebankcheck = np.sum(results.activebank[:,indexcheck],axis=0) 
    inactivebankcheck = np.sum(results.inactivebank[:,indexcheck],axis=0)
    
#Figure 2
    fig = plt.figure(figsize=(4,3), dpi=300)
    ax1 = fig.add_subplot()
    cregion = np.sum(consump_bymarket,axis=2)
    cregion_total = np.sum(cregion,axis=0)
    v = np.zeros(len(cregion[:,0]))
    vold = np.zeros(len(cregion[:,0]))
    csorttmp = np.argsort(cregion_total)
    regname = np.copy(regions)
    regname[5], regname[6], regname[7] = 'Middle East/\nNorth Africa','North America','Northeast Asia'
    y1 = np.zeros(5)
    x1 = [2008,1994,1996,2006,2008]
    y1 = np.array([20000,30000,82000,62000,90000])
    y1 = y1/1000.
    ctuple=[(1,0.,0.),(0,1,1),(.75,0,1.),(0,1,0),(1,.5,0),(1,1,0),(.5,.5,.5)]
    csort = csorttmp[::-1]
    for i in range(0,5):
        v = v + cregion[:,csort[i]]
        ax1.fill_between(years,vold/1000,v/1000,color=ctuple[i],alpha=.75)
        if (i != 4 and i !=2):
            plt.text(x1[i],y1[i],regname[csort[i]],fontsize=8)
        elif (i == 4):
            plt.annotate('Middle East/\nNorth Africa',fontsize=8,xy=(2011.2,95),xytext=(2014,110),\
                arrowprops=dict(color='black',headwidth=3.5,headlength=2,shrink=.01,width=.3))
        elif (i == 2):
            plt.text(x1[2],y1[2],'Europe\n+Japan',fontsize=8)
        vold = v
    ax1.fill_between(years,vold/1000,np.sum(cregion,axis=1)/1000,color=ctuple[5])
    plt.text(2008,93,'Others',fontsize=8,rotation=45)
    plt.title('Annual Regional Consumption'),plt.xlabel('Year'),plt.ylabel('Gg/yr')
    plt.xlim([1989,2027]),plt.ylim([0,160])
    plt.fill_between(years[34::],0,np.sum(np.sum(consump_bymarket[34::,:,:],axis=1),axis=1)/1000,color='white',hatch="||",edgecolor=(.1,.1,.1))
    plt.annotate('Projection\n(NE Asia)',fontsize=8,xy=(2025.,8),xytext=(2022,70),\
                 arrowprops=dict(color='black',headwidth=3.5,headlength=2,\
                                 width=.3),horizontalalignment='center')
    ax1.tick_params(axis='both', direction='in', top=True, right=True)
    ax1.tick_params(axis='x', pad = 5)
    plt.savefig('figure3.png',bbox_inches='tight')
    plt.show()
    
 #Figure 3
    plt.figure(figsize=(4,3), dpi=300)
    fig, ax1 = plt.subplots(figsize=(4,3), dpi=300)
    index=[5,0,10,3]
    index=[8,9,3]
    co = ['b','black','g','y']
    wyears = np.zeros(len(weibull[0,:]))
    for i in range(0,len(weibull[0,:])):
        wyears[i] = i
    for i in range(0,len(index)):
        plt.plot(wyears,weibull[index[i],:],color=co[i])
        plt.plot(wyears,np.cumsum(weibullem[index[i],:]),color=co[i],linestyle='dashed')
    plt.xlim([0,150])
    plt.ylim([-.0,1])
    plt.text(100,.49,'Block Pipe',fontsize=8)
    plt.text(100,.15,'Block Foam Slab',fontsize=8)
    plt.text(90,.27,'Discontinuous Panels',fontsize=8)
    plt.xlabel('Years After Installation'); plt.ylabel('Fraction of Installed Amount'); plt.title('Bank Remaining and Cumulative Emission')
    ax1.tick_params(axis='both', direction='in', top=True, right=True)
    ax1.tick_params(axis='x', pad = 5)
    plt.savefig('figure2.png',bbox_inches='tight')
    plt.show()

# # Figure 4
    fig = plt.figure(figsize=(4,3), dpi=300)
    fig, ax1 = plt.subplots(figsize=(4,3), dpi=300)
    c_market = np.sum(consump_bymarket,axis=1)
    v = np.zeros(len(consump_bymarket))
    vold = np.zeros(len(consump_bymarket))
    vtot = np.sum(c_market,axis=1)+.1
    index = [5,0,7,4,3,11]
    for i in range(0,len(index)):
        v = v + c_market[:,index[i]]
        ax1.fill_between(years,vold/vtot,v/vtot,color=ctuple[i],alpha=.75)
        vold = v
    ax1.fill_between(years,vold/vtot,1,color=ctuple[6])
    plt.title('HCFC-141b Markets')
    plt.ylabel('Fraction'); plt.xlabel('Year')
    plt.text(1994,.12,'Boardstock', fontsize=8)
    plt.text(1994,.31,'Domestic\nRefrigeration', fontsize=8)
    plt.text(2015,.13,'Pipe \nin Pipe', fontsize=8)
    plt.text(2014,.5,'Spray Foam', fontsize=8)
    plt.text(1995,.65,'Discontinuous Panels', fontsize=8)
    plt.text(1991,.78,'Continuous Panels', fontsize=8)
    plt.text(2000,.9,'Others', fontsize=8)
    plt.ylim([0,1.]); plt.xlim([1990,2021])
    ax1.tick_params(axis='both', direction='in', top=True, right=True)
    ax1.tick_params(axis='x', pad = 5)
    plt.savefig('figure4.png',bbox_inches='tight')
    plt.show()
        
#Figure 5
    percentile = 0.165
    #percentile = 0.025
    plt.figure(figsize=(4,3), dpi=300)
    fig, ax = plt.subplots(figsize=(4,3), dpi=300)
    w_agage,w_noaa = getWestern()
    itop = int((1-percentile) * len(resultsmc))
    ibot = int(percentile * len(resultsmc))
    index1 = [0]*dim[1]
    index2 = [0]*dim[1]
    yrout = [0.]*dim[1]
    for i in range(0,dim[1]-1):
        index1[i]=i
        index2[i]=i+1
        yrout[i] = years[0] + float(i) + .5
    yrout[-1] = yrout[-2] + 1
    attrUse = 'totalemiss'
    emisslow, emisshi = getRange(ibot,itop,dim[1],resultsmc,attrUse)
    plt.fill_between(yrout,emisslow/1000,emisshi/1000,color='black',alpha=.4)
    plt.fill_between(w_agage[:,2],w_agage[:,3]-w_agage[:,4],w_agage[:,3]+w_agage[:,4],color=(.7,0,.1),alpha=.5)
    plt.fill_between(w_noaa[:,2],w_noaa[:,3]-w_noaa[:,4],w_noaa[:,3]+w_noaa[:,4],color=(.1,.1,.7),alpha=.5)
    plt.plot(yrout,np.sum(resultsmc[0].totalemiss/1000,axis=0),color='black',linewidth = 2)
    plt.xlim([1990,2022])
    plt.ylim([0,80])
    plt.title('Annual Emission')
    plt.xlabel('Year')
    plt.ylabel('Emission (Gg/yr)')
    plt.fill_between([2000,2003],[20,20],[24,24],color='black',alpha=.4)
    plt.plot([2000,2003],[22,22],color='black')
    plt.text(2003.5,21,'This Work',fontsize=7)
    plt.fill_between([2000,2003],[14,14],[18,18],color=(.1,.1,.7),alpha=.5)
    plt.text(2003.5,15,'From NOAA, Western (2022)',fontsize=7)
    plt.fill_between([2000,2003],[8,8],[12,12],color=(.7,.0,.1),alpha=.5)
    plt.text(2003.5,9,'From AGAGE, Western (2022)',fontsize=7)
    ax.tick_params(axis='both', direction='in', top=True, right=True)
    ax.tick_params(axis='x', pad = 5)
    plt.savefig('figure5.png',bbox_inches='tight')
    plt.show()

    vals = np.zeros((len(resultsmc),200))
    low1, hi1 = 1998, 2002
    low2, hi2 = 2011, 2013
    low3, hi3 = 2004, 2007
    ir1 = np.where((w_agage[:,2] >= low1) & (w_agage[:,2] <= hi1))
    ir2 = np.where((w_agage[:,2] >= low2) & (w_agage[:,2] <= hi2))
    ir3 = np.where((w_agage[:,2] >= low3) & (w_agage[:,2] <= hi3))
    sumlow1 = np.mean(w_agage[ir1,3]-w_agage[ir1,4])
    sumlow2 = np.mean(w_agage[ir2,3]-w_agage[ir2,4])
    sumlow3 = np.mean(w_agage[ir3,3]-w_agage[ir3,4])
    sumhi1 = np.mean(w_agage[ir1,3]+w_agage[ir1,4])
    sumhi2 = np.mean(w_agage[ir2,3]+w_agage[ir2,4])
    sumhi3 = np.mean(w_agage[ir3,3]+w_agage[ir3,4])
    is1 = np.where((years >= low1) & (years <= hi1))
    is2 = np.where((years >= low2) & (years <= hi2))
    is3 = np.where((years >= low3) & (years <= hi3))
    iuse = []#np.array([])
    for i in range(0,len(resultsmc)):
        vals[i,:] = np.sum(getattr(resultsmc[i],'totalemiss'),axis=0)/1000.
        vsum1, vsum2, vsum3 = np.mean(vals[i,is1]),np.mean(vals[i,is2]),\
                    np.mean(vals[i,is3])
        if ((vsum1 >= sumlow1) & (vsum1 <= sumhi1) & (vsum2 >= sumlow2) & (vsum2 <= sumhi2) &\
            (vsum3 >= sumlow3) & (vsum3 <= sumhi3)):
            iuse = np.append(iuse,[int(i)])
    
    s = len(iuse)
    for i in range(0,len(iuse)):
        if (i == 0):
            iuseint = [int(iuse[i])]
        else:
            iuseint = np.append(iuseint,int(iuse[i]))
    

 
#Figure 6
    plt.figure(figsize=(4,3), dpi=300)
    fig, ax = plt.subplots(figsize=(4,3), dpi=300)
    fact=1000
    a00 = (resultsmc[0].activebank[0,:])/fact
    a01 = a00 + (resultsmc[0].activebank[3,:])/fact
    a02 = a01 + (resultsmc[0].activebank[4,:])/fact
    a03 = a02 + (resultsmc[0].activebank[5,:])/fact
    a04 = a03 + (resultsmc[0].activebank[7,:])/fact
    plt.fill_between(yrout,0,a00,color=(.85,.85,.85))
    plt.fill_between(yrout,a00,a01,color=(.75,.75,.75))
    plt.fill_between(yrout,a01,a02,color=(.65,.65,.65))
    plt.fill_between(yrout,a02,a03,color=(.55,.55,.55))
    plt.fill_between(yrout,a03,a04,color=(.45,.45,.45))
    a1 = np.sum(resultsmc[0].activebank,axis=0)/fact
    plt.fill_between(yrout,a04,a1,color=(.3,.3,.3))
    a2 = np.sum(resultsmc[0].inactivebank,axis=0)/fact
    a3 = np.sum(resultsmc[0].cactiveemiss,axis=0)/fact
    a4 = np.sum(resultsmc[0].cinactiveemiss,axis=0)/fact+np.sum(resultsmc[0].cdecomemiss,axis=0)/fact
    a5 = np.sum(resultsmc[0].cinstallemiss,axis=0)/fact
    a6 = np.sum(resultsmc[0].csolventemiss,axis=0)/fact
    a7 = np.sum(resultsmc[0].cprodemiss,axis=0)/fact
    plt.fill_between(yrout,a1,a1+a2,color=(0,.7,.2),alpha=.8)
    plt.fill_between(yrout,a1+a2,a1+a2+a3,color=(.8,.8,0),alpha=.8)
    plt.fill_between(yrout,a1+a2+a3,a1+a2+a3+a4,color=(.8,0,0),alpha=.8)
    plt.fill_between(yrout,a1+a2+a3+a4,a1+a2+a3+a4+a5,color=(0.75,0,1.),alpha=.8)
    plt.fill_between(yrout,a1+a2+a3+a4+a5,a1+a2+a3+a4+a5+a6,color=(.9,.45,0),alpha=.7)
    plt.fill_between(yrout,a1+a2+a3+a4+a5+a6,a1+a2+a3+a4+a5+a6+a7,color=(0,.45,.9),alpha=.8)
    plt.xlim([1990,2050]), plt.ylim([0,3.1E+3])
    plt.title('Life Cycle Evolution')
    plt.xlabel('Year'), plt.ylabel('Cumulative Emission, Bank (Gg)')
    plt.text(2005,50,'AB Domestic Refrig.',fontsize=7)
    plt.text(2020,375,'AB Spray Foam',fontsize=7,rotation=-12)
    plt.text(2025,30,'AB Discontinuous Panels',fontsize=7,rotation=-5)
    plt.text(2002,620,'Active Bank (AB) Other',fontsize=7,rotation=30,color='white')
    plt.text(2020,680,'AB Pipe-in-Pipe', fontsize=7,rotation=-20)
    plt.text(2003,460,'AB Boardstock', fontsize=7,rotation=25)
    plt.text(2030,1100,'Inactive Bank',fontsize=7)
    plt.text(2030,1660,'Active Bank Emission',fontsize=7)
    plt.text(2049.5,1900,'Decommissioning and Inactive\nBank Emissions',fontsize=7,horizontalalignment='right')
    plt.text(2030,2300,'Installation Emission',fontsize=7)
    plt.text(2030,2680,'Solvent Emission',fontsize=7)
    plt.text(2030,2900,'Production Emission',fontsize=7)
    ax.tick_params(axis='y', direction='in', right=True)
    plt.savefig('figure6.png',bbox_inches='tight')
    plt.show()
    
    for j in range(0,len(markets)):
        akeep = np.zeros((len(iuse),200))
        for i in range(0,len(iuse)):
            ii = int(iuse[i])
            a1 = resultsmc[ii].activebank[j,:]/fact
            akeep[i,:] = a1
        activekeep = np.zeros((len(resultsmc),200))
        inactivekeep = np.copy(activekeep)
    for i in range(0,len(resultsmc)):
        activekeep[i,:] = np.sum(resultsmc[i].activebank,axis=0)/fact
        inactivekeep[i,:] = np.sum(resultsmc[i].inactivebank,axis=0)/fact
    
    ii = np.where(np.min(activekeep,axis=1) < 0)
    jj = np.where(np.min(inactivekeep,axis=1) < 0)
    print('Indices when active bank is negative')
    print(ii,np.shape(ii))
    print('Indices when inactive bank is negative')
    print(jj,np.shape(jj))
    print(inactivekeep[jj])
    kk = np.where(activekeep < 0)
    print(kk,np.shape(kk))
    
#Figure 7    
    percentile = 0.05/2
#    percentile = .165
    itop = int((1-percentile) * len(resultsmc))
    ibot = int(percentile * len(resultsmc))
    fig, (ax,ax2) = plt.subplots(2,1,figsize=(4,6), dpi=300, constrained_layout=True)
    ax = plt.subplot2grid((5,1),(0,0),rowspan=3)
    ax2 = plt.subplot2grid((5,1),(3,0),rowspan=2)
    lickley_low, lickley_high = getLickley()
    OA_bank = getOABank()
    attrUse = 'activebank'
    activearraybot, activearraytop = getRange(ibot,itop,dim[1],resultsmc,attrUse)
    attrUse = 'totalbank'
    inactivearraybot, inactivearraytop = getRange(ibot,itop,dim[1],resultsmc,attrUse)
    ax.plot(yrout,np.sum(resultsmc[0].activebank/1000,axis=0),color='orange')
    tmpbank = np.sum(resultsmc[0].activebank/1000,axis=0)+np.sum(resultsmc[0].inactivebank/1000,axis=0)
    ax.plot(yrout,tmpbank,color='black')
    ax.fill_between(yrout,inactivearraybot/1000,inactivearraytop/1000,alpha=.5,color='black')
    ax.fill_between(yrout,activearraybot/1000,activearraytop/1000,alpha=.5,color='orange')
    ax.fill_between(lickley_low[:,0],lickley_low[:,1],lickley_high,alpha=.3)
    ax.plot(OA_bank[70::,0],OA_bank[70::,7]/1000,color='blue')
    ax.set_xlim([1990,2050]); ax.set_title('Bank Sizes')
    ax.text(1991,1870,"(a)", weight = 'bold')
    ax2.text(1991,1110,"(b)", weight = 'bold')
    ax.set_ylim([0,2000])
    ax.set_ylabel('Bank Size (Gg)')
    ax.annotate('Active Bank\nThis Work',fontsize=8,xy=(2025,1170),xytext=(2022,750),horizontalalignment='right',\
        arrowprops=dict(color='black',headwidth=3.5,headlength=2,shrink=.01,width=.3))
    ax.annotate('Total Bank\nThis Work',fontsize=8,xy=(2045,1450),xytext=(2040,1800),horizontalalignment='center',\
        arrowprops=dict(color='black',headwidth=3.5,headlength=2,shrink=.01,width=.3))
    ax.annotate('Total Bank,\nLickley et al. (2022)',fontsize=8,xy=(2015,1800),xytext=(2000,1600),horizontalalignment='center',\
        arrowprops=dict(color='black',headwidth=3.5,headlength=2,shrink=.01,width=.3))
    ax.annotate('Total Bank\nWMO (2022)',fontsize=8,xy=(2045,825),xytext=(2040.,300),horizontalalignment='right',\
        arrowprops=dict(color='black',headwidth=3.5,headlength=2,shrink=.01,width=.3))
    ax.tick_params(axis='both', direction='in', top=True, right=True)
    ax.tick_params(axis='x', pad = 5)
    print(iuse,np.shape(iuse))
    tmp = np.zeros([len(iuse),len(resultsmc[0].activebank[0,:])])
    tmp1 = np.zeros([len(iuse),len(resultsmc[0].activebank[0,:])])
    tmp2 = np.zeros([len(iuse),len(resultsmc[0].inactivebank[0,:])])
    for i in range(0,len(iuse)):
        tmp[i,:] = np.sum(resultsmc[iuseint[i]].activebank,axis=0)/1000+\
            np.sum(resultsmc[iuseint[i]].inactivebank,axis=0)/1000
        tmp1[i,:] = np.sum(resultsmc[iuseint[i]].activebank,axis=0)/1000
        tmp2[i,:] = np.sum(resultsmc[iuseint[i]].inactivebank,axis=0)/1000
    print('Total Bank numbers in 2020: %12.2f %12.2f%12.2f' % (np.mean(tmp[:,31]),tmpbank[31],OA_bank[70,7]/1000))
    ax.plot(yrout,np.mean(tmp,axis=0),color='lightgrey')
    ax.plot(yrout,np.mean(tmp1,axis=0),color='lightgrey')
    ax2.plot(yrout,np.mean(tmp2,axis=0),color='lightgrey')
    attrUse = 'inactivebank'
    inactivearraybot, inactivearraytop = getRange(ibot,itop,dim[1],resultsmc,attrUse)
    ax2.plot(yrout,np.sum(resultsmc[0].inactivebank/1000,axis=0),color='black')
    ax2.fill_between(yrout,inactivearraybot/1000,inactivearraytop/1000,alpha=.5,color='black')
    ax2.set_xlim(1990,2050);ax2.set_ylim(0,1200)
    ax2.tick_params(axis='both', direction='in', top=True, right=True)
    ax2.tick_params(axis='x', pad = 5)
    ax2.set_ylabel('Bank Size (Gg)')
    ax2.set_xlabel('Year')
    ax2.annotate('Only Incl. "accurate"\nEmission Pathways',fontsize=8,xy=(2045,750),xytext=(2043.,200),horizontalalignment='right',\
        arrowprops=dict(color='black',headwidth=3.5,headlength=2,shrink=.01,width=.3))
    ax2.annotate('Inactive Bank\nThis Work',fontsize=8,xy=(2032,800),xytext=(2017.,1000),horizontalalignment='right',\
        arrowprops=dict(color='black',headwidth=3.5,headlength=2,shrink=.01,width=.3))
    plt.savefig('figure7.png',bbox_inches='tight')
    plt.show()    

    return(iuseint)