%pylab inline

def calculate_phase(t, FSS, phase = None):
    """ 
        Calculated the exciton - biexciton phase
    """
    
    if phase == None:
        phase = 0
    
    hbar = 6.56e-16 
    return np.exp((1.0j*(FSS*1e-6*t*1e-9 + phase))/hbar)

def phase_with_dephasing(xlifetime, crosstau, FSS):
    """
        Calculated the phase including dephasing.
    """
    
    t = np.linspace(0, xlifetime, 100)
    sin = []
    phase = 0
    dephasing_event = np.random.exponential(crosstau, 1)[0]    

    for i in xrange(t.size):
    
        p = calculate_phase(t[i], FSS, phase)
        sin.append(p)
        
        if t[i] > dephasing_event:

            if dephasing_event < xlifetime:   
                phase = np.random.random_sample()*1e-9
                dephasing_event += np.random.exponential(crosstau, 1)[0]        
                
    sin = np.array(sin)
    return sin[-1], sin

xtau = 1.
crosstau = 1.
FSS = 2.

xlifetime = np.random.exponential(1, 1)[0]
plt.figure(figsize = (16./1.5, 9./1.5))

t = np.linspace(0, xlifetime, 100)
plt.plot(t, phase_with_dephasing(xlifetime, crosstau, FSS)[1],'bo', markersize = 2)
plt.plot(t, phase_with_dephasing(xlifetime, 1e20, FSS)[1],'r--', markersize = 2)

plt.xlabel('$\\tau_x (ns)$', fontsize = 16) ; plt.ylabel('$Re(\\phi)$', fontsize = 16)
plt.xlim([0, xlifetime]) ; plt.ylim([-1.1, 1.1])
plt.legend(['No decoherence', 'With decoherence'])

test_phases = np.array([phase_with_dephasing(xl, crosstau, 0)[0] for xl in np.random.exponential(xtau, 500)])

real_t_p = np.real(test_phases)
1.0 - float(real_t_p[real_t_p != 1.].size) / real_t_p.size

ghv = 1./(1 + xtau / crosstau)
ghv