#!/usr/bin/env python
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from scipy import interpolate
from math import erf

# v=0, J=2, K=1, mJ sampled randomly
# Ei, sticking of all mk states, mK=-1, 0, 1, seed 11, 10000 trajectories
HSE03_13X_random = np.array(
[[0.05, 0.02160, 0.02370, 0.01680, 0.02430],
[0.075, 0.11428, 0.12418, 0.08308, 0.13557],
[0.1, 0.26647, 0.30539, 0.16168, 0.33234],#1000 trajectories
[0.125, 0.45907, 0.52231, 0.32208, 0.53281],#seed 8, 8000 trajectories
[0.15, 0.62275, 0.68862, 0.45210, 0.72754],#seed 10, 1000 trajectories
[0.175, 0.76462, 0.83058, 0.61169, 0.85157],#seed 10, 2000 trajectories
[0.2, 0.88922, 0.92216, 0.81138, 0.93413],#1000 trajectories
[0.225, 0.95908, 0.96707, 0.92814, 0.98204],#seed 8, 1000 trajectories
[0.25, 0.98703, 0.99401, 0.97605, 0.99102]]#1000 trajectories
)
HSE03_13X_random[:,0] *= 96.4853075

HSE03_13X_random_cartwheel = (3.*HSE03_13X_random[:,1]-HSE03_13X_random[:,4])/2.
HSE03_13X_random_perpendicular = 2.*HSE03_13X_random_cartwheel - HSE03_13X_random[:,4]

# Phi is 0 degrees
HSE03_13X_0180 = np.array(
[[0., 0.79733, 0.87800, 0.66300, 0.85100],
[5., 0.78058, 0.86331, 0.64269, 0.83573],
[10., 0.76898, 0.82734, 0.62230, 0.85731],
[15., 0.72382, 0.80695, 0.56355, 0.80096],
[20., 0.69664, 0.79017, 0.52038, 0.77938],
[25., 0.59133, 0.67500, 0.41600, 0.68300],
[30., 0.48715, 0.56424, 0.33084, 0.56638]]
#[35., ],
#[40., ]]
)
HSE03_13X_0180_cartwheel = (3.*HSE03_13X_0180[:,1]-HSE03_13X_0180[:,4])/2.
HSE03_13X_0180_perpendicular = 2.*HSE03_13X_0180_cartwheel - HSE03_13X_0180[:,4]

# v=0, J=2, K=1, mJ sampled randomly
#Ei=0.115 eV, angle / degrees, sticking of all mk states, mK=-1, 0, 1
HSE03_13X_0115_phirandom = np.array(
[[0., 0.36000, 0.43600, 0.24200, 0.40200],
[10., 0.33200, 0.39500, 0.21500, 0.38600],
[15., 0.30600, 0.36100, 0.19700, 0.36000],
[20., 0.24400, 0.27900, 0.16200, 0.29100],
[25., 0.19867, 0.23000, 0.13500, 0.23100],
[30., 0.14267, 0.17300, 0.09900, 0.15600]]
)

# Phi is 30 degrees
HSE03_13X_0115_phi30 = np.array(
[[0., 0.36000, 0.43600, 0.24200, 0.40200],
[10., 0.32467, 0.38300, 0.21500, 0.37600],
[15., 0.30650, 0.34100, 0.22300, 0.35550],
[20., 0.24900, 0.28200, 0.17900, 0.28600],
[25., 0.19700, 0.22200, 0.15300, 0.21600],
[30., 0.14100, 0.16100, 0.11300, 0.14900],
[35., 0.08067, 0.08800, 0.05900, 0.09500],
[40., 0.03933, 0.04200, 0.02800, 0.04800]]
)

# Phi is 0 degrees
HSE03_13X_0115 = np.array(
[[0., 0.36000, 0.43600, 0.24200, 0.40200],
[5., 0.35967, 0.40800, 0.24800, 0.42300],
[10., 0.33833, 0.39600, 0.22900, 0.39000],
[15., 0.31500, 0.37500, 0.20700, 0.36300],
[20., 0.26600, 0.32500, 0.16900, 0.30400],
[25., 0.20300, 0.25200, 0.13000, 0.22700],
[30., 0.13033, 0.16800, 0.08000, 0.14300],
[35., 0.08733, 0.09600, 0.06300, 0.10300],
[40., 0.03933, 0.03500, 0.02800, 0.05500]]
)
HSE03_13X_0115_cartwheel = (3.*HSE03_13X_0115[:,1]-HSE03_13X_0115[:,4])/2.
HSE03_13X_0115_perpendicular = 2.*HSE03_13X_0115_cartwheel - HSE03_13X_0115[:,4]

HSE03_13X_0085_phi30 = np.array(
[[0., 0.16200, 0.19000, 0.11500, 0.18100],
[10, 0.15567, 0.17100, 0.10800, 0.18800],
[15., 0.13367, 0.14900, 0.10100, 0.15100],
[20., 0.11000, 0.11700, 0.08900, 0.12400]]
)

# Phi is 0 degrees
HSE03_13X_0085 = np.array(
[[0., 0.16200, 0.19000, 0.11500, 0.18100],
[5., 0.16133, 0.16600, 0.12400, 0.19400],
[10, 0.15467, 0.18100, 0.10300, 0.18000],
[15., 0.13533, 0.16400, 0.09300, 0.14900],
[20., 0.10333, 0.13100, 0.07600, 0.10300],
[25., 0.07533, 0.08100, 0.05200, 0.09300],
[30., 0.04533, 0.06100, 0.02700, 0.04800],
[35., 0.01833, 0.01600, 0.01500, 0.02400]]
)
HSE03_13X_0085_cartwheel = (3.*HSE03_13X_0085[:,1]-HSE03_13X_0085[:,4])/2.
HSE03_13X_0085_perpendicular = 2.*HSE03_13X_0085_cartwheel - HSE03_13X_0085[:,4]

HSE03_13X_0090_phi30 = np.array(
[[0., 0.18100, 0.21200, 0.13000, 0.20100],
[5., 0.18267, 0.20800, 0.13000, 0.21000],
[10., 0.17767, 0.20100, 0.13100, 0.20100],
[15., 0.15200, 0.17400, 0.11700, 0.16500],
[20., 0.12200, 0.14200, 0.09500, 0.12900],
[25., 0.08833, 0.09500, 0.07300, 0.09700],
[30., 0.05533, 0.05900, 0.04200, 0.06500],
[35., 0.02833, 0.03600, 0.02000, 0.02900],
[40., 0.03933, 0.04200, 0.02800, 0.04800]]
)
#HSE03_13X_0090_cartwheel = (3.*HSE03_13X_0090[:,1]-HSE03_13X_0090[:,4])/2.
#HSE03_13X_0090_perpendicular = 2.*HSE03_13X_0090_cartwheel - HSE03_13X_0090[:,4]

#angle / degrees, S0 helicopter, S0 cartwheel
ECW_023 = np.array(
[[0., 0.42329, 0.32440],
[10., 0.39556, 0.32625],
[20., 0.34750, 0.29945],
[30., 0.29390, 0.25323],
[40., 0.23013, 0.19039]]
)
ECW_023_NES_heli = np.array(
[[1.03156, 0.41999],
[1.61142, 0.41949],
[2.19128, 0.41907],
[3.35091, 0.41712],
[3.93075, 0.41646],
[5.09037, 0.41435],
[6.24992, 0.41135],
[6.82966, 0.40955],
[7.98920, 0.40642],
[8.56892, 0.40428],
[11.25644, 0.39118],
[11.59897, 0.38950],
[12.46842, 0.38481],
[12.94272, 0.38278],
[14.10200, 0.37660],
[14.54990, 0.37418],
[15.47204, 0.36917],
[15.89357, 0.36659],
[17.15808, 0.35797],
[18.58048, 0.34616],
[18.89659, 0.34380],
[19.95024, 0.33526],
[20.29268, 0.33248],
[21.37262, 0.32319],
[21.68876, 0.32117],
[22.74231, 0.31147],
[23.08477, 0.30903],
[24.16463, 0.29875],
[24.48072, 0.29605],
[25.53423, 0.28592],
[25.85037, 0.28390],
[26.85115, 0.27355],
[27.16724, 0.27092],
[28.22071, 0.26030],
[28.56316, 0.25773],
[29.59029, 0.24725],
[29.93273, 0.24454],
[30.90715, 0.23419],
[31.24956, 0.23115],
[32.27670, 0.22076],
[32.59282, 0.21859],
[33.59354, 0.20744],
[33.90962, 0.20477],
[34.91034, 0.19370],
[35.22638, 0.19054],
[36.22710, 0.17943],
[37.49114, 0.16512],
[38.80781, 0.14983],
[39.91379, 0.13663]]
)
ECW_023_NES_cart = np.array(
[[1.78382, 0.32581],
[2.30117, 0.32562],
[3.29057, 0.32571],
[3.87042, 0.32518],
[5.03012, 0.32402],
[5.45179, 0.32310],
[6.61124, 0.31893],
[6.95384, 0.31820],
[8.13972, 0.31487],
[8.50865, 0.31370],
[9.72081, 0.30947],
[10.06340, 0.30852],
[11.24913, 0.30351],
[11.59171, 0.30242],
[12.77742, 0.29717],
[13.12000, 0.29612],
[14.30568, 0.29038],
[14.64823, 0.28898],
[15.83387, 0.28286],
[16.17644, 0.28164],
[17.30933, 0.27515],
[17.65184, 0.27325],
[18.83749, 0.26722],
[19.15367, 0.26569],
[20.31290, 0.25887],
[20.65542, 0.25709],
[21.78829, 0.25037],
[22.13082, 0.24870],
[23.31637, 0.24141],
[23.63253, 0.23969],
[24.73899, 0.23230],
[25.08149, 0.23026],
[26.21427, 0.22246],
[26.55679, 0.22061],
[27.63682, 0.21237],
[27.97930, 0.21014],
[29.05931, 0.20166],
[29.40177, 0.19925],
[30.48174, 0.19024],
[30.82420, 0.18773],
[31.90412, 0.17821],
[32.22021, 0.17558],
[33.22110, 0.16659],
[33.56356, 0.16406],
[34.64342, 0.15380],
[34.95949, 0.15087],
[36.01303, 0.14113],
[36.32914, 0.13873],
[37.38263, 0.12830],
[37.69876, 0.12616],
[38.75224, 0.11558],
[39.06834, 0.11318],
[39.93748, 0.10470]]
)

ECW_015 = np.array(
[[0., 0.17652, 0.14880],
[10., 0.16266, 0.13494],
[20., 0.12662, 0.10628],
[30., 0.09982, 0.08133],
[40., 0.05915, 0.05176]]
)

#Ei=0.18 eV, angle / degrees, S0 helicopter, S0 cartwheel
exp_018 = np.array(
[[1.54186, 0.38593, 0.29347],
[9.66381, 0.37119, 0.27538],
[16.58333, 0.33233, 0.24523],
[22.65457, 0.29012, 0.20838],
[31.51559, 0.24054, 0.16750],
[39.22170, 0.15745, 0.10854],
[47.27960, 0.08576, 0.05628]]
)
exp_018_NES_heli = np.array(
[[0.8674101610904597, 0.4041811235026848],
[1.7161240512701372, 0.4037830409879757],
[2.3342057811567773, 0.40359371149295753],
[2.8706540750206546, 0.40300019488957106],
[3.4479190868959133, 0.4023360618301099],
[4.025184098771172, 0.4015115111175553],
[4.602449110646431, 0.4005907098131446],
[5.179714122521689, 0.3997340755699713],
[5.756979134396948, 0.3984924389593738],
[6.334244146272207, 0.3972187188181576],
[6.9115091581474655, 0.3959449986769414],
[7.488774170022724, 0.39447877735201303],
[8.066039181897983, 0.39288422190461],
[8.643304193773242, 0.3911292488041134],
[9.2205692056485, 0.38927802511176085],
[9.71911626135895, 0.38787465888213385],
[10.637492416615043, 0.3844891109058073],
[11.214757428490302, 0.3821887177847931],
[11.79202244036556, 0.37979207407192284],
[12.369287452240819, 0.3773954303590526],
[12.946552464116078, 0.3748704525237076],
[13.523817475991336, 0.37221714056588784],
[14.101082487866595, 0.36937132742435597],
[14.678347499741854, 0.36655759781344277],
[15.255612511617112, 0.36361553408005487],
[15.832877523492371, 0.36054513622419215],
[16.41014253536763, 0.3572822371846174],
[16.98740754724289, 0.35408350520627996],
[17.564672559118147, 0.35056393792175566],
[18.141937570993406, 0.3471727047597061],
[18.719202582868665, 0.34355688688332575],
[19.296467594743923, 0.33990898547632664],
[19.873732606619182, 0.33606858288561553],
[20.45099761849444, 0.3321319297030483],
[21.0282626303697, 0.32822736005109976],
[21.605527642244958, 0.3242907068685325],
[22.182792654120217, 0.32009738544101585],
[22.760057665995475, 0.3158398969522617],
[23.337322677870734, 0.3115824084635076],
[23.914587689745993, 0.3071965858522788],
[24.49185270162125, 0.3027145126491939],
[25.06911771349651, 0.2981361888542529],
[25.64638272537177, 0.29352578152869324],
[26.223647737247028, 0.28885120714189627],
[26.80091274912229, 0.2840482986326245],
[27.378177760997545, 0.279213306592734],
[27.955442772872807, 0.27424998043036874],
[28.532707784748062, 0.2692224872067662],
[29.109972796623325, 0.26419499398316365],
[29.68723780849858, 0.25919958429017975],
[30.264502820373842, 0.25397958988286506],
[32.888434692534105, 0.22954358821867632],
[33.46569970440937, 0.2241631761582682],
[34.04296471628462, 0.218686513506004],
[34.620229728159885, 0.21314568379250243],
[35.19749474003514, 0.20750860348714484],
[35.7747597519104, 0.20177527258993116],
[36.35202476378566, 0.1960740252233361],
[36.92928977566092, 0.19043694491797847],
[37.506554787536174, 0.18463944695952744],
[38.08381979941144, 0.17880986547045763],
[38.63484549256509, 0.17357018362279303],
[40.60279439668529, 0.15312085682422538],
[41.180059408560545, 0.14725919180453695],
[41.49493123321978, 0.14489025840249847]]
)
exp_018_NES_cart = np.array(
[[0.82399, 0.28114],
[1.40125, 0.28080],
[1.97852, 0.28062],
[2.55578, 0.28008],
[3.13305, 0.27948],
[3.71031, 0.27882],
[4.28758, 0.27812],
[4.86484, 0.27704],
[5.44211, 0.27599],
[6.01937, 0.27488],
[6.59664, 0.27357],
[7.17390, 0.27224],
[7.75117, 0.27067],
[8.32843, 0.26905],
[8.69578, 0.26798],
[9.69288, 0.26438],
[10.27014, 0.26285],
[10.84741, 0.26075],
[11.42467, 0.25864],
[12.00194, 0.25637],
[12.57920, 0.25404],
[13.15647, 0.25161],
[13.73373, 0.24899],
[14.31100, 0.24643],
[14.88826, 0.24375],
[15.43929, 0.24109],
[16.14775, 0.23757],
[16.98741, 0.23316],
[17.56467, 0.23038],
[18.14194, 0.22731],
[18.71920, 0.22424],
[19.29647, 0.22104],
[19.87373, 0.21790],
[20.45100, 0.21458],
[21.02826, 0.21125],
[21.57929, 0.20812],
[22.31399, 0.20380],
[23.54724, 0.19631],
[24.12450, 0.19305],
[24.70177, 0.18959],
[25.27903, 0.18607],
[25.85630, 0.18271],
[26.43356, 0.17929],
[27.01083, 0.17580],
[27.58809, 0.17228],
[28.16536, 0.16873],
[28.74262, 0.16544],
[29.31989, 0.16198],
[29.89715, 0.15849],
[30.47442, 0.15510],
[31.05168, 0.15174],
[31.62895, 0.14842],
[32.20621, 0.14512],
[32.78348, 0.14183],
[33.36074, 0.13869],
[33.93801, 0.13546],
[34.51527, 0.13242],
[35.09254, 0.12935],
[35.66980, 0.12622],
[36.24707, 0.12331],
[36.82433, 0.12046],
[37.40160, 0.11762],
[37.97886, 0.11477],
[38.55613, 0.11220],
[39.13339, 0.10946],
[39.52698, 0.10762],
[40.39288, 0.10390],
[40.97014, 0.10154],
[41.38997, 0.09986]]
)

#Ei=0.1 eV, angle / degrees, S0 helicopter, S0 cartwheel
exp_010 = np.array(
[[1.473886548060996, 0.18358458961473983, 0.11323283082077007],
[9.594499583327071, 0.1554438860971517, 0.09313232830820717],
[16.522041914457617, 0.1407035175879392, 0.08643216080401939],
[22.600456586688843, 0.1172529313232823, 0.07303182579564438],
[31.468826351765465, 0.08375209380234472, 0.052931323283081366],
[39.19280683911416, 0.05896147403685037, 0.04154103852596247]]
)

exp_osterlund_0183 = np.array(
[[0.10017, 0.41850],
[0.10017, 0.44700],
[10.16694, 0.41792],
[20.18364, 0.37783],
[25.19199, 0.39779],
[30.15025, 0.33725],
[35.10851, 0.16671],
[40.06678, 0.05567],
[45.12521, -0.00388]]
)

exp_osterlund_0183_NES = np.array(
[[0.10017, 0.42900],
[5.15860, 0.42046],
[10.16694, 0.39992],
[15.17529, 0.39037],
[20.13356, 0.35083],
[25.14190, 0.31029],
[30.15025, 0.26025],
[35.10851, 0.21021],
[40.11686, 0.17967],
[45.12521, 0.11862],
[50.08347, 0.07808],
[55.09182, 0.01804],
[57.99666, -0.00098]]
)

exp_osterlund_0112 = np.array(
[[0.15025, 0.20150],
[0.15025, 0.17100],
[10.11686, 0.17092],
[10.16694, 0.15892],
[20.13356, 0.19083],
[25.14190, 0.24079],
[30.15025, 0.19025],
[35.10851, 0.06821],
[40.01669, -0.00083]]
)

exp_osterlund_0112_NES = np.array(
[[0.20033, 0.19200],
[5.15860, 0.18596],
[10.16694, 0.18092],
[15.12521, 0.17037],
[20.13356, 0.15033],
[25.14190, 0.12479],
[30.10017, 0.09875],
[35.10851, 0.07821],
[40.06678, 0.03817]]
)

Nrows = 2
Ncols = 1
fig, ax = plt.subplots(nrows=Nrows, ncols=Ncols, figsize=(3.69,4.))
#fig = plt.figure(figsize=(3.69,3))

ax2 = plt.subplot(Nrows,Ncols,1)

Order = 3
Knots = 2
xaxis = np.linspace( min(exp_018[:,0]), max(exp_018[:,0]), 50 )
fexp_018_heli = interpolate.UnivariateSpline( exp_018[:,0], exp_018[:,1], s=Knots, k=Order)
fexp_018_perp = interpolate.UnivariateSpline( exp_018[:,0], exp_018[:,2], s=Knots, k=Order)
#plt.plot( xaxis, fexp_018_heli(xaxis), c='r', linestyle='--' )
#plt.plot( xaxis, fexp_018_perp(xaxis), c='r', linestyle='--' )

xaxis = np.linspace( min(exp_010[:,0]), max(exp_010[:,0]), 50 )
fexp_010_heli = interpolate.UnivariateSpline( exp_010[:,0], exp_010[:,1], s=Knots, k=Order)
fexp_010_perp = interpolate.UnivariateSpline( exp_010[:,0], exp_010[:,2], s=Knots, k=Order)
#plt.plot( xaxis, fexp_010_heli(xaxis), c='b', linestyle='--' )
#plt.plot( xaxis, fexp_010_perp(xaxis), c='b', linestyle='--' )

Knots = 0
xaxis = np.linspace( 0., 40., 40 )
fHSE03_13X_heli = interpolate.UnivariateSpline( HSE03_13X_random[:,0], HSE03_13X_random[:,4], s=Knots, k=Order)
fHSE03_13X_cart = interpolate.UnivariateSpline( HSE03_13X_random[:,0], HSE03_13X_random_cartwheel[:], s=Knots, k=Order)
S0115_heli_NES = []
S0115_cart_NES = []
for alpha in xaxis:
	EN = 0.115 * 96.485307499 * np.cos( np.radians(alpha) )**2
	S0115_heli_NES.append( fHSE03_13X_heli(EN) )
	S0115_cart_NES.append( fHSE03_13X_cart(EN) )
plt.plot( xaxis, S0115_heli_NES, color='r', linestyle='-' )
plt.plot( xaxis, S0115_cart_NES, color='r', linestyle='-' )

S0085_heli_NES = []
S0085_cart_NES = []
for alpha in xaxis:
	EN = 0.085 * 96.485307499 * np.cos( np.radians(alpha) )**2
	S0085_heli_NES.append( fHSE03_13X_heli(EN) )
	S0085_cart_NES.append( fHSE03_13X_cart(EN) )
plt.plot( xaxis, S0085_heli_NES, color='b', linestyle='-' )
plt.plot( xaxis, S0085_cart_NES, color='b', linestyle='-' )

plt.plot( exp_018_NES_heli[:,0], exp_018_NES_heli[:,1], c='r', linestyle='--' )
plt.plot( exp_018_NES_cart[:,0], exp_018_NES_cart[:,1], c='r', linestyle='--' )

plt.scatter( exp_010[:,0], exp_010[:,1], edgecolor='b', facecolor='None', marker='v', label='')
plt.scatter( exp_010[:,0], exp_010[:,2], edgecolor='b', facecolor='None', marker='^', label='')
plt.scatter( exp_018[:,0], exp_018[:,1], edgecolor='r', facecolor='None', marker='v', label='')
plt.scatter( exp_018[:,0], exp_018[:,2], edgecolor='r', facecolor='None', marker='^', label='')

#plt.plot( ECW_015[:,0], ECW_015[:,1], color='b', markerfacecolor='k', marker='v', label='')
#plt.plot( ECW_015[:,0], ECW_015[:,2], color='b', markerfacecolor='k', marker='^', label='')
#plt.plot( ECW_023[:,0], ECW_023[:,1], color='r', markerfacecolor='k', marker='v', label='')
#plt.plot( ECW_023[:,0], ECW_023[:,2], color='r', markerfacecolor='k', marker='^', label='')

plt.scatter( HSE03_13X_0115[:,0], HSE03_13X_0115[:,4], c='r', marker='v', label='helicopter' )
plt.scatter( HSE03_13X_0115[:,0], HSE03_13X_0115_cartwheel[:], c='r', marker='^', label='cartwheel' )
plt.scatter( HSE03_13X_0085[:,0], HSE03_13X_0085[:,4], c='b', marker='v', label='' )
plt.scatter( HSE03_13X_0085[:,0], HSE03_13X_0085_cartwheel[:], c='b', marker='^', label='' )

plt.annotate('(a)', xy=(43, 0.22), size=12)

plt.legend(loc='upper right', numpoints=1, scatterpoints=1)
plt.xlim(0.,50.)
plt.ylim(0.,0.5)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel('Reaction probability')
#plt.xlabel('Incidence angle (deg.)')

ax2 = plt.subplot(Nrows,Ncols,2)

#plt.plot( xaxis, fexp_018_heli(xaxis), c='r', linestyle='--' )
#plt.plot( xaxis, fexp_018_perp(xaxis), c='r', linestyle='--' )

#plt.plot( xaxis, fexp_010_heli(xaxis), c='b', linestyle='--' )
#plt.plot( xaxis, fexp_010_perp(xaxis), c='b', linestyle='--' )

#plt.scatter( exp_010[:,0], exp_010[:,1], edgecolor='b', facecolor='None', marker='v', label='')
#plt.scatter( exp_010[:,0], exp_010[:,2], edgecolor='b', facecolor='None', marker='^', label='')
#plt.scatter( exp_018[:,0], exp_018[:,1], edgecolor='r', facecolor='None', marker='v', label='')
#plt.scatter( exp_018[:,0], exp_018[:,2], edgecolor='r', facecolor='None', marker='^', label='')

plt.plot( ECW_023_NES_heli[:,0], ECW_023_NES_heli[:,1], c='r', linestyle='--' )
plt.plot( ECW_023_NES_cart[:,0], ECW_023_NES_cart[:,1], c='r', linestyle='--' )

plt.scatter( ECW_015[:,0], ECW_015[:,1], edgecolor='b', facecolor='None', marker='v', label='' )
plt.scatter( ECW_015[:,0], ECW_015[:,2], edgecolor='b', facecolor='None', marker='^', label='' )
plt.scatter( ECW_023[:,0], ECW_023[:,1], edgecolor='r', facecolor='None', marker='v', label='' )
plt.scatter( ECW_023[:,0], ECW_023[:,2], edgecolor='r', facecolor='None', marker='^', label='' )

plt.plot( xaxis, S0115_heli_NES, color='r', linestyle='-' )
plt.plot( xaxis, S0115_cart_NES, color='r', linestyle='-' )
plt.plot( xaxis, S0085_heli_NES, color='b', linestyle='-' )
plt.plot( xaxis, S0085_cart_NES, color='b', linestyle='-' )

plt.scatter( HSE03_13X_0115[:,0], HSE03_13X_0115[:,4], c='r', marker='v', label='helicopter' )
plt.scatter( HSE03_13X_0115[:,0], HSE03_13X_0115_cartwheel[:], c='r', marker='^', label='cartwheel' )
plt.scatter( HSE03_13X_0085[:,0], HSE03_13X_0085[:,4], c='b', marker='v', label='' )
plt.scatter( HSE03_13X_0085[:,0], HSE03_13X_0085_cartwheel[:], c='b', marker='^', label='' )

plt.annotate('(b)', xy=(43, 0.03), size=12)

#plt.legend(loc='upper right', numpoints=1, scatterpoints=1)
plt.xlim(0.,50.)
plt.ylim(0.,0.47)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel('Reaction probability')
plt.xlabel('Incidence angle (deg.)')

plt.tight_layout()
plt.subplots_adjust(wspace=0, hspace=0)
plt.savefig('reactionprobability_angle.pdf')
#plt.show()