!pmuxi_lim pmuxilim

pmuxilim = ?1

n_ckm = 11
generate cth_r 0 ,, 1 n_ckm
cth_l = 0.9735
cth_ratio = cth_r/cth_l

! parameters in quadratic eqn

epsi = `(m_l*m_l/(m_r*m_r))'

b_manif = `(epsi)'
c_manif = `(epsi*epsi + (pmuxilim - 1.)/4.)'

b_nonmanif = `(epsi*rckm)'
c_nonmanif = `(0.5*(epsi*rckm)**2 + 0.5*(epsi)**2 + (pmuxilim - 1.)/4.)'

! constants
m_l = 80.423
lim = pmuxilim

! for non-manifest LR symmetry, with 0<=V_ud^R<=1, using cth_ratio
m_low_nonmanif =  m_l*sqrt(sqrt(2/(1.-lim)))
m_r = [m_low_nonmanif:1500:5]

display ` '
display `+++++++++++++++++++++++++++++++++++++++++++++++++++'
display `Minimum non-manifest LRS W_R mass is '//rchar(m_low_nonmanif)
display `for P_mu*xi limit = '//rchar(pmuxilim)
display `+++++++++++++++++++++++++++++++++++++++++++++++++++'
display ` '

zeta_upper = m_r - m_r  ! start with zeros

do jj=[1:n_ckm]

if (exist(`zeta_save')) then destroy zeta_save
zeta_save = m_r - m_r  ! start with zeros
rckm = cth_ratio[jj]

b = eval(b_nonmanif)
c = eval(c_nonmanif)
determ = b*b/4. - c
realroot = where(determ>1.e-10)
! note positive b/2., since rckm can be <0 and we want maximum
zeta_save[realroot] = abs(-b[realroot])/2. + sqrt(determ[realroot]) 

!save the maximum values of zeta
zeta_upper = max(zeta_upper,zeta_save)

enddo

m_r_rev = m_r
zeta_lower = -zeta_upper
sort\down m_r_rev zeta_lower

massvar_nonmanif = m_r_rev//m_r
zetavar_nonmanif = zeta_lower//zeta_upper

destroy zeta_lower zeta_upper



! for manifest LR symmetry, with V_ud^R = V_ud^L
m_low_manif = m_l*sqrt(sqrt(3/(1.-lim))) + 1.e-4
display ` '
display `+++++++++++++++++++++++++++++++++++++++++++++++++++'
display `Minimum pseudomanifest LRS W_R mass is '//rchar(m_low_manif)
display `for P_mu*xi limit = '//rchar(pmuxilim)
display `+++++++++++++++++++++++++++++++++++++++++++++++++++'
display ` '
destroy m_r
m_r = [m_low_manif:1500:5]
b = eval(b_manif)
c = eval(c_manif)
determ = b*b/4. - c
zeta_lower = -b/2. - sqrt(determ)
zeta_upper = -b/2. + sqrt(determ)

m_r_rev = m_r
sort\down m_r_rev zeta_lower

massvar_manif = m_r_rev//m_r
zetavar_manif = zeta_lower//zeta_upper