Here is a notebook for homogeneous gas model.
Here we are talking about a homogeneous gas bulk of neutrinos with single energy. The EoM is $$ i \partial_t \rho_E = \left[ \frac{\delta m^2}{2E}B ,\rho_E \right] $$
while the EoM for antineutrinos is $$ i \partial_t \bar\rho_E = \left[- \frac{\delta m^2}{2E}B ,\bar\rho_E \right] $$
Initial: Homogeneous, Isotropic, Monoenergetic $\nu_e$ and $\bar\nu_e$
The equations becomes $$ i \partial_t \rho_E = \left[ \frac{\delta m^2}{2E} B ,\rho_E \right] $$ $$ i \partial_t \bar\rho_E = \left[- \frac{\delta m^2}{2E}B,\bar\rho_E \right] $$
Define $\omega=\frac{\delta m^2}{2E}$, $\omega = \frac{\delta m^2}{-2E}$, $\mu=\sqrt{2}G_F n_\nu$ $$ i \partial_t \rho_E = \left[ \omega B ,\rho_E \right] $$ $$ i \partial_t \bar\rho_E = \left[\bar\omega B,\bar\rho_E \right] $$
where
$$ B = \frac{1}{2} \begin{pmatrix} -\cos 2\theta_v & \sin 2\theta_v \\ \sin 2\theta_v & \cos 2\theta_v \end{pmatrix} $$or just use theta =0.2rad
$$ L = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} $$Initial condition $$ \rho(t=0) = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} $$
$$ \bar\rho(t=0) =\begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} $$define the following quantities
%2. delm2E$= \delta m^2/2E$ %3. lamb $= \lambda$, lambb $= \bar\lambda$ %4. gF $= G_F$ %5. mu $=\mu$ 6. omega $=\omega$, omegab $=-\bar\omega$
# This line configures matplotlib to show figures embedded in the notebook,
# instead of opening a new window for each figure. More about that later.
# If you are using an old version of IPython, try using '%pylab inline' instead.
%matplotlib inline
%load_ext snakeviz
import numpy as np
from scipy.optimize import minimize
from scipy.special import expit
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import timeit
import pandas as pd
import plotly.plotly as py
from plotly.graph_objs import *
import plotly.tools as tls
# hbar=1.054571726*10**(-34)
hbar=1.0
delm2E=1.0
# lamb=1.0 ## lambda for neutrinos
# lambb=1.0 ## lambda for anti neutrinos
# gF=1.0
# nd=1.0 ## number density
# ndb=1.0 ## number density
omega=1.0
omegab=-1.0
## Here are some matrices to be used
elM = np.array([[1.0,0.0],[0.0,0.0]])
#bM = 1.0/2*np.array( [ [ - 0.38729833462,0.31622776601] , [0.31622776601,0.38729833462] ] )
bM = 1.0/2*np.array( [ [ -np.cos(0.2),np.sin(0.2)] , [np.sin(0.2),np.cos(0.2)] ] )
print bM
## sqareroot of 2
sqrt2=np.sqrt(2.0)
[[-0.49003329 0.09933467] [ 0.09933467 0.49003329]]
I am going to substitute all density matrix elements using their corrosponding network expressions.
So first of all, I need the network expression for the unknown functions.
A function is written as
$$ y_i= initial+t_i v_k f(t_i w_k+u_k) ,$$while it's derivative is
$$v_k f(t w_k+u_k) + t v_k f(tw_k+u_k) (1-f(tw_k+u_k)) w_k .$$Now I can write down the equations using these two forms.
def trigf(x):
#return 1/(1+np.exp(-x)) # It's not bad to define this function here for people could use other functions other than expit(x).
return expit(x)
rho = [[r11, a+ib ], [a-ib,r22]], the following rho4 function is rho4 = [r11,a, b, r22]
## The time derivative part
### Here are the initial conditions
init = np.array( [1.0, 0.0, 0.0, 0.0] )
### For neutrinos
def rho4(x,ti,initialCondition): # x is the input structure arrays, ti is a time point
elem0 = np.sum(ti * x[0] * trigf( ti*x[1] + x[2] ) )
elem1 = np.sum(ti * x[3] * trigf( ti*x[4] + x[5] ) )
elem2 = np.sum(ti * x[6] * trigf( ti*x[7] + x[8] ) )
elem3 = np.sum(ti * x[9] * trigf( ti*x[10] + x[11] ) )
return initialCondition + np.array([ elem0 , elem1, elem2, elem3 ])
## Test
xtemp=np.ones(12)
rho4(xtemp,1,init)
array([ 1.88079708, 0.88079708, 0.88079708, 0.88079708])
## Define Hamiltonians for both
hamilv = delm2E*bM.flatten()
print hamilv
[-0.49003329 0.09933467 0.09933467 0.49003329]
## Test
print bM
print hamilv
[[-0.49003329 0.09933467] [ 0.09933467 0.49003329]] [-0.49003329 0.09933467 0.09933467 0.49003329]
A function is written as
$$ y_i= initial+t_i v_k f(t_i w_k+u_k) ,$$while it's derivative is
$$v_k f(t w_k+u_k) + t v_k f(tw_k+u_k) (1-f(tw_k+u_k)) w_k .$$## The COST of the eqn set
#remember that commutator is[1.0j* e0, e1 + 1.0j*e2 , - e1 + 1.0j*e2 ,1.0j*e3 ], where exs are the elements of the commv
regularization = 0.0001
npsum = np.sum
def costvTi(x,ti,initialCondition): # l is total length of x
list = np.linspace(0,3,4)
fvec = []
rho4i = rho4(x,ti,initialCondition)
a0 = 2.0*rho4i[2]*hamilv[1]
a1 = -2.0*rho4i[2]*hamilv[0]
a2 = 2.0*rho4i[1]*hamilv[0] + hamilv[1]*( rho4i[3] - rho4i[0] )
a3 = -2.0*rho4i[2]*hamilv[2]
added = np.array([a0,a1,a2,a3])
fvecappend = fvec.append
for i in list:
fvecappend(np.asarray(trigf(ti*1.0*x[i*3+1] + 1.0*x[i*3+2]) ) )
fvec = np.array(fvec)
rhoprime = np.zeros(4)
for i in list:
rhoprime[i] = ( np.sum (1.0*x[i*3]*fvec[i] + 1.0*ti * x[i*3]* fvec[i] * ( 1.0 - fvec[i] ) * x[i*3+1] ) )
costi = rhoprime + added
costiTemp = np.sum(costi**2)
return costiTemp
print costvTi(xtemp,2,init)
11.7835898849
## Calculate the total cost
def costv(x,t,initialCondition):
t = np.array(t)
costvTotal = np.sum( costvTList(x,t,initialCondition) )
return costvTotal
def costvTList(x,t,initialCondition): ## This is the function WITHOUT the square!!!
t = np.array(t)
costvList = np.asarray([])
for temp in t:
tempElement = costvTi(x,temp,initialCondition)
costvList = np.append(costvList, tempElement)
return np.array(costvList)
ttemp = np.linspace(0,10)
print ttemp
[ 0. 0.20408163 0.40816327 0.6122449 0.81632653 1.02040816 1.2244898 1.42857143 1.63265306 1.83673469 2.04081633 2.24489796 2.44897959 2.65306122 2.85714286 3.06122449 3.26530612 3.46938776 3.67346939 3.87755102 4.08163265 4.28571429 4.48979592 4.69387755 4.89795918 5.10204082 5.30612245 5.51020408 5.71428571 5.91836735 6.12244898 6.32653061 6.53061224 6.73469388 6.93877551 7.14285714 7.34693878 7.55102041 7.75510204 7.95918367 8.16326531 8.36734694 8.57142857 8.7755102 8.97959184 9.18367347 9.3877551 9.59183673 9.79591837 10. ]
ttemp = np.linspace(0,10)
print costvTList(xtemp,ttemp,init)
print costv(xtemp,ttemp,init)
[ 2.00241504 2.52486601 3.12957558 3.82314867 4.61583765 5.52007582 6.548817 7.71415413 9.0264291 10.49383832 12.12242401 13.91630864 15.87804541 18.00899333 20.30966179 22.77999828 25.41961232 28.22793954 31.20435522 34.34824788 37.65906304 41.13632579 44.77964879 48.58873118 52.56335186 56.70335997 61.00866425 65.47922231 70.11503074 74.91611611 79.88252726 85.01432875 90.31159553 95.77440856 101.40285145 107.19700781 113.15695936 119.28278448 125.57455727 132.03234697 138.65621761 145.4462279 152.4024313 159.52487616 166.81360596 174.26865961 181.89007178 189.67787322 197.63209112 205.75274943] 3592.25842934
Here is the minimization
endpoint = 10
tlin = np.linspace(0,endpoint,20)
# tlinTest = np.linspace(0,14,10) + 0.5
# initGuess = np.ones(120)
# initGuess = np.asarray(np.split(np.random.rand(1,720)[0],12))
initGuess = np.asarray(np.split(np.random.rand(1,60)[0],12))
costvF = lambda x: costv(x,tlin,init)
#costvFTest = lambda x: costv(x,tlinTest,init)
print costv(initGuess,tlin,init)#, costv(initGuess,tlinTest,init)
10483.2460234
#%%snakeviz
startSLSQP = timeit.default_timer()
costvFResultSLSQP = minimize(costvF,initGuess,method="SLSQP")
stopSLSQP = timeit.default_timer()
print stopSLSQP - startSLSQP
print costvFResultSLSQP
22.0414690971 status: 0 success: True njev: 44 nfev: 2745 fun: 0.029594802355524336 x: array([ -1.59252118e-03, 3.51014982e+00, 1.65251562e+00, -1.74724687e-01, -2.04199807e-01, -9.39817008e-01, 1.19918520e-01, -1.04374519e-01, -7.08835998e+00, 1.56467928e-03, 6.42148180e+00, 2.63351999e+00, 2.74194493e-01, 6.04660148e-01, 8.27430694e-01, 9.41868026e-01, 9.55032063e-01, 2.71249180e-01, 6.63824188e-01, 6.47188705e-01, 8.24388861e-01, 4.96686894e-01, 8.53520987e-01, 6.36443930e-02, 3.37829263e-01, 5.43065972e-02, 3.08361127e-01, 8.36703782e-01, 3.90895320e-01, 1.91154610e-02, 5.54548230e-02, 5.61579367e-01, 8.69265422e-01, 3.13275701e-01, 2.83225135e-02, 5.01541153e-01, 4.92425726e-01, 6.16425625e-01, 9.44842566e-02, 9.27639463e-01, 7.72354571e-01, 8.12257532e-01, 6.47191868e-01, 7.61174793e-01, 8.18460476e-01, 1.37300048e-01, 5.09062206e-01, 8.67575780e-01, 1.77022215e-02, 8.47841558e-01, 4.94487180e-01, 5.93556947e-01, 9.99917064e-01, 7.53008878e-01, 2.76437081e-01, 8.67706376e-01, 5.50429901e-01, 6.06363724e-01, 8.78098759e-01, 2.77274115e-01]) message: 'Optimization terminated successfully.' jac: array([ -1.17622549e-03, -3.59956175e-07, -8.19563866e-08, 3.20923305e-03, 4.90120612e-04, 5.18076122e-05, -5.90027310e-04, 8.09070189e-05, -7.06948340e-05, 1.03940256e-04, -2.79396772e-08, 2.23983079e-07, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]) nit: 44 *** Profile stats marshalled to file u'/var/folders/mj/1sl30v6x2g5_lnlgdnngtd3c0000gn/T/tmpbh57FC'.
#np.savetxt('./assets/homogen/optimize_ResultSLSQPT2120_Vac.txt', costvFResultSLSQP.get('x'), delimiter = ',')
Find the solutions to each elements.
# costvFResultSLSQPx = np.genfromtxt('./assets/homogen/optimize_ResultSLSQP.txt', delimiter = ',')
## The first element of neutrino density matrix
xresult = np.asarray(costvFResultSLSQP.get('x'))
#xresult = np.asarray(costvFResultNM.get('x'))
#xresult = np.asarray(costvFResultBFGS.get('x'))
print xresult
plttlin=np.linspace(0,endpoint,100)
pltdata11 = np.array([])
pltdata11Test = np.array([])
pltdata22 = np.array([])
pltdata12 = np.array([])
pltdata21 = np.array([])
for i in plttlin:
pltdata11 = np.append(pltdata11 ,rho4(xresult,i,init)[0] )
print pltdata11
for i in plttlin:
pltdata12 = np.append(pltdata12 ,rho4(xresult,i,init)[1] )
print pltdata12
for i in plttlin:
pltdata21 = np.append(pltdata21 ,rho4(xresult,i,init)[2] )
print pltdata21
#for i in plttlin:
# pltdata11Test = np.append(pltdata11Test ,rho(xresultTest,i,init)[0,0] )
#
#print pltdata11Test
for i in plttlin:
pltdata22 = np.append(pltdata22 ,rho4(xresult,i,init)[3] )
print pltdata22
print rho4(xresult,0,init)
[ -1.59252118e-03 3.51014982e+00 1.65251562e+00 -1.74724687e-01 -2.04199807e-01 -9.39817008e-01 1.19918520e-01 -1.04374519e-01 -7.08835998e+00 1.56467928e-03 6.42148180e+00 2.63351999e+00 2.74194493e-01 6.04660148e-01 8.27430694e-01 9.41868026e-01 9.55032063e-01 2.71249180e-01 6.63824188e-01 6.47188705e-01 8.24388861e-01 4.96686894e-01 8.53520987e-01 6.36443930e-02 3.37829263e-01 5.43065972e-02 3.08361127e-01 8.36703782e-01 3.90895320e-01 1.91154610e-02 5.54548230e-02 5.61579367e-01 8.69265422e-01 3.13275701e-01 2.83225135e-02 5.01541153e-01 4.92425726e-01 6.16425625e-01 9.44842566e-02 9.27639463e-01 7.72354571e-01 8.12257532e-01 6.47191868e-01 7.61174793e-01 8.18460476e-01 1.37300048e-01 5.09062206e-01 8.67575780e-01 1.77022215e-02 8.47841558e-01 4.94487180e-01 5.93556947e-01 9.99917064e-01 7.53008878e-01 2.76437081e-01 8.67706376e-01 5.50429901e-01 6.06363724e-01 8.78098759e-01 2.77274115e-01] [ 1. 0.9998582 0.99970599 0.99954735 0.99938508 0.99922104 0.99905637 0.99889172 0.99872741 0.99856357 0.99840023 0.99823737 0.99807491 0.99791279 0.99775096 0.99758935 0.99742792 0.99726663 0.99710544 0.99694434 0.9967833 0.9966223 0.99646134 0.99630041 0.99613949 0.99597859 0.9958177 0.99565682 0.99549594 0.99533507 0.9951742 0.99501333 0.99485247 0.9946916 0.99453074 0.99436988 0.99420902 0.99404816 0.99388729 0.99372643 0.99356557 0.99340471 0.99324385 0.99308299 0.99292213 0.99276127 0.99260041 0.99243955 0.99227869 0.99211782 0.99195696 0.9917961 0.99163524 0.99147438 0.99131352 0.99115266 0.9909918 0.99083094 0.99067008 0.99050922 0.99034836 0.9901875 0.99002664 0.98986577 0.98970491 0.98954405 0.98938319 0.98922233 0.98906147 0.98890061 0.98873975 0.98857889 0.98841803 0.98825717 0.98809631 0.98793545 0.98777458 0.98761372 0.98745286 0.987292 0.98713114 0.98697028 0.98680942 0.98664856 0.9864877 0.98632684 0.98616598 0.98600512 0.98584426 0.9856834 0.98552253 0.98536167 0.98520081 0.98503995 0.98487909 0.98471823 0.98455737 0.98439651 0.98423565 0.98407479] [ 0. -0.00488505 -0.00962502 -0.01422196 -0.01867792 -0.02299496 -0.02717519 -0.03122071 -0.03513365 -0.03891614 -0.04257033 -0.04609838 -0.04950246 -0.05278473 -0.05594737 -0.05899256 -0.06192248 -0.0647393 -0.06744519 -0.07004234 -0.0725329 -0.07491904 -0.0772029 -0.07938664 -0.08147237 -0.08346222 -0.08535831 -0.08716271 -0.08887752 -0.09050479 -0.09204656 -0.09350487 -0.09488172 -0.09617911 -0.09739899 -0.09854332 -0.09961403 -0.10061301 -0.10154215 -0.10240331 -0.1031983 -0.10392896 -0.10459704 -0.10520433 -0.10575253 -0.10624337 -0.10667851 -0.10705962 -0.10738832 -0.1076662 -0.10789484 -0.10807578 -0.10821053 -0.1083006 -0.10834743 -0.10835247 -0.10831712 -0.10824275 -0.10813073 -0.10798237 -0.10779897 -0.10758181 -0.10733212 -0.10705113 -0.10674002 -0.10639995 -0.10603207 -0.10563749 -0.10521729 -0.10477253 -0.10430424 -0.10381345 -0.10330112 -0.10276824 -0.10221572 -0.10164449 -0.10105543 -0.10044942 -0.09982729 -0.09918987 -0.09853795 -0.09787233 -0.09719374 -0.09650294 -0.09580063 -0.0950875 -0.09436424 -0.09363149 -0.0928899 -0.09214007 -0.09138262 -0.09061811 -0.08984711 -0.08907016 -0.0882878 -0.08750052 -0.08670884 -0.08591322 -0.08511412 -0.084312 ] [ 0.00000000e+00 9.99719455e-06 1.97848692e-05 2.93663146e-05 3.87447756e-05 4.79234516e-05 5.69054976e-05 6.56940242e-05 7.42920987e-05 8.27027454e-05 9.09289461e-05 9.89736409e-05 1.06839729e-04 1.14530067e-04 1.22047474e-04 1.29394728e-04 1.36574568e-04 1.43589694e-04 1.50442767e-04 1.57136414e-04 1.63673220e-04 1.70055737e-04 1.76286479e-04 1.82367925e-04 1.88302519e-04 1.94092669e-04 1.99740750e-04 2.05249102e-04 2.10620033e-04 2.15855817e-04 2.20958695e-04 2.25930878e-04 2.30774543e-04 2.35491838e-04 2.40084877e-04 2.44555746e-04 2.48906502e-04 2.53139168e-04 2.57255743e-04 2.61258194e-04 2.65148460e-04 2.68928452e-04 2.72600053e-04 2.76165121e-04 2.79625484e-04 2.82982944e-04 2.86239278e-04 2.89396237e-04 2.92455544e-04 2.95418901e-04 2.98287981e-04 3.01064434e-04 3.03749888e-04 3.06345943e-04 3.08854179e-04 3.11276151e-04 3.13613391e-04 3.15867409e-04 3.18039692e-04 3.20131707e-04 3.22144895e-04 3.24080681e-04 3.25940464e-04 3.27725626e-04 3.29437525e-04 3.31077502e-04 3.32646876e-04 3.34146947e-04 3.35578994e-04 3.36944278e-04 3.38244043e-04 3.39479511e-04 3.40651886e-04 3.41762357e-04 3.42812091e-04 3.43802240e-04 3.44733937e-04 3.45608299e-04 3.46426426e-04 3.47189400e-04 3.47898287e-04 3.48554137e-04 3.49157985e-04 3.49710849e-04 3.50213730e-04 3.50667616e-04 3.51073479e-04 3.51432275e-04 3.51744947e-04 3.52012422e-04 3.52235614e-04 3.52415421e-04 3.52552729e-04 3.52648408e-04 3.52703317e-04 3.52718299e-04 3.52694186e-04 3.52631795e-04 3.52531933e-04 3.52395390e-04] [ 0. 0.00015233 0.00031001 0.00046933 0.00062882 0.00078803 0.0009469 0.00110549 0.00126388 0.00142214 0.00158031 0.00173843 0.00189652 0.0020546 0.00221266 0.00237072 0.00252877 0.00268682 0.00284487 0.00300292 0.00316097 0.00331902 0.00347706 0.00363511 0.00379316 0.00395121 0.00410926 0.00426731 0.00442536 0.0045834 0.00474145 0.0048995 0.00505755 0.0052156 0.00537365 0.00553169 0.00568974 0.00584779 0.00600584 0.00616389 0.00632194 0.00647998 0.00663803 0.00679608 0.00695413 0.00711218 0.00727023 0.00742828 0.00758632 0.00774437 0.00790242 0.00806047 0.00821852 0.00837657 0.00853461 0.00869266 0.00885071 0.00900876 0.00916681 0.00932486 0.0094829 0.00964095 0.009799 0.00995705 0.0101151 0.01027315 0.0104312 0.01058924 0.01074729 0.01090534 0.01106339 0.01122144 0.01137949 0.01153753 0.01169558 0.01185363 0.01201168 0.01216973 0.01232778 0.01248582 0.01264387 0.01280192 0.01295997 0.01311802 0.01327607 0.01343412 0.01359216 0.01375021 0.01390826 0.01406631 0.01422436 0.01438241 0.01454045 0.0146985 0.01485655 0.0150146 0.01517265 0.0153307 0.01548874 0.01564679] [ 1. 0. 0. 0.]
rho4(xresult,1,init)
array([ 9.98416546e-01, -4.22106245e-02, 9.01145420e-05, 1.56449654e-03])
#np.savetxt('./assets/homogen/optimize_pltdatar11.txt', pltdata11, delimiter = ',')
#np.savetxt('./assets/homogen/optimize_pltdatar22.txt', pltdata22, delimiter = ',')
plt.figure(figsize=(16,9.36))
plt.ylabel('rho11')
plt.xlabel('Time')
plt11=plt.plot(plttlin,pltdata11,"b4-",label="vac_rho11")
#plt.plot(plttlin,pltdata11Test,"m4-",label="vac_rho11Test")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="vac_HG-rho11")
# tls.embed("https://plot.ly/~emptymalei/73/")
plt.figure(figsize=(16,9.36))
plt.ylabel('rho12')
plt.xlabel('Time')
plt12=plt.plot(plttlin,pltdata12,"b4-",label="vac_rho12")
#plt.plot(plttlin,pltdata11Test,"m4-",label="vac_rho11Test")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="vac_HG-rho11")
# tls.embed("https://plot.ly/~emptymalei/73/")
plt.figure(figsize=(16,9.36))
plt.ylabel('rho21')
plt.xlabel('Time')
plt11=plt.plot(plttlin,pltdata21,"b4-",label="vac_rho21")
#plt.plot(plttlin,pltdata11Test,"m4-",label="vac_rho11Test")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="vac_HG-rho11")
# tls.embed("https://plot.ly/~emptymalei/73/")
plt.figure(figsize=(16,9.36))
plt.ylabel('Time')
plt.xlabel('rho22')
plt22=plt.plot(plttlin,pltdata22,"r4-",label="vac_rho22")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="vac_HG-rho22")
MMA_optmize_Vac_pltdata = np.genfromtxt('./assets/homogen/MMA_optmize_Vac_pltdata.txt', delimiter = ',')
plt.figure(figsize=(16,9.36))
plt.ylabel('MMArho11')
plt.xlabel('Time')
plt.plot(np.linspace(0,15,4501),MMA_optmize_Vac_pltdata,"r-",label="MMAVacrho11")
plt.plot(plttlin,pltdata11,"b4-",label="vac_rho11")
plt.show()
#py.iplot_mpl(plt.gcf(),filename="MMA-rho11-Vac-80-60")
xtemp1 = np.arange(4)
xtemp1.shape = (2,2)
print xtemp1
xtemp1[0,1]
np.dot(xtemp1,xtemp1)
xtemp1[0,1]