"""
Istogrammi con trasparenza
Fit gaussiano
"""

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

def fitpeak(nbins, x1, x2, xs1, xs2, data):
    selected = data[np.flatnonzero( (data>xs1) & (data<xs2) )] #data in the relevant peak
    mu, std = norm.fit(selected)
    plt.figure()
    n_norm, bins_norm, patches = plt.hist(selected, bins=nbins, range=(x1, x2), normed=1) 
    plt.close()
    plt.figure()
    n_all, bins_all, patches = plt.hist(data, bins=nbins, range=(x1, x2))
    plt.close()
    normfactor = np.amax(n_all)/np.amax(n_norm)
    y = mlab.normpdf(bins_norm, mu, std)
    return mu, std, y*normfactor

fff=np.load('data_5c.npz')

mhd_f = fff['mhd_f']   # mhd frequency in time slices around MSE data points
f2 = fff['f2']         # reference freq. aligned with MSE data with Ti>0 and non-empty slices
Fion = fff['Ftor'] + 5./6. * fff['Fsd']

lowm = np.flatnonzero( mhd_f > f2 ) # frequencies higher than the one at q=3/2
highm = np.flatnonzero( mhd_f <= f2 ) # frequencies lower than the one at q=3/2

# SCATTER PLOT OF FREQUENCIES 
fig = plt.figure(facecolor='w')    # plot raw frequencies and separators
plt.scatter(1.e-3*Fion[lowm], 1.e-3*mhd_f[lowm], 20, 'b', marker='o')  # can be picked by mouse
plt.scatter(1.e-3*Fion[highm], 1.e-3*mhd_f[highm], 20, 'r', marker='o')  # can be picked by mouse
plt.plot([0.0, 70.0], [0.0, 70.0], '--k')
plt.xlabel('ion frequency  (kHz)')
plt.ylabel('mhd frequency  (kHz)')
plt.xlim([0.0, 15.0])
plt.ylim([0.0, 30.0])

# HISTOGRAMS OF FREQUENCY RATIO FOR HIGH AND LOW M
nbins = 100   # 
x1 = 0.0    # lower range of the histogram bins
x2 = 5.0    # upper range of the histogram bins
ratios = mhd_f / Fion

plt.figure(facecolor='w')
n_i, bins, patches = plt.hist(ratios[lowm], 
    bins=nbins, range=(x1, x2), 
    alpha=0.5)  # normed=1 
nn, bb, pp = plt.hist(ratios[highm], bins=nbins, range=(x1, x2), facecolor='r', alpha=0.5) 
# Gaussian fit to high-m modes
mu, std, y = fitpeak(nbins, x1, x2, 0.8, 1.2, ratios[highm])
print ('mean and std of frequency ratio', mu, std)
plt.plot(bins, y, '-sm')
plt.xlabel('frequency ratio')
plt.ylabel('Counts')
plt.title('n = 1')

plt.show()