#!/usr/bin/env python

###############################################################################
## This script reads in 250-hPa zonal wind EOF patterns and projects the     ##
## real-time data/forecasts from a particular case onto the EOF phase space  ##
###############################################################################

import matplotlib as mpl
mpl.use('pdf')
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import time
import datetime
import netCDF4
from numpy import *
from pylab import *

day = ['01','02','03','04','05','06','07','08','09','10','11','12','13','14','15','16','17','18','19','20','21','22','23','24','25','26','27','28','29','30','31'] 
hour = ['00','06','12','18'] 
mon = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'] 
monnum = ['1','2','3','4','5','6','7','8','9','10','11','12'] 
mtimes = [124,112,124,120,124,120,124,124,120,124,120,124] 
tstart = [0,124,236,360,480,604,724,848,972,1092,1216,1336] 
tstartleap = [0,124,240,364,484,608,728,852,976,1096,1220,1340] 
mondays = [31,28,31,30,31,30,31,31,30,31,30,31] 
mondaysleap = [31,29,31,30,31,30,31,31,30,31,30,31] 
lonpts = 720 
latpts = 361 
lonpts1deg = 360
latpts1deg = 181

####################### Specify the start times #########################

hr = 0 
#year = 2017 
#monindex = 7
#dy = 8

year = datetime.date.today().strftime("%Y") 
monindex = datetime.date.today().strftime("%m") 
dy = datetime.date.today().strftime("%d") 

year = int(year)
monindex = int(monindex)
dy = int(dy)

monstr = monnum[monindex-1] 
daystr = day[dy-1]  
dy = dy - 11      ### to go back 11 days ###

labelanl = linspace(dy,dy+11,12)   ### create labels for plots ###

if (dy <= 0) :     ### Adjust the dates if dy <= 0 ### 
	monindex = monindex - 1 
	if (monindex == 0) :
		monindex = 12 
		year = year - 1 
	if ((year%4) == 0) :
		if (monindex == 2) :
			dy = dy + mondaysleap[monindex-1] 
		else :
			dy = dy + mondaysleap[monindex-1]
	else :
		dy = dy + mondays[monindex-1] 

for i in range(len(labelanl)) :   ### Adjust labels for plots if straddling two months ###
	if (labelanl[i]) <= 0 :
		if ((year%4) == 0) :
			if (monindex == 2) :
				labelanl[i] = labelanl[i] + mondaysleap[monindex-1] 
			else :
				labelanl[i] = labelanl[i] + mondaysleap[monindex-1] 
		else :
			labelanl[i] = labelanl[i] + mondays[monindex-1] 

print ('Year: %d' % year)
print ('Month: %d' % monindex)
print ('Analysis start day: %d' % dy)
print labelanl

if ((year%4) == 0) :
	tcount = tstartleap[monindex-1] + (4*(dy-1) + hr + 1) 
else :
	tcount = tstart[monindex-1] + (4*(dy-1) + hr + 1) 

mcount = (4*(dy-1) + hr + 1) 
month = mon[monindex-1] 
montimes = mtimes[monindex-1] 

### Domain bounds over the N. Pacific Basin
latS = 10    
latN = 80  
lonW = 100 
lonE = 240   

####################### Initialize some indices #########################

dlat = abs(latN - latS) 
dlon = abs(lonE - lonW) 
resolution = 1.0 
steplat = int((dlat/resolution)+1.0) 
steplon = int((dlon/resolution)+1.0)  
totpts = steplat*steplon 
Lat = linspace(latS,latN,steplat)  
Lon = linspace(lonW,lonE,steplon) 

latSind = ((latS-(-90)))+1 
latNind = ((latN-(-90)))+1 
lonWind = ((lonW-(0)))+1 
lonEind = ((lonE-(0)))+1 
numpts = (((latN-latS))+1)*(((lonE-lonW))+1) 

####################### Read in EOF patterns ###############################

eof1 = []
eof2 = []
eof3 = []
for j in range(steplon):
    eof1.append([])
    eof2.append([])
    eof3.append([])
    for i in range(steplat):
        eof1[j].append(0.)
        eof2[j].append(0.)
        eof3[j].append(0.)

file = open("250uwnd_eof_9mo1deg.out", "r")
data = file.readlines()

count = 0 
for i in range(steplat) :
	for j in range(steplon) :
		eof1[j][(steplat-1)-i] = float(data[count]) 
		count = count + 1 
for i in range(steplat) :
	for j in range(steplon) :
		eof2[j][(steplat-1)-i] = -1.0*float(data[count]) 
		count = count + 1 
for i in range(steplat) :
	for j in range(steplon) :
		eof3[j][(steplat-1)-i] = float(data[count])
		count = count + 1 

###################### Read in mean U-wind values  ###############################

filinstr = 'uwnd'+month+'.out' 
file = open(filinstr,"r") 
data = file.readlines()

meanpre = []
meanwind = []
meanwind1deg = []
for j in range(lonpts) :
	meanpre.append([])
	meanwind.append([])
	for i in range(latpts) :
		meanpre[j].append([]) 
		meanwind[j].append([]) 
		for t in range(montimes) :
			meanpre[j][i].append(0.)
			meanwind[j][i].append(0.)
			
for j in range(lonpts1deg) :
	meanwind1deg.append([])
	for i in range(latpts1deg) :
		meanwind1deg[j].append([])
		for t in range(montimes) :
			meanwind1deg[j][i].append(0.) 
			
count = 0 ;
for t in range(montimes) :
	for i in range(latpts) :
		for j in range(lonpts) :
			meanpre[j][(latpts-1)-i][t] = float(data[count]) 
			count = count + 1 

if ((year%4) == 0) :    ### If during February of a leap year, grab leap year data ###
	if (monindex == 2) :
		montimes = 116 
		filinstr = 'uwnd0229.out' 
		file = open(filinstr,"r")
		data = file.readlines() 
		
		for j in range(lonpts) :
			for i in range(latpts) :
				for l in range(4) :
					meanpre[j][i].append(0.) 
					meanwind[j][i].append(0.)
					
		for j in range(lonpts1deg) :
			for i in range(latpts1deg) :
				for l in range(4) :
					meanwind1deg[j][i].append(0.) 
					
		count = 0 
		for l in range(4) :
			for i in range(latpts) :
				for j in range(lonpts) :
					meanpre[j][(latpts-1)-i][t+(l+1)] = float(data[count]) 
					count = count + 1 

for t in range(montimes) :    ### Change so that 0 lon index corresponds to OE longitude ###
	for i in range(latpts) :
		for j in range(360) :
			meanwind[j+360][i][t] = meanpre[j][i][t] 
		for j in range(360) :
			meanwind[j][i][t] = meanpre[j+360][i][t] 

for t in range(montimes) :    ### Coarsen to 1 degree resolution ###
	for i in range(latpts1deg) :
		for j in range(lonpts1deg) :
			meanwind1deg[j][i][t] = meanwind[(j*2)][(i*2)][t] 

################# Read in GFS analysis data from today's forecast ###########################	

upert = []
for j in range(steplon) :
	upert.append([])
	for i in range(steplat) :
		upert[j].append([])
		for t in range(85) :
			upert[j][i].append(0.)

for t in range(85) :  ## 11-days of GFS analysis + 10 days of GFS forecast in this file ##

	filestr = '/nfs/winterslab_rit/realtime/rawfiles/gfsanl_'+monstr+daystr+'.nc' 
	ufile = netCDF4.Dataset(filestr,"r") 
	uwnd = ufile.variables['UGRD_P0_L100_GLL0'][t,:,:] 
	lat = ufile.variables['lat_0'][:]
	lon = ufile.variables['lon_0'][:]
	ufile.close()
	
	uwndnew = [] ;    ## because GFS data runs from 90N lat to 90S lat
	for j in range(lonpts1deg) :
		uwndnew.append([])
		for i in range(latpts1deg) :
			uwndnew[j].append(0.) 
			
	for i in range(latpts1deg) :
		for j in range(lonpts1deg) : 
			uwndnew[j][(latpts1deg-1)-i] = uwnd[i][j] 

	ilat = 0     ### Grab data only over the N. Pacific Basin
	for i in range((latSind-1),latNind) :
		ilon = 0 
		for j in range((lonWind-1),lonEind) :
			upert[ilon][ilat][t] = uwndnew[j][i] - meanwind1deg[j][i][mcount-1] ;
			ilon = ilon + 1 
		ilat = ilat + 1 
	tcount = tcount + 1 
	mcount = mcount + 1 	
	
	if (mcount > montimes) :    ##### If advancing to next month update climo file ######
		newmonindex = monindex + 1 
		if (newmonindex > 12) :
			newmonindex = 1 
			year = year + 1 
		newmontimes = mtimes[newmonindex-1] 
		newmonth = mon[newmonindex-1] 
	
		filinstr = 'uwnd'+newmonth+'.out' 
		file = open(filinstr,"r") 
		data = file.readlines()
		
		meanpre = []
		meanwind = []
		meanwind1deg = []
		for j in range(lonpts) :
			meanpre.append([])
			meanwind.append([])
			for i in range(latpts) :
				meanpre[j].append([]) 
				meanwind[j].append([]) 
				for k in range(newmontimes) :
					meanpre[j][i].append(0.)
					meanwind[j][i].append(0.)
			
		for j in range(lonpts1deg) :
			meanwind1deg.append([])
			for i in range(latpts1deg) :
				meanwind1deg[j].append([])
				for k in range(newmontimes) :
					meanwind1deg[j][i].append(0.) 
					
		count = 0 ;
		for k in range(newmontimes) :
			for i in range(latpts) :
				for j in range(lonpts) :
					meanpre[j][(latpts-1)-i][k] = float(data[count]) 
					count = count + 1 
					
		if ((year%4) == 0) :
			if (newmonindex == 2) :
				newmontimes = 116 
				filinstr = 'uwnd0229.out' 
				file = open(filinstr,"r")
				data = file.readlines() 
		
				for j in range(lonpts) :
					for i in range(latpts) :
						for l in range(4) :
							meanpre[j][i].append(0.) 
							meanwind[j][i].append(0.)
					
				for j in range(lonpts1deg) :
					for i in range(latpts1deg) :
						for l in range(4) :
							meanwind1deg[j][i].append(0.) 
					
				count = 0 
				for l in range(4) :
					for i in range(latpts) :
						for j in range(lonpts) :
							meanpre[j][(latpts-1)-i][k+(l+1)] = float(data[count]) 
							count = count + 1 
		
		for k in range(newmontimes) :
			for i in range(latpts) :
				for j in range(360) :
					meanwind[j+360][i][k] = meanpre[j][i][k] 
				for j in range(360) :
					meanwind[j][i][k] = meanpre[j+360][i][k] 

		for k in range(newmontimes) :
			for i in range(latpts1deg) :
				for j in range(lonpts1deg) :
					meanwind1deg[j][i][k] = meanwind[(j*2)][(i*2)][k] 
		mcount = 1

######################### Weight data by cos(lat)  ################################

numdates = len(upert[0][0]) 

for k in range(numdates) :
	for i in range(steplat) :
		for j in range(steplon) :
			# weight by sqrt(cos(lat)) #
			upert[j][i][k] = upert[j][i][k]*sqrt(cos(Lat[i]*(pi/180.))) 

########################## Two-dimensionalize data ##################################

upert2D = []
eof12D = []
eof22D = []
eof32D = []
count = 0 
for i in range(steplat) :
	for j in range(steplon) :
		upert2D.append([]) 
		eof12D.append(0.)
		eof22D.append(0.)
		eof32D.append(0.)
		for k in range(numdates) :
			upert2D[count].append(0.) 
		count = count + 1		

for k in range(numdates) :
	count = 0 
	for i in range(steplat) :
		for j in range(steplon) :
			upert2D[count][k] = upert[j][i][k] 
			eof12D[count] = eof1[j][i] 
			eof22D[count] = eof2[j][i] 
			eof32D[count] = eof3[j][i] 
			count = count + 1 

##################### Project Uwnd Perts onto EOF patterns #############################

#### need to divide by sqrt(eigenvalue) (i.e., lamda)) ########

for i in range(len(eof12D)) :     
	eof12D[i] = eof12D[i]/(sqrt(1.3253e5)) 
	eof22D[i] = eof22D[i]/(sqrt(1.0022e5)) 
	eof32D[i] = eof32D[i]/(sqrt(0.8083e5)) 
	
#### transpose upert2D ####
upert2D_trans = []
pc1 = []
pc2 = []
pc3 = []
for k in range(numdates) :
	upert2D_trans.append([]) 
	for i in range(len(upert2D)) :
		upert2D_trans[k].append(upert2D[i][k]) 

for k in range(numdates) :
	sumpc1 = 0 
	sumpc2 = 0
	sumpc3 = 0
	for i in range(len(eof12D)) :    ### upert2D'*eof12D - matrix multiplication ###
		sumpc1 = sumpc1 + (upert2D_trans[k][i]*eof12D[i])    ### Creates an 85x1 array of PCs ###
		sumpc2 = sumpc2 + (upert2D_trans[k][i]*eof22D[i])
		sumpc3 = sumpc3 + (upert2D_trans[k][i]*eof32D[i])
	pc1.append(sumpc1)    
	pc2.append(sumpc2)
	pc3.append(sumpc3)

for k in range(numdates) :
	pc1[k] = pc1[k]/(sqrt(1.3253e5)) 
	pc2[k] = pc2[k]/(sqrt(1.0022e5)) 
	pc3[k] = pc3[k]/(sqrt(0.8083e5)) 

#### smooth u-pert over the previous/next 1 day ####

upertsmoo = []
for j in range(steplon) :
	upertsmoo.append([]) 
	for i in range(steplat) :
		upertsmoo[j].append([])
		for k in range(numdates-8) :
			upertsmoo[j][i].append(0.)		

for i in range(steplat) :
	for j in range(steplon) :
		count = 0 
		for k in range(4,(numdates-4)) :
			weight = 5.0     
			upertsum = upert[j][i][k]*weight 
			for l in range(1,5) :   ## do a weighted average that decays in a pyramidal fashion over +/- 24h ##
				weight = weight - 1.0 
				upertsum = upertsum + upert[j][i][k-l]*weight + upert[j][i][k+l]*weight 
			upertsmoo[j][i][count] = upertsum/25.0 
			count = count + 1 

###### smooth PCs over the previous/next 1 day ####

pc1smoo = []
pc2smoo = []
pcs = []
count = 0 
for k in range(4,(numdates-4)) :     ## 5:(numdates-4) for +/- 1 day ##
	weight = 5.0                     ## max weight = 5 for +/- 1 day ##
	pc1sum = pc1[k]*weight 
	pc2sum = pc2[k]*weight 
	for l in range(1,5) :
		weight = weight - 1.0 
		pc1sum = pc1sum + pc1[k-l]*weight + pc1[k+l]*weight 
		pc2sum = pc2sum + pc2[k-l]*weight + pc2[k+l]*weight 
	pc1smoo.append(pc1sum/25.0)      ##  25.0 for sum of weights +/- 1 day ##
	pc2smoo.append(pc2sum/25.0) 
	count = count + 1 

for i in range(len(pc1smoo)) :
	pcs.append([])
	pcs[i].append(pc1smoo[i]) 
	pcs[i].append(pc2smoo[i])

########################################################################
######### Now calculate PCs for all the GEFS ensemble members ##########
########################################################################

ensnum = ['0','1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20'] 

################### Read in mean values ################################

#hr = 0 
#year = 2017 
#monindex = 7
#dy = 8
year = datetime.date.today().strftime("%Y") 
monindex = datetime.date.today().strftime("%m") 
dy = datetime.date.today().strftime("%d") 

year = int(year)
monindex = int(monindex)
dy = int(dy)

labelfcast = linspace(dy,dy+9,10) 
for i in range(len(labelfcast)) :
	if ((year%4) == 0) :
		if (labelfcast[i] > mondaysleap[monindex-1]) :
			labelfcast[i] = labelfcast[i] - mondaysleap[monindex-1] 
	else :
		if (labelfcast[i] > mondays[monindex-1]) :
			labelfcast[i] = labelfcast[i] - mondays[monindex-1] 

print labelfcast

pcsens = []
for e in range(21) :   ## change to 21 after completing script ##
	print ('Ensemble Number: %d' % e)
	#year = 2017 
	#monindex = 7
	#dy = 8
	
	year = datetime.date.today().strftime("%Y") 
	monindex = datetime.date.today().strftime("%m") 
	dy = datetime.date.today().strftime("%d") 

	year = int(year)
	monindex = int(monindex)
	dy = int(dy)
	hr = 1 ;   ## to start at F06 ##
	
	print ('Year: %d' % year)
	print ('Month: %d' % monindex)
	print ('Day: %d' % dy) 
	
	####################### Read in Mean U-wind data #######################
	
	mcount = (4*(dy-1) + hr + 1) 
	month = mon[monindex-1] ;
	montimes = mtimes[monindex-1] 
	
	filinstr = 'uwnd'+month+'.out' 
	file = open(filinstr,"r") 
	data = file.readlines()

	meanpre = []
	meanwind = []
	meanwind1deg = []
	for j in range(lonpts) :
		meanpre.append([])
		meanwind.append([])
		for i in range(latpts) :
			meanpre[j].append([]) 
			meanwind[j].append([]) 
			for t in range(montimes) :
				meanpre[j][i].append(0.)
				meanwind[j][i].append(0.)
				
	for j in range(lonpts1deg) :
		meanwind1deg.append([])
		for i in range(latpts1deg) :
			meanwind1deg[j].append([])
			for t in range(montimes) :
				meanwind1deg[j][i].append(0.) 
			
	count = 0 ;
	for t in range(montimes) :
		for i in range(latpts) :
			for j in range(lonpts) :
				meanpre[j][(latpts-1)-i][t] = float(data[count]) 
				count = count + 1 
	
	if ((year%4) == 0) :    ### If during February of a leap year, grab leap year data ###
		if (monindex == 2) :
			montimes = 116 
			filinstr = 'uwnd0229.out' 
			file = open(filinstr,"r")
			data = file.readlines() 
			
			for j in range(lonpts) :
				for i in range(latpts) :
					for l in range(4) :
						meanpre[j][i].append(0.) 
						meanwind[j][i].append(0.)
						
			for j in range(lonpts1deg) :
				for i in range(latpts1deg) :
					for l in range(4) :
						meanwind1deg[j][i].append(0.) 
						
			count = 0 
			for l in range(4) :
				for i in range(latpts) :
					for j in range(lonpts) :
						meanpre[j][(latpts-1)-i][t+(l+1)] = float(data[count]) 
						count = count + 1 
	
	for t in range(montimes) :    ### Change so that 0 lon index corresponds to OE longitude ###
		for i in range(latpts) :
			for j in range(360) :
				meanwind[j+360][i][t] = meanpre[j][i][t] 
			for j in range(360) :
				meanwind[j][i][t] = meanpre[j+360][i][t] 
	
	for t in range(montimes) :    ### Coarsen to 1 degree resolution ###
		for i in range(latpts1deg) :
			for j in range(lonpts1deg) :
				meanwind1deg[j][i][t] = meanwind[(j*2)][(i*2)][t] 
				
	######################## Read in GEFS ensemble data for today's forecast ##############################
	
	upertens = []
	for j in range(steplon) :
		upertens.append([])
		for i in range(steplat) :
			upertens[j].append([])
			for t in range(40) :
				upertens[j][i].append(0.)
	
	for t in range(40) : ## one-day buffer on either side of nine-day forecast ##
		
		filestr = '/nfs/winterslab_rit/realtime/rawfiles/gefsens'+ensnum[e]+'_'+monstr+daystr+'.nc' 
		ufile = netCDF4.Dataset(filestr,"r") 
		uwnd = ufile.variables['UGRD_P1_L100_GLL0'][t,:,:] 
		lat = ufile.variables['lat_0'][:]
		lon = ufile.variables['lon_0'][:]
		ufile.close()
	
		uwndnew = [] ;    ## because GEFS data runs from 90N lat to 90S lat
		for j in range(lonpts1deg) :
			uwndnew.append([])
			for i in range(latpts1deg) :
				uwndnew[j].append(0.) 
			
		for i in range(latpts1deg) :
			for j in range(lonpts1deg) : 
				uwndnew[j][(latpts1deg-1)-i] = uwnd[i][j] 
	
		ilat = 0     ### Grab data only over the N. Pacific Basin
		for i in range((latSind-1),latNind) :
			ilon = 0 
			for j in range((lonWind-1),lonEind) :
				upertens[ilon][ilat][t] = uwndnew[j][i] - meanwind1deg[j][i][mcount-1] ;
				ilon = ilon + 1 
			ilat = ilat + 1 
		tcount = tcount + 1 
		mcount = mcount + 1 
		
		if (mcount > montimes) :    ##### If advancing to next month update climo file ######
			newmonindex = monindex + 1 
			if (newmonindex > 12) :
				newmonindex = 1 
				year = year + 1 
			newmontimes = mtimes[newmonindex-1] 
			newmonth = mon[newmonindex-1] 
		
			filinstr = 'uwnd'+newmonth+'.out' 
			file = open(filinstr,"r") 
			data = file.readlines()
			
			meanpre = []
			meanwind = []
			meanwind1deg = []
			for j in range(lonpts) :
				meanpre.append([])
				meanwind.append([])
				for i in range(latpts) :
					meanpre[j].append([]) 
					meanwind[j].append([]) 
					for k in range(newmontimes) :
						meanpre[j][i].append(0.)
						meanwind[j][i].append(0.)
				
			for j in range(lonpts1deg) :
				meanwind1deg.append([])
				for i in range(latpts1deg) :
					meanwind1deg[j].append([])
					for k in range(newmontimes) :
						meanwind1deg[j][i].append(0.) 
						
			count = 0 ;
			for k in range(newmontimes) :
				for i in range(latpts) :
					for j in range(lonpts) :
						meanpre[j][(latpts-1)-i][k] = float(data[count]) 
						count = count + 1 
						
			if ((year%4) == 0) :
				if (newmonindex == 2) :
					newmontimes = 116 
					filinstr = 'uwnd0229.out' 
					file = open(filinstr,"r")
					data = file.readlines() 
			
					for j in range(lonpts) :
						for i in range(latpts) :
							for l in range(4) :
								meanpre[j][i].append(0.) 
								meanwind[j][i].append(0.)
						
					for j in range(lonpts1deg) :
						for i in range(latpts1deg) :
							for l in range(4) :
								meanwind1deg[j][i].append(0.) 
						
					count = 0 
					for l in range(4) :
						for i in range(latpts) :
							for j in range(lonpts) :
								meanpre[j][(latpts-1)-i][k+(l+1)] = float(data[count]) 
								count = count + 1 
			
			for k in range(newmontimes) :
				for i in range(latpts) :
					for j in range(360) :
						meanwind[j+360][i][k] = meanpre[j][i][k] 
					for j in range(360) :
						meanwind[j][i][k] = meanpre[j+360][i][k] 
	
			for k in range(newmontimes) :
				for i in range(latpts1deg) :
					for j in range(lonpts1deg) :
						meanwind1deg[j][i][k] = meanwind[(j*2)][(i*2)][k] 
			mcount = 1
			
	######################## Weight data by sqrt(cos(lat)) ###############################
	
	numdates = len(upertens[0][0]) 

	for k in range(numdates) :
		for i in range(steplat) :
			for j in range(steplon) :
				# weight by sqrt(cos(lat)) #
				upertens[j][i][k] = upertens[j][i][k]*sqrt(cos(Lat[i]*(pi/180.))) 
				
	########################## Two-dimensionalize data ##################################

	upert2Dens = []
	count = 0 
	for i in range(steplat) :
		for j in range(steplon) :
			upert2Dens.append([]) 
			for k in range(numdates) :
				upert2Dens[count].append(0.) 
			count = count + 1		

	for k in range(numdates) :
		count = 0 
		for i in range(steplat) :
			for j in range(steplon) :
				upert2Dens[count][k] = upertens[j][i][k] 
				count = count + 1 
				
	##################### Project Uwnd Perts onto EOF patterns #############################

	#### transpose upert2Dens ####
	upert2Dens_trans = []
	temppc1ens = []
	temppc2ens = []
	temppc3ens = []
	pc1ens = []
	pc2ens = []
	pc3ens = []
	for k in range(numdates) :
		upert2Dens_trans.append([]) 
		for i in range(len(upert2Dens)) :
			upert2Dens_trans[k].append(upert2Dens[i][k]) 
	
	for k in range(numdates) :
		sumpc1ens = 0 
		sumpc2ens = 0
		sumpc3ens = 0
		for i in range(len(eof12D)) :    ### upert2Dens'*eof12D - matrix multiplication ###
			sumpc1ens = sumpc1ens + (upert2Dens_trans[k][i]*eof12D[i])    ### Creates an 40x1 array of PCs ###
			sumpc2ens = sumpc2ens + (upert2Dens_trans[k][i]*eof22D[i])
			sumpc3ens = sumpc3ens + (upert2Dens_trans[k][i]*eof32D[i])
		temppc1ens.append(sumpc1ens)    
		temppc2ens.append(sumpc2ens)
		temppc3ens.append(sumpc3ens)
	
	for k in range(numdates) :
		temppc1ens[k] = temppc1ens[k]/(sqrt(1.3253e5)) 
		temppc2ens[k] = temppc2ens[k]/(sqrt(1.0022e5)) 
		temppc3ens[k] = temppc3ens[k]/(sqrt(0.8083e5)) 
	
	####### Add the last 4 PCs from GFS analysis data to the ensemble member PC list #########
		
	for k in range(4) :
		pc1ens.append(pc1[41+k])  
		pc2ens.append(pc2[41+k])  
		
	for k in range(40) :
		pc1ens.append(temppc1ens[k]) 
		pc2ens.append(temppc2ens[k]) 
	
	###### smooth PCs over the previous/next 1 day ####

	pc1enssmoo = []
	pc2enssmoo = []
	pcsens.append([])
	numdates = numdates + 4 ;
	count = 0 
	for k in range(4,(numdates-4)) :     ## 5:(numdates-4) for +/- 1 day ##
		weight = 5.0                     ## max weight = 5 for +/- 1 day ##
		pc1enssum = pc1ens[k]*weight 
		pc2enssum = pc2ens[k]*weight 
		for l in range(1,5) :
			weight = weight - 1.0 
			pc1enssum = pc1enssum + pc1ens[k-l]*weight + pc1ens[k+l]*weight 
			pc2enssum = pc2enssum + pc2ens[k-l]*weight + pc2ens[k+l]*weight 
		pc1enssmoo.append(pc1enssum/25.0)      ##  25.0 for sum of weights +/- 1 day ##
		pc2enssmoo.append(pc2enssum/25.0) 
		count = count + 1 
	
	for i in range(len(pc1enssmoo)) :
		pcsens[e].append([])
		pcsens[e][i].append(pc1enssmoo[i]) 
		pcsens[e][i].append(pc2enssmoo[i])

############### Figure 1 - Probabilistic forecast trajectory w/ analysis and ens. mean  ##################

plt.figure(1)

probgrid = []
probgridfinal = []
for j in range(161) :
	probgrid.append([])
	probgridfinal.append([])
	for i in range(161) :
		probgrid[j].append([])
		probgridfinal[j].append(0.)
		for e in range(21) :
			probgrid[j][i].append(0.)

xpts = linspace(-4,4,161) 
ypts = linspace(-4,4,161) 

for e in range(21) :   ## change this to 21 once done testing ##
	###### Compile the trajectory frequencies into one variable ######
	for t in range(36) :
		for i in range(161) :
			for j in range(161) : 
				xcoord = xpts[j] 
				ycoord = ypts[i] 
				xcoordpc = pcsens[e][t][0]
				ycoordpc = pcsens[e][t][1]
				if (sqrt((xcoord - xcoordpc)**2. + (ycoord - ycoordpc)**2.) <= 0.25) :
					if probgrid[i][j][e] == 1 :
						probgrid[i][j][e] = 1. 
					else :
						probgrid[i][j][e] = probgrid[i][j][e] + 1. 

for e in range(21) :
	for i in range(161) : 
		for j in range(161) :
			probgridfinal[i][j] = probgridfinal[i][j] + probgrid[i][j][e] 

for i in range(161) :
	for j in range(161) :
		probgridfinal[i][j] = probgridfinal[i][j]/(e+1) ;


X,Y = meshgrid(xpts,ypts)
cp = plt.contourf(X,Y,probgridfinal,[0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0])
cbar = colorbar(ticks=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0],drawedges='True')
cbar.set_label('Probability')

######### Calculate ensemble mean #######

pc1ensmean = []
pc2ensmean = []
for i in range(36) :
	pc1enssum = 0.
	pc2enssum = 0.
	for e in range(21) :  ## change to 21 when run for all ens. mems. ##
		pc1enssum = pc1enssum + pcsens[e][i][0]
		pc2enssum = pc2enssum + pcsens[e][i][1] 
	pc1ensmean.append(pc1enssum/(e+1))
	pc2ensmean.append(pc2enssum/(e+1)) 

pcsensmean = []
for i in range(36) :
	pcsensmean.append([]) 
	pcsensmean[i].append(pc1ensmean[i])
	pcsensmean[i].append(pc2ensmean[i])

########## plot ensemble mean ###########

count00z = 1 
for i in range(36) :
	if (i == (4*count00z-1)) :
		if i == 35 :
			s1 = plt.plot(pcsensmean[i][0],pcsensmean[i][1],'wd',markersize=10,markerfacecolor='w',markeredgewidth=2,markeredgecolor='k')
			plt.text((pcsensmean[i][0]+0.10),(pcsensmean[i][1]+0.10),str(int(labelfcast[count00z])),fontsize=8,fontweight='bold',color='k')
	  		count00z = count00z + 1 
		else :
 			s2 = plot(pcsensmean[i][0],pcsensmean[i][1],'wd',markersize=10,markerfacecolor='w',markeredgewidth=2,markeredgecolor='k') 
 			plt.text((pcsensmean[i][0]+0.10),(pcsensmean[i][1]+0.10),str(int(labelfcast[count00z])),fontsize=8,fontweight='bold',color='k')
	  		count00z = count00z + 1 
	else :
		s3 = plt.plot(pcsensmean[i][0],pcsensmean[i][1],'wo',markersize=5,markerfacecolor='w',markeredgewidth=1.5,markeredgecolor='k')
	
######### overlay 00-h analysis ##########

s4 = plt.plot(pcs[40][0],pcs[40][1],'gd',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')
plt.text((pcs[40][0]+0.10),(pcs[40][1]+0.10),str(int(labelfcast[0])),fontsize=8,fontweight='bold',color='k')

gefs = mlines.Line2D([], [], color='w', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='GEFS Mean',linestyle='None')
initial = mlines.Line2D([], [], color='g', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Initial.',linestyle='None')	
plt.legend(handles=[gefs,initial],loc=1,fontsize=10,numpoints=1)

######## Labels and axes for the figure #######

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('NPJ Phase Diagram initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('probdiagram.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/probdiagram.pdf') 
#plt.show()

################ Figure 2 - GFS analysis, GFS forecast, GEFS mean trajectory  %%%%%%%%%%%%%%%

plt.figure(2)

count00z = 1 
countfcast = 1 
for i in range(77) :   
	if i < 40 :  ## For GFS analysis times ##
		if i == (4*count00z - 4) :
			s1 = plt.plot(pcs[i][0],pcs[i][1],'kd',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')
			plt.text((pcs[i][0]+0.10),(pcs[i][1]+0.10),str(int(labelanl[count00z])),fontsize=8,fontweight='bold')
			count00z = count00z + 1 
		else :
			s2 = plt.plot(pcs[i][0],pcs[i][1],'ko',markersize=5,markerfacecolor='k',markeredgewidth=1.5,markeredgecolor='k')
	else :     ## For GFS forecast times ##
		if i == (4*count00z - 4) :
			if i == 40 :
				s3 = plt.plot(pcs[i][0],pcs[i][1],'gd',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')
				plt.text((pcs[i][0]+0.10),(pcs[i][1]+0.10),str(int(labelanl[count00z])),fontsize=8,fontweight='bold')
				count00z = count00z + 1 
			else :
				s4 = plt.plot(pcs[i][0],pcs[i][1],'rd',markersize=10,markerfacecolor='r',markeredgewidth=2,markeredgecolor='k')
				plt.text((pcs[i][0]+0.10),(pcs[i][1]+0.10),str(int(labelfcast[countfcast])),fontsize=8,fontweight='bold')
				count00z = count00z + 1 
				countfcast = countfcast + 1 
		else :
			s5 = plt.plot(pcs[i][0],pcs[i][1],'ro',markersize=5,markerfacecolor='r',markeredgewidth=1.5,markeredgecolor='k')

count00z = 1 
for i in range(36) :
	if (i == (4*count00z-1)) :
		if i == 35 :
			s6 = plt.plot(pcsensmean[i][0],pcsensmean[i][1],'bd',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')
			plt.text((pcsensmean[i][0]+0.10),(pcsensmean[i][1]+0.10),str(int(labelfcast[count00z])),fontsize=8,fontweight='bold')
	  		count00z = count00z + 1 
		else :
 			s7 = plot(pcsensmean[i][0],pcsensmean[i][1],'bd',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k') 
 			plt.text((pcsensmean[i][0]+0.10),(pcsensmean[i][1]+0.10),str(int(labelfcast[count00z])),fontsize=8,fontweight='bold')
	  		count00z = count00z + 1 
	else :
		s8 = plt.plot(pcsensmean[i][0],pcsensmean[i][1],'bo',markersize=5,markerfacecolor='b',markeredgewidth=1.5,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='k', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='GFS Analysis',linestyle='None')
gfs = mlines.Line2D([], [], color='r', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='GFS Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='b', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='GEFS Mean',linestyle='None')
initial = mlines.Line2D([], [], color='g', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Initial.',linestyle='None')	
plt.legend(handles=[analysis,gfs,gefs,initial],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('NPJ Phase Diagram initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('deterministic.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/deterministic.pdf') 
#plt.show()

################## Figure 3 - 1-day Ensemble Scatterplot Forecast  ######################

plt.figure(3)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][3][0],pcsens[e][3][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[3][0],pcsensmean[3][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('1-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D1ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D1ens.pdf') 
#plt.show()

################## Figure 4 - 2-day Ensemble Scatterplot Forecast  ######################

plt.figure(4)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][7][0],pcsens[e][7][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[7][0],pcsensmean[7][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('2-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D2ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D2ens.pdf') 
#plt.show()

################## Figure 5 - 3-day Ensemble Scatterplot Forecast  ######################

plt.figure(5)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][11][0],pcsens[e][11][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[11][0],pcsensmean[11][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('3-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D3ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D3ens.pdf') 
#plt.show()

################## Figure 6 - 4-day Ensemble Scatterplot Forecast  ######################

plt.figure(6)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][15][0],pcsens[e][15][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[15][0],pcsensmean[15][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('4-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D4ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D4ens.pdf') 
#plt.show()

################## Figure 7 - 5-day Ensemble Scatterplot Forecast  ######################

plt.figure(7)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][19][0],pcsens[e][19][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[19][0],pcsensmean[19][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('5-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D5ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D5ens.pdf') 
#plt.show()

################## Figure 8 - 6-day Ensemble Scatterplot Forecast  ######################

plt.figure(8)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][23][0],pcsens[e][23][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[23][0],pcsensmean[23][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('6-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D6ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D6ens.pdf') 
#plt.show()

################## Figure 9 - 7-day Ensemble Scatterplot Forecast  ######################

plt.figure(9)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][27][0],pcsens[e][27][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[27][0],pcsensmean[27][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('7-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D7ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D7ens.pdf') 
#plt.show()

################## Figure 10 - 8-day Ensemble Scatterplot Forecast  ######################

plt.figure(10)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][31][0],pcsens[e][31][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[31][0],pcsensmean[31][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('8-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D8ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D8ens.pdf') 
#plt.show()

################## Figure 11 - 9-day Ensemble Scatterplot Forecast  ######################

plt.figure(11)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

for e in range(21) :      
	s2 = plt.plot(pcsens[e][35][0],pcsens[e][35][1],'rx',markersize=8,markeredgewidth=1.5)
	s3 = plt.plot(pcsensmean[35][0],pcsensmean[35][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('9-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D9ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D9ens.pdf') 
#plt.show()

####################### Figure 12 - 0-day Forecast  ########################
	
plt.figure(12)

s1 = plt.plot(pcs[40][0],pcs[40][1],'go',markersize=10,markerfacecolor='g',markeredgewidth=2,markeredgecolor='k')

analysis = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='0-h Forecast',linestyle='None')	
gefs = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,markeredgecolor='r',label='Ens. Mems.',linestyle='None')
mean = mlines.Line2D([], [], color='b', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Ens. Mean',linestyle='None')	
plt.legend(handles=[analysis,gefs,mean],loc=1,fontsize=10,numpoints=1)

xaxe = plt.plot([-4,4],[-4,4],'k',linewidth=3) 
yaxe = plt.plot([-4,4],[4,-4],'k',linewidth=3) 

axes = plt.axis([-4, 4, -4, 4])
xlab = plt.xlabel('PC 1',fontsize=12,fontweight='bold')
ylab = plt.ylabel('PC 2',fontsize=12,fontweight='bold') 
plt.text(-3.75,0,'Jet Retraction',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(3.75,0,'Jet Extension',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=90)
plt.text(0,-3.75,'Equatorward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
plt.text(0,3.75,'Poleward Shift',fontsize=14,fontweight='bold',color='r',horizontalalignment='center',verticalalignment='center',rotation=0)
titl = plt.title('0-Day Forecast initialized 0000 UTC '+daystr+' '+month+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('D0ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/D0ens.pdf') 
#plt.show()

############################# Write Outfiles  ###################################

pcsfile = '/winterslab_rit/realtime/verification/gfspcs_'+monstr+daystr+'.out'
pcsensfile = '/winterslab_rit/realtime/verification/gefspcs_'+monstr+daystr+'.out'
pcsensmeanfile = '/winterslab_rit/realtime/verification/gefsmeanpcs_'+monstr+daystr+'.out'

file_pcs = open(pcsfile,'w')
for i in range(77) :
	file_pcs.write(str(pcs[i][0])+','+str(pcs[i][1])+'\n')
file_pcs.close()

file_mean = open(pcsensmeanfile,'w') 
for i in range(36) :
	file_mean.write(str(pcsensmean[i][0])+','+str(pcsensmean[i][1])+'\n') 
file_mean.close()

file_gefs = open(pcsensfile,'w') 
for e in range(21) :
	for i in range(36) :
		file_gefs.write(str(pcsens[e][i][0])+','+str(pcsens[e][i][1])+'\n')
file_gefs.close()