# Plots all the clouds on the Tully-Fisher relation

import os
import matplotlib
import matplotlib.pyplot as plt
import numpy
import math
import random
from random import random as rand
from math import *


# LaTeX like fonts in the figure
#plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.rc('font', serif='Palatino')
plt.rc('font', cursive='Zapf Chancery')
plt.rc('font', monospace='Courier')

matplotlib.rcParams['xtick.labelsize'] = 14 
matplotlib.rcParams['ytick.labelsize'] = 14

from matplotlib.ticker import AutoMinorLocator
XminorLocator = AutoMinorLocator()
YminorLocator = AutoMinorLocator()


abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)


# Add the simulated dark galaxies
#dgdir = '/home/rhys/Software/Gf/IsoDiscCluster/8.65kpc/1500K/StreamParameterFiles_AGES_test'
#os.chdir(dgdir)

#files = os.listdir(dgdir)

#goodfiles = []

#for eachfile in files:
#	if 'StreamEvolution' in str(eachfile):
#		goodfiles.append(eachfile)

#for eachfile in goodfiles:
#	eachfile = str(eachfile)
			
#	data= numpy.genfromtxt(eachfile, skip_header=1)
	
	#frame = data[:,0]
	#nc		= data[:,1]
	#peaksn= data[:,2]
#	mass	= data[:,6]
#	width	= data[:,9] + (20.0*rand() - 10.0) #+ (50.0*rand() - 10.0)
	
#	plt.figure(1)
#	plt.scatter(numpy.log10(width/2.0),numpy.log10(1.4*mass), marker='.', s=0.01, color='black')

#os.chdir(dname)



# Add the simulated clouds
files = os.listdir(dname)

goodfiles = []

for eachfile in files:
	if 'AllCloudEvolution' in str(eachfile):
		goodfiles.append(eachfile)

for eachfile in goodfiles:
	eachfile = str(eachfile)
			
	data= numpy.genfromtxt(eachfile, skip_header=1)
	
	frame = data[:,0]
	nc		= data[:,1]
	peaksn= data[:,2]
	mass	= data[:,3]
	width	= data[:,4]
	
	plt.figure(1)
	#plt.scatter(numpy.log10(width/2.0),numpy.log10(1.4*mass), marker='^', s=30.0, color='green')


# Add the median dark discs
data = numpy.genfromtxt('DarkGalaxyEvolution.txt',skip_header=1)
m = data[:,0]
v = data[:,1]

#plt.scatter(numpy.log10(v/2.0),numpy.log10(1.4*m), s=40.0, color='grey')

# Add the initial position
m1 = data[0,0]
v1 = data[0,1]
#plt.scatter(numpy.log10(v1/2.0),numpy.log10(1.4*m1), marker='o', s=50.0, edgecolor='yellow', facecolors='none')

# Add the final position
m2 = data[len(data)-1,0]
v2 = data[len(data)-1,1]
#plt.scatter(numpy.log10(v2/2.0),numpy.log10(1.4*m2), marker='s', s=50.0, edgecolor='pink',facecolor='none')




# Add the background Virgo objects
data = numpy.genfromtxt('BackgroundIDL.txt')
v = data[:,0]
m = data[:,1]

plt.scatter(v,m*0.975, s=30.0, linewidth=0.0)


# Add the low-deficiency Virgo objects
data = numpy.genfromtxt('VirgoLDIDL.txt')
v = data[:,0]
m = data[:,1]

plt.scatter(v,m*0.975, s=30.0, linewidth=0.0)


# Add the (optically faint) non-VCC objects
data = numpy.genfromtxt('VirgoNonVCC.txt', skip_header=1)
v = data[:,1]
m = data[:,2]

#plt.scatter(v,m*0.975, marker='s', edgecolor='red',facecolor='none', s=50.0, linewidth=1.0)



# Add the Virgo dark clouds
xvirgo = numpy.array([152.0, 33.0, 157.0, 120.0, 146.0, 173.0, 35.0, 164.0])
yvirgo = numpy.array([4.2E7, 2.3E7,1.4E7, 1.4E7, 2.0E7, 3.2E7, 7.3E6,4.4E7])

plt.scatter(numpy.log10(xvirgo/2.0),numpy.log10(yvirgo*1.4), s=40.0, marker='s', color='black', linewidth=0.0)


# Add the ALFALFA UDGs
data = numpy.genfromtxt('ALFALFAUDGs.txt',skip_header=1)
v = data[:,0]
m = data[:,1]

# Correct for helium to be compatible with existing obs
for i in range(0,len(m)):
	m[i] = log10((10.0**m[i]) * 1.4)

plt.scatter(numpy.log10(v/(2.0*sin(radians(45.0)))), m, s=40.0, marker='s', color='green', linewidth=0.0)

	
# Get data range to set axes for a square plot
#xmin = 1.2	#0.5	# min(plt.xlim())
#xmax = 2.6	#2.25	#max(plt.xlim())
#ymin = 6.8	#6.25	#min(plt.ylim())
#ymax = 11.0	#8.5	#max(plt.ylim())

xmin = 1.1
xmax = 2.6
ymin = 6.5
ymax = 11.0

xlim = xmax - xmin
ylim = ymax - ymin


# Add the baryonic Tully Fisher from McGaugh 2000
Vw = [xmin, xmax]
Mr = [ymin, ymax]
plt.plot(Vw,Mr, color='r', linewidth=1.0)


# *** McGaugh 2000 relation ****
# M = 35.0*V_rot^4
# M in (linear) solar masses when V_rot in km/s
# But the best-fit relation used above appears to be M = 10^3.2 * V_rot^3


# Set aspect ratio
plt.gca().set_aspect(xlim/ylim, adjustable='box')

axes = plt.gca()
axes.set_xlim([xmin,xmax])
axes.set_ylim([ymin,ymax])

axes.xaxis.set_minor_locator(XminorLocator)
axes.yaxis.set_minor_locator(YminorLocator)


# Set font
csfont = {'fontname':'Liberation Serif','fontsize':15}
plt.xlabel('log(W20) / km s$^{-1}$', **csfont)
plt.ylabel('log ($\mathregular{M_{baryons}}$) / M$_{\odot}$' , **csfont)

	
plt.savefig('TEST.png',dpi=100,facecolor='white',bbox_inches='tight')