#!/usr/bin/env python

###############################################################################
## This script reads 9-day old PCs and verifies them against the most recent ##
## analysis  																 ##
###############################################################################

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 current time to grab correct PC files #################

hr = 0 
year = 2017 
monindex = 6
dy = 9
#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]  
dyold = dy - 9      ### to go back 11 days ###

labelanl = linspace(dyold,dy,10)   ### create labels for plots ###

if (dyold <= 0) :     ### Adjust the dates if dyold <= 0 ### 
	monindex = monindex - 1 
	if (monindex == 0) :
		monindex = 12 
		year = year - 1 
	if ((year%4) == 0) :
		if (monindex == 2) :
			dyold = dyold + mondaysleap[monindex-1] 
		else :
			dyold = dyold + mondaysleap[monindex-1]
	else :
		dyold = dyold + 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] 

monoldstr = monnum[monindex-1] 
dayoldstr = day[dyold-1]  

print ('Verify Year: %d' % year)
print ('Verify Month: %d' % monindex)
print ('Verify Start Day: %d' % dyold)
print labelanl

######################## Read in PC text files ###############################

### old files ###
filin = '/winterslab_rit/realtime/verification/gfspcs_'+monoldstr+dayoldstr+'.out'
filin2 = '/winterslab_rit/realtime/verification/gefspcs_'+monoldstr+dayoldstr+'.out'
filin3 = '/winterslab_rit/realtime/verification/gefsmeanpcs_'+monoldstr+dayoldstr+'.out'

### today files ###
filin4 = '/winterslab_rit/realtime/verification/gfspcs_'+monstr+daystr+'.out'
filin5 = '/winterslab_rit/realtime/verification/gefspcs_'+monstr+daystr+'.out'
filin6 = '/winterslab_rit/realtime/verification/gefsmeanpcs_'+monstr+daystr+'.out'

### read in data ###

pcsold = []
count = 0
file = open(filin, "r")
totlines = file.readlines()
for line in totlines :
	pcsold.append([])
	words = line.split('  ')
	pc1 = words[1]
	pc2 = words[2]
	pcsold[count].append(float(pc1)) 
	pcsold[count].append(float(pc2))
	count = count + 1 

rawgefsold = []
count = 0
file = open(filin2, "r")
totlines = file.readlines()
for line in totlines :
	rawgefsold.append([])
	words = line.split('  ')
	pc1 = words[1]
	pc2 = words[2]
	rawgefsold[count].append(float(pc1)) 
	rawgefsold[count].append(float(pc2))
	count = count + 1 
	
pcsensmeanold = []
count = 0
file = open(filin3, "r")
totlines = file.readlines()
for line in totlines :
	pcsensmeanold.append([])
	words = line.split('  ')
	pc1 = words[1]
	pc2 = words[2]
	pcsensmeanold[count].append(float(pc1)) 
	pcsensmeanold[count].append(float(pc2))
	count = count + 1 
	
pcs = []
count = 0
file = open(filin4, "r")
totlines = file.readlines()
for line in totlines :
	pcs.append([])
	words = line.split('  ')
	pc1 = words[1]
	pc2 = words[2]
	pcs[count].append(float(pc1)) 
	pcs[count].append(float(pc2))
	count = count + 1 
	
rawgefs = []
count = 0
file = open(filin5, "r")
totlines = file.readlines()
for line in totlines :
	rawgefs.append([])
	words = line.split('  ')
	pc1 = words[1]
	pc2 = words[2]
	rawgefs[count].append(float(pc1)) 
	rawgefs[count].append(float(pc2))
	count = count + 1
	
pcsensmean = []
count = 0
file = open(filin6, "r")
totlines = file.readlines()
for line in totlines :
	pcsensmean.append([])
	words = line.split('  ')
	pc1 = words[1]
	pc2 = words[2]
	pcsensmean[count].append(float(pc1)) 
	pcsensmean[count].append(float(pc2))
	count = count + 1

pcsens = []
pcsensold = []
count = 0 
for e in range(21) :
	pcsens.append([])
	pcsensold.append([])
	for i in range(36) :
		pcsens[e].append([])
		pcsensold[e].append([])
		pcsens[e][i].append(rawgefs[count][0])
		pcsens[e][i].append(rawgefs[count][1])
		pcsensold[e][i].append(rawgefsold[count][0])
		pcsensold[e][i].append(rawgefsold[count][1]) 
		count = count + 1 
		 
############# Figure 1 - Old Probabilistic and ens. mean forecast trajectory w/ analysis  #####################

plt.figure(1)

probgrid = []
probgridfinal = []
probgridtime = []
for j in range(161) :
	probgrid.append([])
	probgridfinal.append([])
	probgridtime.append([])
	for i in range(161) :
		probgrid[j].append([])
		probgridfinal[j].append(0.)
		probgridtime[j].append([])
		for e in range(21) :
			probgrid[j][i].append(0.)
		for t in range(36) :
			probgridtime[j][i].append(0.)
			
xpts = linspace(-4,4,161) 
ypts = linspace(-4,4,161) 

##### Probability of the analysis being in a certain part of the phase diagram at any time in the forecast  ######

for e in range(21) :   
	###### Compile the trajectory frequencies into one variable (old) ######
	for t in range(36) :
		for i in range(161) :
			for j in range(161) : 
				xcoord = xpts[j] 
				ycoord = ypts[i] 
				xcoordpc = pcsensold[e][t][0]
				ycoordpc = pcsensold[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) ;

########## For individual forecast time probabilities ###########

for t in range(36) :
	for e in range(21) :
		for i in range(161) :
			for j in range(161) :
				xcoord = xpts[j] 
				ycoord = ypts[i] 
				xcoordpc = pcsensold[e][t][0]
				ycoordpc = pcsensold[e][t][1] 
				if (sqrt( (xcoord - xcoordpc)**2. + (ycoord - ycoordpc)**2. ) <= 0.25) :
					probgridtime[i][j][t] = probgridtime[i][j][t] + 1 
					
for t in range(36) :
	for i in range(161) :
		for j in range(161) :
			probgridtime[i][j][t] = probgridtime[i][j][t]/21. 

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')

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

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

############ overlay analysis trajectory (new) ############

count00z = 2 
for i in range(4,41) :   ## for D-9 to D0 ##
	if (i == 4) :
		s5 = 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[0])),fontsize=8,fontweight='bold',color='k')
		count00z = count00z + 1 
	elif i == (4*count00z - 4) :
		s3 = plt.plot(pcs[i][0],pcs[i][1],'yd',markersize=10,markerfacecolor='y',markeredgewidth=2,markeredgecolor='k')
		plt.text((pcs[i][0]+0.10),(pcs[i][1]+0.10),str(int(labelanl[count00z-2])),fontsize=8,fontweight='bold',color='k')
		count00z = count00z + 1 
	else :
		s4 = plt.plot(pcs[i][0],pcs[i][1],'yo',markersize=5,markerfacecolor='y',markeredgewidth=1.5,markeredgecolor='k')

gefs = mlines.Line2D([], [], color='w', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='GEFS Mean',linestyle='None')
gfs = mlines.Line2D([], [], color='y', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='GFS Analysis',linestyle='None')	
initial = mlines.Line2D([], [], color='g', marker='d',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Initial.',linestyle='None')	
plt.legend(handles=[gefs,gfs,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('Verified NPJ Phase Diagram from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_probdiagram.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_probdiagram.pdf') 
#plt.show()

############### Figure 2 - Old GFS forecast, old GEFS ens mean and GFS analysis  #################

plt.figure(2)

count00z = 12 
for i in range(41,77) :   ## For F06-F216 in old GFS forecast
	if i == (4*count00z - 4) :
		s4 = plt.plot(pcsold[i][0],pcsold[i][1],'rd',markersize=10,markerfacecolor='r',markeredgewidth=2,markeredgecolor='k')
		plt.text((pcsold[i][0]+0.10),(pcsold[i][1]+0.10),str(int(labelanl[count00z-11])),fontsize=8,fontweight='bold',color='k')
		count00z = count00z + 1 
	else :
		s5 = plt.plot(pcsold[i][0],pcsold[i][1],'ro',markersize=5,markerfacecolor='r',markeredgewidth=1.5,markeredgecolor='k')

count00z = 1 
for i in range(36) :    ## For F06-F216 in old GEFS forecast
	if (i == (4*count00z-1)) :
		if i == 35 :
			s1 = plt.plot(pcsensmeanold[i][0],pcsensmeanold[i][1],'bd',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')
			plt.text((pcsensmeanold[i][0]+0.10),(pcsensmeanold[i][1]+0.10),str(int(labelanl[count00z])),fontsize=8,fontweight='bold',color='k')
	  		count00z = count00z + 1 
		else :
 			s2 = plot(pcsensmeanold[i][0],pcsensmeanold[i][1],'bd',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k') 
 			plt.text((pcsensmeanold[i][0]+0.10),(pcsensmeanold[i][1]+0.10),str(int(labelanl[count00z])),fontsize=8,fontweight='bold',color='k')
	  		count00z = count00z + 1 
	else :
		s3 = plt.plot(pcsensmeanold[i][0],pcsensmeanold[i][1],'bo',markersize=5,markerfacecolor='b',markeredgewidth=1.5,markeredgecolor='k')

count00z = 2 
for i in range(4,41) :   ## for D-9 to D0 in GFS analysis ##
	if (i == 4) :
		s5 = 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[0])),fontsize=8,fontweight='bold',color='k')
		count00z = count00z + 1 
	elif i == (4*count00z - 4) :
		s3 = 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-2])),fontsize=8,fontweight='bold',color='k')
		count00z = count00z + 1 
	else :
		s4 = plt.plot(pcs[i][0],pcs[i][1],'ko',markersize=5,markerfacecolor='k',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('Verified NPJ Phase Diagram from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_deterministic.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_deterministic.pdf') 
#plt.show()

######################## Figure 3 - Verified 1-day Forecast ###################################

plt.figure(3)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][3][0],pcsensold[e][3][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[3][0],pcsensmeanold[3][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[8][0],pcs[8][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 1 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 1-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D1ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D1ens.pdf') 
#plt.show()

######################## Figure 4 - Verified 2-day Forecast ###################################

plt.figure(4)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][7][0],pcsensold[e][7][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[7][0],pcsensmeanold[7][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[12][0],pcs[12][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 2 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 2-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D2ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D2ens.pdf') 
#plt.show()

######################## Figure 5 - Verified 3-day Forecast ###################################

plt.figure(5)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][11][0],pcsensold[e][11][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[11][0],pcsensmeanold[11][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[16][0],pcs[16][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 3 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 3-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D3ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D3ens.pdf') 
#plt.show()

######################## Figure 6 - Verified 4-day Forecast ###################################

plt.figure(6)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][15][0],pcsensold[e][15][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[15][0],pcsensmeanold[15][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[20][0],pcs[20][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 4 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 4-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D4ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D4ens.pdf') 
#plt.show()

######################## Figure 7 - Verified 5-day Forecast ###################################

plt.figure(7)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][19][0],pcsensold[e][19][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[19][0],pcsensmeanold[19][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[24][0],pcs[24][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 5 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 5-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D5ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D5ens.pdf') 
#plt.show()

######################## Figure 8 - Verified 6-day Forecast ###################################

plt.figure(8)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][23][0],pcsensold[e][23][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[23][0],pcsensmeanold[23][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[28][0],pcs[28][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 6 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 6-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D6ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D6ens.pdf') 
#plt.show()

######################## Figure 9 - Verified 7-day Forecast ###################################

plt.figure(9)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][27][0],pcsensold[e][27][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[27][0],pcsensmeanold[27][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[32][0],pcs[32][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 7 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 7-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D7ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D7ens.pdf') 
#plt.show()

######################## Figure 10 - Verified 8-day Forecast ###################################

plt.figure(10)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][31][0],pcsensold[e][31][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[31][0],pcsensmeanold[31][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

s2 = plt.plot(pcs[36][0],pcs[36][1],'ko',markersize=10,markerfacecolor='k',markeredgewidth=2,markeredgecolor='k')

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 8 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 8-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D8ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D8ens.pdf') 
#plt.show()

######################## Figure 11 - Verified 9-day Forecast ###################################

plt.figure(11)

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

for e in range(21) :    
	s3 = plt.plot(pcsensold[e][35][0],pcsensold[e][35][1],'rx',markersize=8,markeredgewidth=1.5) 

s4 = plt.plot(pcsensmeanold[35][0],pcsensmeanold[35][1],'bo',markersize=10,markerfacecolor='b',markeredgewidth=2,markeredgecolor='k')

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

D0 = mlines.Line2D([], [], color='g', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 0 Forecast',linestyle='None')
D9 = mlines.Line2D([], [], color='k', marker='o',markersize=8,markeredgewidth=1.5,markeredgecolor='k',label='Day 9 Forecast',linestyle='None')
mems = mlines.Line2D([], [], color='r', marker='x',markersize=8,markeredgewidth=1.5,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=[D0,D9,mems,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('Verified 9-Day Forecast from 0000 UTC '+dayoldstr+' '+mon[monindex-1]+' '+str(year),fontsize=14,fontweight='bold')
plt.grid(True)
#plt.savefig('verif_D9ens.pdf')
plt.savefig('/winterslab_rit/realtime/images/phase_diagrams/verif_D9ens.pdf') 
#plt.show()

########################### Calculate Verification Metrics ##################################

	########### Determine what quadrant of the diagram the analysis is in ##############
	
JetMode = 0 ;     ## 1 = extension, 2 = retraction, 3 = poleward shift, 4 = equatorward shift, 5 = origin
if (pcs[4][0] > 0 and abs(pcs[4][1]) < pcs[4][0] and (pcs[4][0]**2. + pcs[4][1]**2.) >= 1.00 ) :   ### Jet Extension
	JetMode = 1 

if (pcs[4][0] < 0 and abs(pcs[4][1]) < abs(pcs[4][0]) and (pcs[4][0]**2. + pcs[4][1]**2.) >= 1.00 ) :  ### Jet Retraction
	JetMode = 2 
	
if (pcs[4][1] > 0 and abs(pcs[4][0]) < pcs[4][1] and (pcs[4][0]**2. + pcs[4][1]**2.) >= 1.00 ) :    ### Poleward Shift
	JetMode = 3
	
if (pcs[4][1] < 0 and abs(pcs[4][0]) < abs(pcs[4][1]) and (pcs[4][0]**2. + pcs[4][1]**2.) >= 1.00 ) :   ### Equatorward Shift
	JetMode = 4 
	
if ((pcs[4][0]**2. + pcs[4][1]**2.) < 1.00) :    ### Origin
	JetMode = 5
	
  ##################### Reliability diagram info #########################

fcastct = []    ## bins for probabilities
fcasthit = []   ## bins for hits by probabilities
for i in range(9) :
	fcastct.append([])
	fcasthit.append([])
	for j in range(13) :
		fcastct[i].append(0.)
		fcasthit[i].append(0.)
		
countday = 1 
countanl = 2 

for t in range(9) :
	fcastct[t][0] = JetMode
	fcasthit[t][0] = JetMode
	for i in range(161) :
		for j in range(161) :
			xcoord = xpts[j] 
			ycoord = ypts[i] 
			xcoordanl = pcs[4*countanl][0] 
			ycoordanl = pcs[4*countanl][1]
			if (probgridtime[i][j][4*countday-1] == 0) :     ## 0% probability
				fcastct[t][1] = fcastct[t][1] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][1] = fcasthit[t][1] + 1
			if (probgridtime[i][j][4*countday-1] > 0 and probgridtime[i][j][4*countday-1] <= 0.10 ) :   ## 0-10% probability
				fcastct[t][2] = fcastct[t][2] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][2] = fcasthit[t][2] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.10 and probgridtime[i][j][4*countday-1] <= 0.20 ) :   ## 10-20% probability
				fcastct[t][3] = fcastct[t][3] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][3] = fcasthit[t][3] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.20 and probgridtime[i][j][4*countday-1] <= 0.30 ) :   ## 20-30% probability
				fcastct[t][4] = fcastct[t][4] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][4] = fcasthit[t][4] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.30 and probgridtime[i][j][4*countday-1] <= 0.40 ) :   ## 30-40% probability
				fcastct[t][5] = fcastct[t][5] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][5] = fcasthit[t][5] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.40 and probgridtime[i][j][4*countday-1] <= 0.50 ) :   ## 40-50% probability
				fcastct[t][6] = fcastct[t][6] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][6] = fcasthit[t][6] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.50 and probgridtime[i][j][4*countday-1] <= 0.60 ) :   ## 50-60% probability
				fcastct[t][7] = fcastct[t][7] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][7] = fcasthit[t][7] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.60 and probgridtime[i][j][4*countday-1] <= 0.70 ) :   ## 60-70% probability
				fcastct[t][8] = fcastct[t][8] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][8] = fcasthit[t][8] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.70 and probgridtime[i][j][4*countday-1] <= 0.80 ) :   ## 70-80% probability
				fcastct[t][9] = fcastct[t][9] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][9] = fcasthit[t][9] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.80 and probgridtime[i][j][4*countday-1] <= 0.90 ) :   ## 80-90% probability
				fcastct[t][10] = fcastct[t][10] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][10] = fcasthit[t][10] + 1 
			if (probgridtime[i][j][4*countday-1] > 0.90 and probgridtime[i][j][4*countday-1] < 1.00 ) :   ## 90-100% probability
				fcastct[t][11] = fcastct[t][11] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][11] = fcasthit[t][11] + 1 
			if (probgridtime[i][j][4*countday-1] == 1.00 ):   ## 100% probability
				fcastct[t][12] = fcastct[t][12] + 1 
				if (sqrt( (xcoord - xcoordanl)**2. + (ycoord - ycoordanl)**2. ) <= 0.25) :
					fcasthit[t][12] = fcasthit[t][12] + 1 
	countday = countday + 1 
	countanl = countanl + 1 
	
	
	###################### Probability of detection (hit/miss) ##########################

hitmiss = []
for t in range(10) :
	hitmiss.append(0.)
hitmiss[0] = JetMode 

for t in range(9) :
	for p in range(11) :
		if (fcasthit[t][p+2] > 0) :
			hitmiss[t+1] = 1

	######################### Day 1 forecast error (PC units) ###############################

D1dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][3][0] - pcs[8][0])**2. + (pcsensold[e][3][1] - pcs[8][1])**2. )
	D1dist = D1dist + tempsum 
	
aveD1dist = D1dist/21. 

D1gfsdist = sqrt( (pcsold[44][0] - pcs[8][0])**2. + (pcsold[44][1] - pcs[8][1])**2. )
D1gefsdist = sqrt( (pcsensmeanold[3][0] - pcs[8][0])**2. + (pcsensmeanold[3][1] - pcs[8][1])**2. )

	######################### Day 2 forecast error (PC units) ###############################

D2dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][7][0] - pcs[12][0])**2. + (pcsensold[e][7][1] - pcs[12][1])**2. )
	D2dist = D2dist + tempsum 
	
aveD2dist = D2dist/21. 

D2gfsdist = sqrt( (pcsold[48][0] - pcs[12][0])**2. + (pcsold[48][1] - pcs[12][1])**2. )
D2gefsdist = sqrt( (pcsensmeanold[7][0] - pcs[12][0])**2. + (pcsensmeanold[7][1] - pcs[12][1])**2. )

	######################### Day 3 forecast error (PC units) ###############################

D3dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][11][0] - pcs[16][0])**2. + (pcsensold[e][11][1] - pcs[16][1])**2. )
	D3dist = D3dist + tempsum 
	
aveD3dist = D3dist/21. 

D3gfsdist = sqrt( (pcsold[52][0] - pcs[16][0])**2. + (pcsold[52][1] - pcs[16][1])**2. )
D3gefsdist = sqrt( (pcsensmeanold[11][0] - pcs[16][0])**2. + (pcsensmeanold[11][1] - pcs[16][1])**2. )

	######################### Day 4 forecast error (PC units) ###############################

D4dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][15][0] - pcs[20][0])**2. + (pcsensold[e][15][1] - pcs[20][1])**2. )
	D4dist = D4dist + tempsum 
	
aveD4dist = D4dist/21. 

D4gfsdist = sqrt( (pcsold[56][0] - pcs[20][0])**2. + (pcsold[56][1] - pcs[20][1])**2. )
D4gefsdist = sqrt( (pcsensmeanold[15][0] - pcs[20][0])**2. + (pcsensmeanold[15][1] - pcs[20][1])**2. )

	######################### Day 5 forecast error (PC units) ###############################

D5dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][19][0] - pcs[24][0])**2. + (pcsensold[e][19][1] - pcs[24][1])**2. )
	D5dist = D5dist + tempsum 
	
aveD5dist = D5dist/21. 

D5gfsdist = sqrt( (pcsold[60][0] - pcs[24][0])**2. + (pcsold[60][1] - pcs[24][1])**2. )
D5gefsdist = sqrt( (pcsensmeanold[19][0] - pcs[24][0])**2. + (pcsensmeanold[19][1] - pcs[24][1])**2. )

	######################### Day 6 forecast error (PC units) ###############################

D6dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][23][0] - pcs[28][0])**2. + (pcsensold[e][23][1] - pcs[28][1])**2. )
	D6dist = D6dist + tempsum 
	
aveD6dist = D6dist/21. 

D6gfsdist = sqrt( (pcsold[64][0] - pcs[28][0])**2. + (pcsold[64][1] - pcs[28][1])**2. )
D6gefsdist = sqrt( (pcsensmeanold[23][0] - pcs[28][0])**2. + (pcsensmeanold[23][1] - pcs[28][1])**2. )

	######################### Day 7 forecast error (PC units) ###############################

D7dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][27][0] - pcs[32][0])**2. + (pcsensold[e][27][1] - pcs[32][1])**2. )
	D7dist = D7dist + tempsum 
	
aveD7dist = D7dist/21. 

D7gfsdist = sqrt( (pcsold[68][0] - pcs[32][0])**2. + (pcsold[68][1] - pcs[32][1])**2. )
D7gefsdist = sqrt( (pcsensmeanold[27][0] - pcs[32][0])**2. + (pcsensmeanold[27][1] - pcs[32][1])**2. )

	######################### Day 8 forecast error (PC units) ###############################

D8dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][31][0] - pcs[36][0])**2. + (pcsensold[e][31][1] - pcs[36][1])**2. )
	D8dist = D8dist + tempsum 
	
aveD8dist = D8dist/21. 

D8gfsdist = sqrt( (pcsold[72][0] - pcs[36][0])**2. + (pcsold[72][1] - pcs[36][1])**2. )
D8gefsdist = sqrt( (pcsensmeanold[31][0] - pcs[36][0])**2. + (pcsensmeanold[31][1] - pcs[36][1])**2. )

	######################### Day 9 forecast error (PC units) ###############################

D9dist = 0 
for e in range(21) :
	tempsum = sqrt( (pcsensold[e][35][0] - pcs[40][0])**2. + (pcsensold[e][35][1] - pcs[40][1])**2. )
	D9dist = D9dist + tempsum 
	
aveD9dist = D9dist/21. 

D9gfsdist = sqrt( (pcsold[76][0] - pcs[40][0])**2. + (pcsold[76][1] - pcs[40][1])**2. )
D9gefsdist = sqrt( (pcsensmeanold[35][0] - pcs[40][0])**2. + (pcsensmeanold[35][1] - pcs[40][1])**2. )

	################### False Alarm Detection - Can we get the phase transitions right?? ##################

for t in range(9) :
	gfsfcastregime = [0,0,0,0,0,0]
	gefsfcastregime = [0,0,0,0,0,0]
	verifregime = [0,0,0,0,0,0]

	gfsfcastregime[0] = JetMode 
	gefsfcastregime[0] = JetMode
	verifregime[0] = JetMode  
	
	###### Forecasted Jet Phase in GEFS Mean Forecast ########
	if (pcsensmeanold[4*t+3][0] > 0 and abs(pcsensmeanold[4*t+3][1]) < pcsensmeanold[4*t+3][0] and (pcsensmeanold[4*t+3][0]**2. + pcsensmeanold[4*t+3][1]**2.) >= 1.00 ) :   ### Jet Extension
		gefsJetRegime_fcast = 1 
	
	if (pcsensmeanold[4*t+3][0] < 0 and abs(pcsensmeanold[4*t+3][1]) < abs(pcsensmeanold[4*t+3][0]) and (pcsensmeanold[4*t+3][0]**2. + pcsensmeanold[4*t+3][1]**2.) >= 1.00 ) :  ### Jet Retraction
		gefsJetRegime_fcast = 2 
	
	if (pcsensmeanold[4*t+3][1] > 0 and abs(pcsensmeanold[4*t+3][0]) < pcsensmeanold[4*t+3][1] and (pcsensmeanold[4*t+3][0]**2. + pcsensmeanold[4*t+3][1]**2.) >= 1.00 ) :    ### Poleward Shift
		gefsJetRegime_fcast = 3
	
	if (pcsensmeanold[4*t+3][1] < 0 and abs(pcsensmeanold[4*t+3][0]) < abs(pcsensmeanold[4*t+3][1]) and (pcsensmeanold[4*t+3][0]**2. + pcsensmeanold[4*t+3][1]**2.) >= 1.00 ) :   ### Equatorward Shift
		gefsJetRegime_fcast = 4 
	
	if ((pcsensmeanold[4*t+3][0]**2. + pcsensmeanold[4*t+3][1]**2.) < 1.00) :    ### Origin
		gefsJetRegime_fcast = 5
	
	###### Forecasted Jet Phase in GFS Deterministic Forecast ########
	if (pcsold[4*t+44][0] > 0 and abs(pcsold[4*t+44][1]) < pcsold[4*t+44][0] and (pcsold[4*t+44][0]**2. + pcsold[4*t+44][1]**2.) >= 1.00 ) :   ### Jet Extension
		gfsJetRegime_fcast = 1 
	
	if (pcsold[4*t+44][0] < 0 and abs(pcsold[4*t+44][1]) < abs(pcsold[4*t+44][0]) and (pcsold[4*t+44][0]**2. + pcsold[4*t+44][1]**2.) >= 1.00 ) :  ### Jet Retraction
		gfsJetRegime_fcast = 2 
	
	if (pcsold[4*t+44][1] > 0 and abs(pcsold[4*t+44][0]) < pcsold[4*t+44][1] and (pcsold[4*t+44][0]**2. + pcsold[4*t+44][1]**2.) >= 1.00 ) :    ### Poleward Shift
		gfsJetRegime_fcast = 3
	
	if (pcsold[4*t+44][1] < 0 and abs(pcsold[4*t+44][0]) < abs(pcsold[4*t+44][1]) and (pcsold[4*t+44][0]**2. + pcsold[4*t+44][1]**2.) >= 1.00 ) :   ### Equatorward Shift
		gfsJetRegime_fcast = 4 
	
	if ((pcsold[4*t+44][0]**2. + pcsold[4*t+44][1]**2.) < 1.00) :    ### Origin
		gfsJetRegime_fcast = 5
		
	###### Verified Jet Phase in GFS Analysis #########
	
	if (pcs[4*t+8][0] > 0 and abs(pcs[4*t+8][1]) < pcs[4*t+8][0] and (pcs[4*t+8][0]**2. + pcs[4*t+8][1]**2.) >= 1.00 ) :   ### Jet Extension
		JetRegime_verif = 1 
	
	if (pcs[4*t+8][0] < 0 and abs(pcs[4*t+8][1]) < abs(pcs[4*t+8][0]) and (pcs[4*t+8][0]**2. + pcs[4*t+8][1]**2.) >= 1.00 ) :  ### Jet Retraction
		JetRegime_verif = 2 
	
	if (pcs[4*t+8][1] > 0 and abs(pcs[4*t+8][0]) < pcs[4*t+8][1] and (pcs[4*t+8][0]**2. + pcs[4*t+8][1]**2.) >= 1.00 ) :    ### Poleward Shift
		JetRegime_verif = 3
	
	if (pcs[4*t+8][1] < 0 and abs(pcs[4*t+8][0]) < abs(pcs[4*t+8][1]) and (pcs[4*t+8][0]**2. + pcs[4*t+8][1]**2.) >= 1.00 ) :   ### Equatorward Shift
		JetRegime_verif = 4 
	
	if ((pcs[4*t+8][0]**2. + pcs[4*t+8][1]**2.) < 1.00) :    ### Origin
		JetRegime_verif = 5
	
	gefsfcastregime[gefsJetRegime_fcast] = 1
	gfsfcastregime[gfsJetRegime_fcast] = 1 
	verifregime[JetRegime_verif] = 1 
	
	if t == 0 :
		D1fcastreg_gfs = gfsfcastregime 
		D1fcastreg_gefs = gefsfcastregime
		D1verifreg = verifregime 
	if t == 1 :
		D2fcastreg_gfs = gfsfcastregime 
		D2fcastreg_gefs = gefsfcastregime
		D2verifreg = verifregime 
	if t == 2 :
		D3fcastreg_gfs = gfsfcastregime 
		D3fcastreg_gefs = gefsfcastregime
		D3verifreg = verifregime 
	if t == 3 :
		D4fcastreg_gfs = gfsfcastregime 
		D4fcastreg_gefs = gefsfcastregime
		D4verifreg = verifregime 
	if t == 4 :
		D5fcastreg_gfs = gfsfcastregime 
		D5fcastreg_gefs = gefsfcastregime
		D5verifreg = verifregime 
	if t == 5 :
		D6fcastreg_gfs = gfsfcastregime 
		D6fcastreg_gefs = gefsfcastregime
		D6verifreg = verifregime 
	if t == 6 :
		D7fcastreg_gfs = gfsfcastregime 
		D7fcastreg_gefs = gefsfcastregime
		D7verifreg = verifregime 
	if t == 7 :
		D8fcastreg_gfs = gfsfcastregime 
		D8fcastreg_gefs = gefsfcastregime
		D8verifreg = verifregime 
	if t == 8 :
		D9fcastreg_gfs = gfsfcastregime 
		D9fcastreg_gefs = gefsfcastregime
		D9verifreg = verifregime 

	####################### Save variables to the relevant verification files ############################

avedist = [JetMode,aveD1dist,aveD2dist,aveD3dist,aveD4dist,aveD5dist,aveD6dist,aveD7dist,aveD8dist,aveD9dist] 
gfsdist = [JetMode,D1gfsdist,D2gfsdist,D3gfsdist,D4gfsdist,D5gfsdist,D6gfsdist,D7gfsdist,D8gfsdist,D9gfsdist] 
gefsdist = [JetMode,D1gefsdist,D2gefsdist,D3gefsdist,D4gefsdist,D5gefsdist,D6gefsdist,D7gefsdist,D8gefsdist,D9gefsdist] 

D1probct = fcastct[0] 
D2probct = fcastct[1] 
D3probct = fcastct[2] 
D4probct = fcastct[3] 
D5probct = fcastct[4] 
D6probct = fcastct[5] 
D7probct = fcastct[6] 
D8probct = fcastct[7] 
D9probct = fcastct[8] 
D1probhit = fcasthit[0] 
D2probhit = fcasthit[1] 
D3probhit = fcasthit[2] 
D4probhit = fcasthit[3] 
D5probhit = fcasthit[4] 
D6probhit = fcasthit[5] 
D7probhit = fcasthit[6] 
D8probhit = fcasthit[7] 
D9probhit = fcasthit[8] 

#### full verification stats ###
avedistfile = 'avedist_gefs.out' 
gfsdistfile = 'dist_gfs.out'
gefsdistfile = 'dist_gefsmean.out'
podfile = 'pod.out'

### For Reliability Diagram ###
D1probctfile = 'D1probct.out'
D2probctfile = 'D2probct.out'
D3probctfile = 'D3probct.out'
D4probctfile = 'D4probct.out'
D5probctfile = 'D5probct.out'
D6probctfile = 'D6probct.out'
D7probctfile = 'D7probct.out'
D8probctfile = 'D8probct.out'
D9probctfile = 'D9probct.out'

D1hitfile = 'D1probhit.out'
D2hitfile = 'D2probhit.out'
D3hitfile = 'D3probhit.out'
D4hitfile = 'D4probhit.out'
D5hitfile = 'D5probhit.out'
D6hitfile = 'D6probhit.out'
D7hitfile = 'D7probhit.out'
D8hitfile = 'D8probhit.out'
D9hitfile = 'D9probhit.out'

### For forecasted/verified verified regime stats ###
D1gefsfcfile = 'D1fcastreg_gefs.out'
D2gefsfcfile = 'D2fcastreg_gefs.out'
D3gefsfcfile = 'D3fcastreg_gefs.out'
D4gefsfcfile = 'D4fcastreg_gefs.out'
D5gefsfcfile = 'D5fcastreg_gefs.out'
D6gefsfcfile = 'D6fcastreg_gefs.out'
D7gefsfcfile = 'D7fcastreg_gefs.out'
D8gefsfcfile = 'D8fcastreg_gefs.out'
D9gefsfcfile = 'D9fcastreg_gefs.out'

D1gfsfcfile = 'D1fcastreg_gfs.out'
D2gfsfcfile = 'D2fcastreg_gfs.out'
D3gfsfcfile = 'D3fcastreg_gfs.out'
D4gfsfcfile = 'D4fcastreg_gfs.out'
D5gfsfcfile = 'D5fcastreg_gfs.out'
D6gfsfcfile = 'D6fcastreg_gfs.out'
D7gfsfcfile = 'D7fcastreg_gfs.out'
D8gfsfcfile = 'D8fcastreg_gfs.out'
D9gfsfcfile = 'D9fcastreg_gfs.out'

D1verffile = 'D1verifreg.out'
D2verffile = 'D2verifreg.out'
D3verffile = 'D3verifreg.out'
D4verffile = 'D4verifreg.out'
D5verffile = 'D5verifreg.out'
D6verffile = 'D6verifreg.out'
D7verffile = 'D7verifreg.out'
D8verffile = 'D8verifreg.out'
D9verffile = 'D9verifreg.out'

##### Save and Write outfiles ######

#### Full verification stats #####
file = open(avedistfile,'a')
file.write(str(avedist[0])+','+str(avedist[1])+','+str(avedist[2])+','+str(avedist[3])+','+str(avedist[4])+','+str(avedist[5])+','+str(avedist[6])+','+str(avedist[7])+','+str(avedist[8])+','+str(avedist[9])+'\n')
file.close()

file = open(gfsdistfile,'a')
file.write(str(gfsdist[0])+','+str(gfsdist[1])+','+str(gfsdist[2])+','+str(gfsdist[3])+','+str(gfsdist[4])+','+str(gfsdist[5])+','+str(gfsdist[6])+','+str(gfsdist[7])+','+str(gfsdist[8])+','+str(gfsdist[9])+'\n')
file.close()

file = open(gefsdistfile,'a')
file.write(str(gefsdist[0])+','+str(gefsdist[1])+','+str(gefsdist[2])+','+str(gefsdist[3])+','+str(gefsdist[4])+','+str(gefsdist[5])+','+str(gefsdist[6])+','+str(gefsdist[7])+','+str(gefsdist[8])+','+str(gefsdist[9])+'\n')
file.close()

file = open(podfile,'a')
file.write(str(hitmiss[0])+','+str(hitmiss[1])+','+str(hitmiss[2])+','+str(hitmiss[3])+','+str(hitmiss[4])+','+str(hitmiss[5])+','+str(hitmiss[6])+','+str(hitmiss[7])+','+str(hitmiss[8])+','+str(hitmiss[9])+'\n')
file.close()

##### Reliability Diagram Stats #######

file = open(D1probctfile,'a')
file.write(str(D1probct[0])+','+str(D1probct[1])+','+str(D1probct[2])+','+str(D1probct[3])+','+str(D1probct[4])+','+str(D1probct[5])+','+str(D1probct[6])+','+str(D1probct[7])+','+str(D1probct[8])+','+str(D1probct[9])+','+str(D1probct[10])+','+str(D1probct[11])+','+str(D1probct[12])+'\n')
file.close()

file = open(D2probctfile,'a')
file.write(str(D2probct[0])+','+str(D2probct[1])+','+str(D2probct[2])+','+str(D2probct[3])+','+str(D2probct[4])+','+str(D2probct[5])+','+str(D2probct[6])+','+str(D2probct[7])+','+str(D2probct[8])+','+str(D2probct[9])+','+str(D2probct[10])+','+str(D2probct[11])+','+str(D2probct[12])+'\n')
file.close()

file = open(D3probctfile,'a')
file.write(str(D3probct[0])+','+str(D3probct[1])+','+str(D3probct[2])+','+str(D3probct[3])+','+str(D3probct[4])+','+str(D3probct[5])+','+str(D3probct[6])+','+str(D3probct[7])+','+str(D3probct[8])+','+str(D3probct[9])+','+str(D3probct[10])+','+str(D3probct[11])+','+str(D3probct[12])+'\n')
file.close()

file = open(D4probctfile,'a')
file.write(str(D4probct[0])+','+str(D4probct[1])+','+str(D4probct[2])+','+str(D4probct[3])+','+str(D4probct[4])+','+str(D4probct[5])+','+str(D4probct[6])+','+str(D4probct[7])+','+str(D4probct[8])+','+str(D4probct[9])+','+str(D4probct[10])+','+str(D4probct[11])+','+str(D4probct[12])+'\n')
file.close()

file = open(D5probctfile,'a')
file.write(str(D5probct[0])+','+str(D5probct[1])+','+str(D5probct[2])+','+str(D5probct[3])+','+str(D5probct[4])+','+str(D5probct[5])+','+str(D5probct[6])+','+str(D5probct[7])+','+str(D5probct[8])+','+str(D5probct[9])+','+str(D5probct[10])+','+str(D5probct[11])+','+str(D5probct[12])+'\n')
file.close()

file = open(D6probctfile,'a')
file.write(str(D6probct[0])+','+str(D6probct[1])+','+str(D6probct[2])+','+str(D6probct[3])+','+str(D6probct[4])+','+str(D6probct[5])+','+str(D6probct[6])+','+str(D6probct[7])+','+str(D6probct[8])+','+str(D6probct[9])+','+str(D6probct[10])+','+str(D6probct[11])+','+str(D6probct[12])+'\n')
file.close()

file = open(D7probctfile,'a')
file.write(str(D7probct[0])+','+str(D7probct[1])+','+str(D7probct[2])+','+str(D7probct[3])+','+str(D7probct[4])+','+str(D7probct[5])+','+str(D7probct[6])+','+str(D7probct[7])+','+str(D7probct[8])+','+str(D7probct[9])+','+str(D7probct[10])+','+str(D7probct[11])+','+str(D7probct[12])+'\n')
file.close()

file = open(D8probctfile,'a')
file.write(str(D8probct[0])+','+str(D8probct[1])+','+str(D8probct[2])+','+str(D8probct[3])+','+str(D8probct[4])+','+str(D8probct[5])+','+str(D8probct[6])+','+str(D8probct[7])+','+str(D8probct[8])+','+str(D8probct[9])+','+str(D8probct[10])+','+str(D8probct[11])+','+str(D8probct[12])+'\n')
file.close()

file = open(D9probctfile,'a')
file.write(str(D9probct[0])+','+str(D9probct[1])+','+str(D9probct[2])+','+str(D9probct[3])+','+str(D9probct[4])+','+str(D9probct[5])+','+str(D9probct[6])+','+str(D9probct[7])+','+str(D9probct[8])+','+str(D9probct[9])+','+str(D9probct[10])+','+str(D9probct[11])+','+str(D9probct[12])+'\n')
file.close()



file = open(D1hitfile,'a')
file.write(str(D1probhit[0])+','+str(D1probhit[1])+','+str(D1probhit[2])+','+str(D1probhit[3])+','+str(D1probhit[4])+','+str(D1probhit[5])+','+str(D1probhit[6])+','+str(D1probhit[7])+','+str(D1probhit[8])+','+str(D1probhit[9])+','+str(D1probhit[10])+','+str(D1probhit[11])+','+str(D1probhit[12])+'\n')
file.close()

file = open(D2hitfile,'a')
file.write(str(D2probhit[0])+','+str(D2probhit[1])+','+str(D2probhit[2])+','+str(D2probhit[3])+','+str(D2probhit[4])+','+str(D2probhit[5])+','+str(D2probhit[6])+','+str(D2probhit[7])+','+str(D2probhit[8])+','+str(D2probhit[9])+','+str(D2probhit[10])+','+str(D2probhit[11])+','+str(D2probhit[12])+'\n')
file.close()

file = open(D3hitfile,'a')
file.write(str(D3probhit[0])+','+str(D3probhit[1])+','+str(D3probhit[2])+','+str(D3probhit[3])+','+str(D3probhit[4])+','+str(D3probhit[5])+','+str(D3probhit[6])+','+str(D3probhit[7])+','+str(D3probhit[8])+','+str(D3probhit[9])+','+str(D3probhit[10])+','+str(D3probhit[11])+','+str(D3probhit[12])+'\n')
file.close()

file = open(D4hitfile,'a')
file.write(str(D4probhit[0])+','+str(D4probhit[1])+','+str(D4probhit[2])+','+str(D4probhit[3])+','+str(D4probhit[4])+','+str(D4probhit[5])+','+str(D4probhit[6])+','+str(D4probhit[7])+','+str(D4probhit[8])+','+str(D4probhit[9])+','+str(D4probhit[10])+','+str(D4probhit[11])+','+str(D4probhit[12])+'\n')
file.close()

file = open(D5hitfile,'a')
file.write(str(D5probhit[0])+','+str(D5probhit[1])+','+str(D5probhit[2])+','+str(D5probhit[3])+','+str(D5probhit[4])+','+str(D5probhit[5])+','+str(D5probhit[6])+','+str(D5probhit[7])+','+str(D5probhit[8])+','+str(D5probhit[9])+','+str(D5probhit[10])+','+str(D5probhit[11])+','+str(D5probhit[12])+'\n')
file.close()

file = open(D6hitfile,'a')
file.write(str(D6probhit[0])+','+str(D6probhit[1])+','+str(D6probhit[2])+','+str(D6probhit[3])+','+str(D6probhit[4])+','+str(D6probhit[5])+','+str(D6probhit[6])+','+str(D6probhit[7])+','+str(D6probhit[8])+','+str(D6probhit[9])+','+str(D6probhit[10])+','+str(D6probhit[11])+','+str(D6probhit[12])+'\n')
file.close()

file = open(D7hitfile,'a')
file.write(str(D7probhit[0])+','+str(D7probhit[1])+','+str(D7probhit[2])+','+str(D7probhit[3])+','+str(D7probhit[4])+','+str(D7probhit[5])+','+str(D7probhit[6])+','+str(D7probhit[7])+','+str(D7probhit[8])+','+str(D7probhit[9])+','+str(D7probhit[10])+','+str(D7probhit[11])+','+str(D7probhit[12])+'\n')
file.close()

file = open(D8hitfile,'a')
file.write(str(D8probhit[0])+','+str(D8probhit[1])+','+str(D8probhit[2])+','+str(D8probhit[3])+','+str(D8probhit[4])+','+str(D8probhit[5])+','+str(D8probhit[6])+','+str(D8probhit[7])+','+str(D8probhit[8])+','+str(D8probhit[9])+','+str(D8probhit[10])+','+str(D8probhit[11])+','+str(D8probhit[12])+'\n')
file.close()

file = open(D9hitfile,'a')
file.write(str(D9probhit[0])+','+str(D9probhit[1])+','+str(D9probhit[2])+','+str(D9probhit[3])+','+str(D9probhit[4])+','+str(D9probhit[5])+','+str(D9probhit[6])+','+str(D9probhit[7])+','+str(D9probhit[8])+','+str(D9probhit[9])+','+str(D9probhit[10])+','+str(D9probhit[11])+','+str(D9probhit[12])+'\n')
file.close()

##### Forecasted/Verified Regime Stats #####

file = open(D1gfsfcfile,'a')
file.write(str(D1fcastreg_gfs[0])+','+str(D1fcastreg_gfs[1])+','+str(D1fcastreg_gfs[2])+','+str(D1fcastreg_gfs[3])+','+str(D1fcastreg_gfs[4])+','+str(D1fcastreg_gfs[5])+'\n')
file.close()

file = open(D2gfsfcfile,'a')
file.write(str(D2fcastreg_gfs[0])+','+str(D2fcastreg_gfs[1])+','+str(D2fcastreg_gfs[2])+','+str(D2fcastreg_gfs[3])+','+str(D2fcastreg_gfs[4])+','+str(D2fcastreg_gfs[5])+'\n')
file.close()

file = open(D3gfsfcfile,'a')
file.write(str(D3fcastreg_gfs[0])+','+str(D3fcastreg_gfs[1])+','+str(D3fcastreg_gfs[2])+','+str(D3fcastreg_gfs[3])+','+str(D3fcastreg_gfs[4])+','+str(D3fcastreg_gfs[5])+'\n')
file.close()

file = open(D4gfsfcfile,'a')
file.write(str(D4fcastreg_gfs[0])+','+str(D4fcastreg_gfs[1])+','+str(D4fcastreg_gfs[2])+','+str(D4fcastreg_gfs[3])+','+str(D4fcastreg_gfs[4])+','+str(D4fcastreg_gfs[5])+'\n')
file.close()

file = open(D5gfsfcfile,'a')
file.write(str(D5fcastreg_gfs[0])+','+str(D5fcastreg_gfs[1])+','+str(D5fcastreg_gfs[2])+','+str(D5fcastreg_gfs[3])+','+str(D5fcastreg_gfs[4])+','+str(D5fcastreg_gfs[5])+'\n')
file.close()

file = open(D6gfsfcfile,'a')
file.write(str(D6fcastreg_gfs[0])+','+str(D6fcastreg_gfs[1])+','+str(D6fcastreg_gfs[2])+','+str(D6fcastreg_gfs[3])+','+str(D6fcastreg_gfs[4])+','+str(D6fcastreg_gfs[5])+'\n')
file.close()

file = open(D7gfsfcfile,'a')
file.write(str(D7fcastreg_gfs[0])+','+str(D7fcastreg_gfs[1])+','+str(D7fcastreg_gfs[2])+','+str(D7fcastreg_gfs[3])+','+str(D7fcastreg_gfs[4])+','+str(D7fcastreg_gfs[5])+'\n')
file.close()

file = open(D8gfsfcfile,'a')
file.write(str(D8fcastreg_gfs[0])+','+str(D8fcastreg_gfs[1])+','+str(D8fcastreg_gfs[2])+','+str(D8fcastreg_gfs[3])+','+str(D8fcastreg_gfs[4])+','+str(D8fcastreg_gfs[5])+'\n')
file.close()

file = open(D9gfsfcfile,'a')
file.write(str(D9fcastreg_gfs[0])+','+str(D9fcastreg_gfs[1])+','+str(D9fcastreg_gfs[2])+','+str(D9fcastreg_gfs[3])+','+str(D9fcastreg_gfs[4])+','+str(D9fcastreg_gfs[5])+'\n')
file.close()



file = open(D1gefsfcfile,'a')
file.write(str(D1fcastreg_gefs[0])+','+str(D1fcastreg_gefs[1])+','+str(D1fcastreg_gefs[2])+','+str(D1fcastreg_gefs[3])+','+str(D1fcastreg_gefs[4])+','+str(D1fcastreg_gefs[5])+'\n')
file.close()

file = open(D2gefsfcfile,'a')
file.write(str(D2fcastreg_gefs[0])+','+str(D2fcastreg_gefs[1])+','+str(D2fcastreg_gefs[2])+','+str(D2fcastreg_gefs[3])+','+str(D2fcastreg_gefs[4])+','+str(D2fcastreg_gefs[5])+'\n')
file.close()

file = open(D3gefsfcfile,'a')
file.write(str(D3fcastreg_gefs[0])+','+str(D3fcastreg_gefs[1])+','+str(D3fcastreg_gefs[2])+','+str(D3fcastreg_gefs[3])+','+str(D3fcastreg_gefs[4])+','+str(D3fcastreg_gefs[5])+'\n')
file.close()

file = open(D4gefsfcfile,'a')
file.write(str(D4fcastreg_gefs[0])+','+str(D4fcastreg_gefs[1])+','+str(D4fcastreg_gefs[2])+','+str(D4fcastreg_gefs[3])+','+str(D4fcastreg_gefs[4])+','+str(D4fcastreg_gefs[5])+'\n')
file.close()

file = open(D5gefsfcfile,'a')
file.write(str(D5fcastreg_gefs[0])+','+str(D5fcastreg_gefs[1])+','+str(D5fcastreg_gefs[2])+','+str(D5fcastreg_gefs[3])+','+str(D5fcastreg_gefs[4])+','+str(D5fcastreg_gefs[5])+'\n')
file.close()

file = open(D6gefsfcfile,'a')
file.write(str(D6fcastreg_gefs[0])+','+str(D6fcastreg_gefs[1])+','+str(D6fcastreg_gefs[2])+','+str(D6fcastreg_gefs[3])+','+str(D6fcastreg_gefs[4])+','+str(D6fcastreg_gefs[5])+'\n')
file.close()

file = open(D7gefsfcfile,'a')
file.write(str(D7fcastreg_gefs[0])+','+str(D7fcastreg_gefs[1])+','+str(D7fcastreg_gefs[2])+','+str(D7fcastreg_gefs[3])+','+str(D7fcastreg_gefs[4])+','+str(D7fcastreg_gefs[5])+'\n')
file.close()

file = open(D8gefsfcfile,'a')
file.write(str(D8fcastreg_gefs[0])+','+str(D8fcastreg_gefs[1])+','+str(D8fcastreg_gefs[2])+','+str(D8fcastreg_gefs[3])+','+str(D8fcastreg_gefs[4])+','+str(D8fcastreg_gefs[5])+'\n')
file.close()

file = open(D9gefsfcfile,'a')
file.write(str(D9fcastreg_gefs[0])+','+str(D9fcastreg_gefs[1])+','+str(D9fcastreg_gefs[2])+','+str(D9fcastreg_gefs[3])+','+str(D9fcastreg_gefs[4])+','+str(D9fcastreg_gefs[5])+'\n')
file.close()



file = open(D1verffile,'a')
file.write(str(D1verifreg[0])+','+str(D1verifreg[1])+','+str(D1verifreg[2])+','+str(D1verifreg[3])+','+str(D1verifreg[4])+','+str(D1verifreg[5])+'\n')
file.close()

file = open(D2verffile,'a')
file.write(str(D2verifreg[0])+','+str(D2verifreg[1])+','+str(D2verifreg[2])+','+str(D2verifreg[3])+','+str(D2verifreg[4])+','+str(D2verifreg[5])+'\n')
file.close()

file = open(D3verffile,'a')
file.write(str(D3verifreg[0])+','+str(D3verifreg[1])+','+str(D3verifreg[2])+','+str(D3verifreg[3])+','+str(D3verifreg[4])+','+str(D3verifreg[5])+'\n')
file.close()

file = open(D4verffile,'a')
file.write(str(D4verifreg[0])+','+str(D4verifreg[1])+','+str(D4verifreg[2])+','+str(D4verifreg[3])+','+str(D4verifreg[4])+','+str(D4verifreg[5])+'\n')
file.close()

file = open(D5verffile,'a')
file.write(str(D5verifreg[0])+','+str(D5verifreg[1])+','+str(D5verifreg[2])+','+str(D5verifreg[3])+','+str(D5verifreg[4])+','+str(D5verifreg[5])+'\n')
file.close()

file = open(D6verffile,'a')
file.write(str(D6verifreg[0])+','+str(D6verifreg[1])+','+str(D6verifreg[2])+','+str(D6verifreg[3])+','+str(D6verifreg[4])+','+str(D6verifreg[5])+'\n')
file.close()

file = open(D7verffile,'a')
file.write(str(D7verifreg[0])+','+str(D7verifreg[1])+','+str(D7verifreg[2])+','+str(D7verifreg[3])+','+str(D7verifreg[4])+','+str(D7verifreg[5])+'\n')
file.close()

file = open(D8verffile,'a')
file.write(str(D8verifreg[0])+','+str(D8verifreg[1])+','+str(D8verifreg[2])+','+str(D8verifreg[3])+','+str(D8verifreg[4])+','+str(D8verifreg[5])+'\n')
file.close()

file = open(D9verffile,'a')
file.write(str(D9verifreg[0])+','+str(D9verifreg[1])+','+str(D9verifreg[2])+','+str(D9verifreg[3])+','+str(D9verifreg[4])+','+str(D9verifreg[5])+'\n')
file.close()