#!/usr/bin/env python

import numpy as np
import math
import matplotlib.pyplot as plt
from   pylab import *

## Create two plots, one for the average density and one for the average diameter
fig1 = plt.figure()
denPlot = plt.subplot(211)
diamPlot = plt.subplot(212)

## Open the data file generated by Xolotl (default is Alloy.dat)
time1, vDen1, vDiam1, voidDen1, voidDiam1, faultedDen1, faultedDiam1, iDen1, iDiam1, perfectDen1, perfectDiam1, frankDen1, frankDiam1 = loadtxt('/home/sophie/Workspace/xolotl-plsm-build/script/Alloy.dat', usecols = (1,2,3,6,7,10,11,14,15,18,19,22,23) , unpack=True) 

## Density data
denPlot.plot(time1, iDen1 * 1.0e27, lw=15, ls='-', color=plt.cm.jet(100), alpha=0.7, label='I')
denPlot.plot(time1, vDen1 * 1.0e27, lw=13, ls='-', color=plt.cm.jet(180), alpha=0.7, label='V')
denPlot.plot(time1, voidDen1 * 1.0e27, lw=11, ls='-', color=plt.cm.jet(225), alpha=0.7, label='void')
denPlot.plot(time1, faultedDen1 * 1.0e27, lw=9, ls='-', color=plt.cm.jet(150), alpha=0.7, label='faulted')
denPlot.plot(time1, perfectDen1 * 1.0e27, lw=7, ls='-', color=plt.cm.jet(80), alpha=0.7, label='perfect')
denPlot.plot(time1, frankDen1 * 1.0e27, lw=5, ls='-', color=plt.cm.jet(0), alpha=0.7, label='frank')

## Diameter data
diamPlot.plot(time1, iDiam1, lw=15, ls='-', color=plt.cm.jet(100), alpha=0.7, label='I')
diamPlot.plot(time1, vDiam1, lw=13, ls='-', color=plt.cm.jet(180), alpha=0.7, label='V')
diamPlot.plot(time1, voidDiam1, lw=11, ls='-', color=plt.cm.jet(225), alpha=0.7, label='void')
diamPlot.plot(time1, faultedDiam1, lw=9, ls='-', color=plt.cm.jet(150), alpha=0.7, label='faulted')
diamPlot.plot(time1, perfectDiam1, lw=7, ls='-', color=plt.cm.jet(80), alpha=0.7, label='perfect')
diamPlot.plot(time1, frankDiam1, lw=5, ls='-', color=plt.cm.jet(0), alpha=0.7, label='frank')

## Plot the legend
l = denPlot.legend(loc='best')
setp(l.get_texts(), fontsize=20)

## Some shaping
time = time1[len(time1)-1]
denPlot.set_xlabel("Time (s)",fontsize=25)
denPlot.set_ylabel("Density (m-3)",fontsize=25)
denPlot.set_xlim([1.0e-1, time])
#denPlot.set_ylim([1.0e11, 1.0e27])
denPlot.set_xscale('log')
denPlot.set_yscale('log')
denPlot.grid()
denPlot.tick_params(axis='both', which='major', labelsize=25)
denPlot.tick_params(axis='both', which='minor', labelsize=25)
diamPlot.set_xlabel("Time (s)",fontsize=25)
diamPlot.set_ylabel("Diameter (nm)",fontsize=25)
diamPlot.set_xlim([1.0e-1, time])
#diamPlot.set_ylim([0.0, 4.0])
diamPlot.set_xscale('log')
diamPlot.set_yscale('log')
diamPlot.grid()
diamPlot.tick_params(axis='both', which='major', labelsize=25)
diamPlot.tick_params(axis='both', which='minor', labelsize=25)

## Show the plots
plt.show()
