from __future__ import division # this is to avoid python 2.x integer division issues

"""This is the Python code and datasets to generate the figures 1 abd 5 of the manuscript entitled:
"An example of low Th/U zircon overgrowths of magmatic origin in a late orogenic Variscan
intrusive: the San Ciprián massif (NW Spain)" by Marco A. Lopez-Sanchez et al. (2015)"""

"""
Requirements:
    Python version 2.7.x or 3.4.x or higher
    Numpy version 1.5 or higher
    Matplotlib version 1.4.2 or higher
""" 

from numpy import array
import matplotlib.pyplot as plt # import the plot library

import matplotlib as mpl
mpl.style.use('ggplot')
# mpl.rcParams['font.family'] = 'Helvetica Neue'
mpl.rcParams['axes.labelsize'] = 15.
mpl.rcParams['xtick.labelsize'] = 15.
mpl.rcParams['ytick.labelsize'] = 15.

####################################################################################################################################################################3
### DATA SETS

# San Ciprian granite

UcoreDGP09 = array([364, 449, 298, 504, 503, 365, 267, 378, 449, 491, 198, 236, 173, 343, 252, 253, 308, 385, 189, 244, 343])
ThcoreDGP09 = array([158, 136, 123, 217, 239, 142, 84, 170, 155, 164, 77, 120, 48, 110, 89, 78, 109, 115, 56, 91, 156])

UrimDGP09 = array([3175, 2151, 1729])
ThrimDGP09 = array([127, 60, 45])

UcorePR12 = array([366, 244, 717, 288])
ThcorePr12 = array([155.92, 215.70, 185.70, 71.424])

UrimPR12 = array([934, 1464, 801, 951, 2426, 1808, 1791, 2018])
ThrimPR12 = array([53.24, 45.38, 16.82, 116.97, 33.96, 77.74, 84.18, 42.38])

UuztPR12 = array([2251, 867, 4874, 5984, 7166])
ThuztPR12 = array([49.52, 37.28, 146.22, 5.98, 100.32])

### OTHER DATASETS

# from Zeck and Whitehouse (1999) Contrib. Min. Pet. 134: 134-149 (only zircon rims)

"""Torrox gneiss; orthogneiss derived from a megacrystic plutonic complex; peraluminous"""
U_zeck = array([3732, 3057, 1659, 1988, 1670, 5344, 2302, 8870, 3604, 5833, 7039, 7214, 1349, 1308, 6315, 3850, 5543, 7798, 7188, 1693, 17614, 1462, 39364, 8378, 13745, 9854, 2167, 1112, 1410])
Th_zeck = array([74.64, 30.57, 16.59, 19.88, 33.4, 53.44, 46.04, 177.4, 72.08, 116.66, 140.78, 144.28, 13.49, 13.08, 126.3, 38.5, 55.43, 77.98, 71.88, 16.93, 880.7, 14.62, 1180.92, 167.56, 137.45, 98.54, 21.67, 11.12, 14.10])

# from Lopez-Sanchez et al. 2014. Peraluminous granites.

"""granite dykes Portomarin - Sarria leucogranite (strongly peraluminous)"""
U_DGP08 = array([5178, 8818, 3556, 7449, 13527, 5073, 4163, 7634, 28273, 4813, 3000, 6321, 7475, 4761, 11187, 5492, 4016, 144381, 6305, 3426, 136830, 4077, 83559])
Th_DGP08 = array([300, 2645, 73, 190, 946, 204, 76, 109, 997, 124, 51, 249, 812, 150, 231, 244, 338, 4492, 342, 363, 4486, 146, 1908])

"""Viveiro monzogranite (slighly peraluminous; Si content 59-70 %; Th/U = 3.32 average of four samples); discarded older (i.e. presumably inherited) ages"""
U_DGP03 = array([494, 1541, 874, 794, 743, 1207, 323, 842, 521, 1256, 590, 704, 821, 1235, 2005, 473, 1104, 586, 436, 732, 1510])
Th_DGP03 = array([160, 390, 279, 332, 1032, 670, 174, 1218, 390, 905, 643, 386, 839, 789, 262, 16, 624, 361, 349, 679, 859])

WR_ThUratio_DGP03 = 3.32

"""Penedo Gordo granite (strongly peraluminous low-Ca; Si content 70-75 %)"""
U_DGP04 = array([413, 296, 242, 944, 302, 170, 172, 149, 1117, 74, 155, 103, 1947, 1499, 271, 294, 190, 354, 142, 280, 173, 645, 601, 2118, 262, 536, 213, 124, 98, 113, 107, 1540, 170, 513, 162])
Th_DGP04 = array([437, 1035, 115, 187, 157, 91, 73, 170, 246, 62, 97, 76, 337, 170, 98, 134, 89, 132, 127, 111, 84, 193, 185, 490, 130, 156, 101, 83, 88, 56, 85, 160, 215, 151, 111])

### Data from Belousova et al. (2002) average Th and U values in zircons from different igneous rocks

"""Kimberlite"""
U_kimb = [9.8]
Th_kimb = [3.6]
U_upperLimit_kimb = [59.2]
U_lowerLimit_kimb = [6.8]
Th_upperLimit_kimb = [44.4]
Th_lowerLimit_kimb = [2.76]

"""Lamproite"""
U_lamp = [74]
Th_lamp = [64]
U_upperLimit_lamp = [274]
U_lowerLimit_lamp = [57]
Th_upperLimit_lamp = [124]
Th_lowerLimit_lamp = [52]

"""Carbonatite"""
U_carb = [0.29]
Th_carb = [212]
U_upperLimit_carb = [130.71]
U_lowerLimit_carb = [0.26]
Th_upperLimit_carb = [633]
Th_lowerLimit_carb = [209.1]

"""Basalt"""
U_basalt = [106]
Th_basalt = [64]
U_upperLimit_basalt = [1341]
U_lowerLimit_basalt = [89]
Th_upperLimit_basalt = [3106]
Th_lowerLimit_basalt = [56]

"""Dolerite"""
U_dol = [274]
Th_dol = [331]
U_upperLimit_dol = [260]
U_lowerLimit_dol = [133]
Th_upperLimit_dol = [320]
Th_lowerLimit_dol = [236]

"""Granitoid"""
U_gr = [764]
Th_gr = [368]
U_upperLimit_gr = [11617]
U_lowerLimit_gr = [688]
Th_upperLimit_gr = [11720]
Th_lowerLimit_gr = [337]

"""Syenite"""
U_syen = [88]
Th_syen = [31]
U_upperLimit_syen = [16]
U_lowerLimit_syen = [17]
Th_upperLimit_syen = [4]
Th_lowerLimit_syen = [7]

"""Larvikite"""
U_larv = [662]
Th_larv = [595]
U_upperLimit_larv = [788]
U_lowerLimit_larv = [358]
Th_upperLimit_larv = [858]
Th_lowerLimit_larv = [337]

"""Syenite pegmatite"""
U_syPeg = [356]
Th_syPeg = [647]
U_upperLimit_syPeg = [658]
U_lowerLimit_syPeg = [352.4]
Th_upperLimit_syPeg = [1871]
Th_lowerLimit_syPeg = [528]

"""Nepheline-syenite pegmatite"""
U_nePeg = [9.0]
Th_nePeg = [2497]
U_upperLimit_nePeg = [12]
U_lowerLimit_nePeg = [7]
Th_upperLimit_nePeg = [8587]
Th_lowerLimit_nePeg = [1868]

### Data from Kirkland et al. (2015). Average Th and U values in zircons from different igneous rocks

"""Alkali granite"""
U_alGR = [145.27] # Tukey's biweight values instead of means
Th_alGR = [111.07]
U_upperLimit_alGR = [1005.82]
U_lowerLimit_alGR = [18.00]
Th_upperLimit_alGR = [906.00]
Th_lowerLimit_alGR = [14.84]
ratio_alGR = [0.72]
IQ_upper_alGR = [0.945] # IQ means Interquartil
IQ_lower_alGR = [0.535]

"""Granite"""
U_GR = [193.46]
Th_GR = [138.93]
U_upperLimit_GR = [6485.00]
U_lowerLimit_GR = [11.00]
Th_upperLimit_GR = [4413.11]
Th_lowerLimit_GR = [4.79]
ratio_GR = [0.74]
IQ_upper_GR = [1.025]
IQ_lower_GR = [0.475]

### Th-U data from Rubatto (2002), table 1. Chemical Geology 184:123-138

"""Metamorphic zircon (granulitic facies)"""
Th_Rubatto = array([8.63, 14.8, 6.01, 10.3, 17.9, 8.21, 41.7, 9.41, 18.0, 11.3, 44.1, 58.0, 12.1, 1.23, 13.6, 1.37, 0.52, 48, 29.8, 183, 62.4, 275, 163, 80.8, 2.66, 1.99, 1.62, 0.918, 0.859])
U_Rubatto = array([712, 688, 415, 811, 634, 155, 581, 469, 441, 819, 1135, 1020, 124, 601, 250, 349, 166, 698, 619, 6477, 2858, 3157, 4300, 2707, 316, 192, 204, 94.9, 95.7])

### Th-U data from Vavra et al. (1999), table 4. Contrib Mineral Petrol 134:380-404

"""Metamorphic zircon (amphibolite to granulite facies)"""
Th_Vavra_STR1 = array([4, 8, 5, 10, 11, 9, 10, 8, 8, 2, 3, 10, 11, 5, 4, 4, 4, 2, 1, 3, 1])
U_Vavra_STR1 = array([1008, 2149, 1428, 2243, 2154, 1961, 2141, 2022, 2279, 527, 793, 2119, 2257, 1289, 721, 676, 905, 579, 487, 660, 427])

Th_Vavra_STR2 = array([10, 13, 10, 11, 13, 29, 13, 8, 12, 40, 18, 20, 22, 12, 27, 15, 17, 48, 23])
U_Vavra_STR2 = array([216, 302, 297, 300, 309, 566, 299, 384, 308, 718, 501, 416, 699, 291, 393, 388, 377, 1900, 370])

Th_Vavra_STR3 = array([37, 22, 33, 54, 53, 40, 60, 542, 281, 113, 166, 55, 38, 181, 153, 104, 116, 76, 141, 107])
U_Vavra_STR3 = array([675, 579, 674, 690, 844, 468, 699, 458, 292, 224, 224, 208, 153, 303, 467, 267, 350, 238, 255, 215])

Th_Vavra_RIM = array([18, 22, 39, 39, 9, 15, 10, 18, 30])
U_Vavra_RIM = array([921, 469, 732, 666, 364, 369, 607, 469, 663])

### Th-U data from Schaltegger et al. (1999), table 1. Contrib Mineral Petrol 134:186-201.

"""Metamorphic zircon high grade facies"""
Th_Schal = array([13, 177, 237, 26, 59, 10, 44, 10, 14])
U_Schal = array([237, 1387, 1612, 232, 447, 323, 252, 381, 188])

### Th-U data from Zhu et al. (2009), table 1. Tectonophysics 469:48-60

"""Igneous zircon in a Permian strongly peraluminous granite (Lhasa terrane, southern Tibet); Si content 73.37, Th/U = 6)"""
Th_Zhu = array([65, 75, 81, 56, 65, 45, 79, 40, 99, 71, 272, 152, 93, 81, 53, 55, 46, 53, 346, 60, 122, 106, 107, 738, 105, 850, 88, 552, 96, 106, 159, 59, 57, 123, 81, 47, 164, 129, 50, 255, 58, 132, 64, 64, 122, 54, 56, 55, 72, 43, 63, 83, 70, 58, 407, 72, 142, 176, 110])
U_Zhu = array([214, 216, 195, 168, 279, 324, 90, 408, 160, 165, 403, 284, 200, 285, 149, 626, 327, 335, 172, 93, 572, 389, 185, 898, 296, 951, 308, 956, 416, 243, 229, 122, 346, 239, 188, 283, 204, 303, 191, 591, 223, 185, 394, 208, 262, 232, 298, 687, 298, 232, 191, 146, 183, 428, 388, 178, 313, 379, 192])

WR_ThUratio_Zhu = 20.3/3.39 # see table 4 sample PK01-1

### Th-U data from Ladenberger et al. (2014), table 2. Geol. Soc. London, Special publications 390. doi: 10.1144/SP390.10

"""sample YBS08-6 Small pegmatite dyke with low Th/U"""
Th_YBS08_6 = array([536, 108, 658, 8, 789, 791, 701, 701, 1552, 371, 389, 569, 670, 682, 805])
U_YBS08_6 = array([12632, 9206, 18751, 1628, 18030, 17589, 14675, 30148, 16011, 9374, 13958, 14487, 14716, 16493])

"""sample YBS08-7 Large pegmatitic dyke with low Th/U"""
Th_YBS08_7 = array([48, 303, 147, 33, 67, 38, 11, 126, 10, 49, 37, 40, 49, 238, 19, 80, 87, 108, 63, 9, 15])
U_YBS08_7 = array([5088, 20044, 3936, 3891, 3787, 9471, 3349, 7793, 5316, 5715, 2014, 7747, 7253, 13329, 2931, 10122, 12451, 9589, 11953, 2902, 4356])

"""sample AL07_7 Leucogranite, bright igneous rims"""
Th_AL07_7 = array([48, 50, 68, 27, 20])
U_AL07_7 = array([272, 239, 600, 199, 249])

"""sample AL07_9 Leucogranite, bright igneous rims"""
Th_AL07_9 = array([17, 18, 21, 18, 21, 15, 18, 17, 11, 26, 15, 16])
U_AL07_9 = array([255, 289, 247, 218, 284, 230, 321, 284, 524, 437, 257, 287])

### Th-U data from Scharer (1984), table 1. Earth and Plannetary Sci. Letters 67:191-204.

"""sample A-141 Makalu leucogranite, three groups of zircons based on their morphologies"""
Th_A141 = array([585.48, 174.35, 47.03])
U_A141 = array([2863, 8674, 4479])

##################################################################################################################################################################
### PLOTS

plt.figure(1, figsize=(15, 5))
plt.gcf().subplots_adjust(bottom=0.15)
    
# Igneous samples Belousova et al. (2002) and Kirkland et al. (2015)
plt.subplot(121)
plt.semilogy()
plt.semilogx()
plt.plot([0.1, 15000], [0.1, 15000], 'k--', linewidth=1.5) # Th/U = 1
plt.plot([1, 15000], [0.1, 1500],'k-', linewidth=2) # Th/U = 0.1
plt.plot([10, 15000], [0.1, 150], 'k--', linewidth=1.5) # Th/U = 0.01
plt.plot([0.135, 15000], [0.1, 11100], 'b-', linewidth=2) # Th/U = 0.74
plt.plot([0.211, 15000], [0.1, 7125], 'b--', linewidth=1.5) # Th/U = 0.475 lower bound
plt.plot([0.098, 15000], [0.1, 15375], 'b--', linewidth=1.5) # Th/U = 1.025 upper bound
#plt.errorbar(U_kimb, Th_kimb, 
#             xerr=[U_lowerLimit_kimb, U_upperLimit_kimb], 
#             yerr=[Th_lowerLimit_kimb, Th_upperLimit_kimb], 
#             fmt='o', color='#8dd3c7',label='kimberlite', markersize=8)
#plt.errorbar(U_lamp, Th_lamp, 
#             xerr=[U_lowerLimit_lamp, U_upperLimit_lamp], 
#             yerr=[Th_lowerLimit_lamp, Th_upperLimit_lamp], 
#             fmt='o', color='#bebada', label='lamproite', markersize=8)
#plt.errorbar(U_carb, Th_carb, 
#             xerr=[U_lowerLimit_carb, U_upperLimit_carb], 
#             yerr=[Th_lowerLimit_carb, Th_upperLimit_carb], 
#             fmt='o', color='#ffed6f', label='carbonatite', markersize=6)
#plt.errorbar(U_basalt, Th_basalt, 
#             xerr=[U_lowerLimit_basalt, U_upperLimit_basalt], 
#             yerr=[Th_lowerLimit_basalt, Th_upperLimit_basalt], 
#             fmt='o', color='#fb8072', label='basalt', markersize=8)
#plt.errorbar(U_dol, Th_dol, 
#             xerr=[U_lowerLimit_dol, U_upperLimit_dol], 
#             yerr=[Th_lowerLimit_dol, Th_upperLimit_dol], 
#             fmt='o', color='#80b1d3', label='dolerite', markersize=8)
plt.errorbar(U_gr, Th_gr,
             xerr=[U_lowerLimit_gr, U_upperLimit_gr], 
             yerr=[Th_lowerLimit_gr, Th_upperLimit_gr], 
             fmt='o', color='#fdb462', label='granitoid', markersize=8)
plt.errorbar(U_syen, Th_syen, 
             xerr=[U_lowerLimit_syen, U_upperLimit_syen], 
             yerr=[Th_lowerLimit_syen, Th_upperLimit_syen], 
             fmt='^', color='#b3de69', label='syenite', markersize=8)
plt.errorbar(U_larv, Th_larv, 
             xerr=[U_lowerLimit_larv, U_upperLimit_larv], 
             yerr=[Th_lowerLimit_larv, Th_upperLimit_larv], 
             fmt='s', color='#fb8072', label='larvikite', markersize=8)
plt.errorbar(U_syPeg, Th_syPeg, 
             xerr=[U_lowerLimit_syPeg, U_upperLimit_syPeg], 
             yerr=[Th_lowerLimit_syPeg, Th_upperLimit_syPeg], 
             fmt='p', color='#bc80bd', label='syenite peg.', markersize=8)
#plt.errorbar(U_nePeg, Th_nePeg, 
#             xerr=[U_lowerLimit_nePeg, U_upperLimit_nePeg], 
#             yerr=[Th_lowerLimit_nePeg, Th_upperLimit_nePeg], 
#             fmt='o', color='#bc80bd', label='ne-sy pegmatite', markersize=6)
plt.errorbar(U_GR, Th_GR, 
             xerr=[U_lowerLimit_GR, U_upperLimit_GR], 
             yerr=[Th_lowerLimit_GR, Th_upperLimit_GR], 
             fmt='*', color='#ffed6f', label='granite', markersize=11)
plt.errorbar(U_alGR, Th_alGR, 
             xerr=[U_lowerLimit_alGR, U_upperLimit_alGR], 
             yerr=[Th_lowerLimit_alGR, Th_upperLimit_alGR], 
             fmt='H', color='#80b1d3', label='alkali granite', markersize=8)
plt.ylabel('Th (ppm)')
plt.xlabel('U (ppm)')
plt.legend(loc=2, fontsize=12, numpoints=1)

# Metamorphic zircons from Rubatto (2002), Vavra et al. (1999) and Schaltegger et al. (1999)
plt.subplot(122)
plt.semilogy()
plt.semilogx()
plt.plot([0.1, 15000], [0.1, 15000], 'k--', linewidth=1.5) # Th/U = 1
plt.plot([1, 15000], [0.1, 1500],'k-', linewidth=2) # Th/U = 0.1
plt.plot([10, 15000], [0.1, 150], 'k--', linewidth=1.5) # Th/U = 0.01
plt.plot(U_Schal, Th_Schal, 'H', color='#bc80bd', label='Schaltegger et al.', markersize=8)
plt.plot(U_Vavra_STR1, Th_Vavra_STR1, '^', color='#ffed6f', label='STR1', markersize=8)
plt.plot(U_Vavra_STR2, Th_Vavra_STR2, 's', color='#80b1d3', label='STR2', markersize=8)
plt.plot(U_Vavra_STR3, Th_Vavra_STR3, 'p', color='#fdb462', label='STR3', markersize=8)
plt.plot(U_Vavra_RIM, Th_Vavra_RIM, '*', color='#b3de69', label='RIM', markersize=11)
plt.plot(U_Rubatto, Th_Rubatto, 'o', color='#fb8072', label='Rubatto', markersize=8)
plt.xlabel('U (ppm)')
plt.legend(loc=2, fontsize=12, numpoints=1)


plt.figure(2) # San Ciprian sample and reference values
plt.gcf().subplots_adjust(bottom=0.15)
plt.semilogy()
plt.semilogx()
##plt.fill([10, 10000], [1, 1000], color='#CCCCCC')
plt.plot([10, 1351.35], [7.4, 1000], 'b-', linewidth=2) # Th/U = 0.74
plt.plot([10, 2105.26], [4.75, 1000], 'b--', linewidth=1.5) # Th/U = 0.475 lower bound
plt.plot([10, 975.61], [10.25, 1000], 'b--', linewidth=1.5) # Th/U = 1.025 upper bound
plt.plot([10, 1000], [10, 1000], 'k--', linewidth=1.5) # Th/U = 1
plt.plot([10, 10000], [1, 1000],'k-', linewidth=1.5) # Th/U = 0.1
plt.plot([100, 10000], [1, 100], 'k--', linewidth=1.5) # Th/U = 0.01
plt.plot(UcoreDGP09, ThcoreDGP09, 'o', color='#fb8072', label='zoned cores', markersize=9)
plt.plot(UrimDGP09, ThrimDGP09, 's', color='#fb8072', label='dark rims', markersize=9)
plt.plot(UcorePR12, ThcorePr12, 'o', color='#80b1d3', label='zoned cores', markersize=9)
plt.plot(UrimPR12, ThrimPR12, 's', color='#80b1d3', label='dark rims', markersize=9)
plt.plot(UuztPR12, ThuztPR12, '^', color='#80b1d3', label='unzoned tips', markersize=9)
plt.errorbar(U_GR, Th_GR, 
             xerr=[U_lowerLimit_GR, U_upperLimit_GR], 
             yerr=[Th_lowerLimit_GR, Th_upperLimit_GR], 
             fmt='o', color='#b3de69', label='granite', markersize=9)
plt.errorbar(U_alGR, Th_alGR, 
             xerr=[U_lowerLimit_alGR, U_upperLimit_alGR], 
             yerr=[Th_lowerLimit_alGR, Th_upperLimit_alGR], 
             fmt='o', color='#fccde5', label='alkali granite', markersize=9)
plt.errorbar(U_gr, Th_gr, 
             xerr=[U_lowerLimit_gr, U_upperLimit_gr], 
             yerr=[Th_lowerLimit_gr, Th_upperLimit_gr], 
             fmt='o', color='#fdb462', label='granitoid', markersize=9)
plt.ylabel('Th (ppm)')
##plt.xlabel('U (ppm)')
##plt.title ('Samples DGP-09 and PR12', fontsize=16, color='#262d35')
plt.legend(loc=2, fontsize=12, numpoints=1)  

plt.figure(3) # Other peraluminous granite samples from NW of Iberian massif and other places
plt.gcf().subplots_adjust(bottom=0.15)
plt.semilogy()
plt.semilogx()
plt.plot([10, 10000], [10, 10000], 'k--', linewidth=1.5) # Th/U = 1
plt.plot([10, 100000], [1, 10000],'k-', linewidth=2) # Th/U = 0.1
plt.plot([100, 100000], [1, 1000], 'k--', linewidth=1.5) # Th/U = 0.01
plt.plot(U_DGP03, Th_DGP03, 'o', color='#8dd3c7', label='sample DGP03', markersize=9)
plt.plot(U_DGP04, Th_DGP04, 'o', color='#ffed6f', label='sample DGP04', markersize=9)
plt.plot(U_Zhu, Th_Zhu, 'o', color='#80b1d3', label='sample PK01-1', markersize=9)
plt.plot(U_AL07_7, Th_AL07_7, 'o', color='#bebada', label='sample AL07_7', markersize=9)
plt.plot(U_AL07_9, Th_AL07_9, 's', color='#bc80bd', label='sample AL07_9', markersize=9)
plt.plot(UrimDGP09, ThrimDGP09, 's', color='#fdb462', label='San Ciprian', markersize=9)
plt.plot(UrimPR12, ThrimPR12, 's', color='#fdb462', markersize=9)
plt.plot(UrimDGP09, ThrimDGP09, 's', color='#fdb462', markersize=9)
plt.plot(UrimPR12, ThrimPR12, 's', color='#fdb462', markersize=9)
plt.plot(U_zeck, Th_zeck, 's', color='#fccde5', label='sample 90Z88', markersize=9)
plt.plot(U_DGP08, Th_DGP08, '^', color='#b3de69', label='sample DGP08', markersize=9)
plt.plot(U_A141, Th_A141, '^', color='#fb8072', label='sample A-141', markersize=9)
plt.errorbar(U_GR, Th_GR, 
             xerr=[U_lowerLimit_GR, U_upperLimit_GR], 
             yerr=[Th_lowerLimit_GR, Th_upperLimit_GR], 
             fmt='o', color='#737373', label='granite', markersize=9)
plt.ylabel('Th (ppm)')
plt.xlabel('U (ppm)')
plt.legend(loc=2, fontsize=12, numpoints=1)


plt.show()