!pmuxi_lim pmuxidelrholim

pmuxidelrholim = ?1

n_ckm = 11

! should lower limit be 0 or -1?
generate cth_r 0 ,, 1 n_ckm
cth_l = 0.9735
cth_ratio = cth_r/cth_l

n_cpv = 11
generate cth_cpv 0 ,, 1 n_cpv



! parameters in quadratic eqn

! square of mass ratio
epsi = `(m_l*m_l/(m_r*m_r))'

! from Herczeg, R = 1 - pmuxidelrholim

! MLRS: zeta^2 + zeta*2*epsi +(2*epsi*epsi -R/2) = 0
b_manif = `(2*epsi)'
c_manif = `(2*epsi*epsi + (pmuxidelrholim - 1.)/2.)'

! NMLRS: ((g_r/g_l)*zeta)^2 + (g_r/g_l)*zeta*2*epsi*(g_r/g_l)*(costh_r/costh_l)*cos(alpha+omega) 
!+(2*[epsi*(g_r/g_l)*(costh_r/costh_l)]^2 - R/2) <~= 0 (eq. 62)
!
! so consider epsi as (g_r**2/g_l**2)*epsi
! and zeta as (g_r/g_l)*zeta
b_nonmanif = `(2*epsi*rckm*cpv)'
c_nonmanif = `(2.*(epsi*rckm)**2  + (pmuxidelrholim - 1.)/2.)'

! constants
m_l = 80.423
lim = pmuxidelrholim


! for non-manifest LR symmetry, with 0<=V_ud^R<=1, using cth_ratio

m_r = [0:1500:5]
zeta_lower = m_r - m_r  ! start with zeros
zeta_upper = m_r - m_r  ! start with zeros

scalar\dummy jj
scalar\dummy kk

do jj=[1:n_ckm]
do kk=[1:n_cpv]

if (exist(`zeta_lower_save')) then destroy zeta_lower_save
zeta_lower_save = m_r - m_r  ! start with zeros

if (exist(`zeta_upper_save')) then destroy zeta_upper_save
zeta_upper_save = m_r - m_r  ! start with zeros

rckm = cth_ratio[jj]
cpv = cth_cpv[kk]

b = eval(b_nonmanif)
c = eval(c_nonmanif)
determ = b*b/4. - c
realroot = where(determ>1.e-10)
norealroot = 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]) 
!zeta_lower_save[realroot] = -b[realroot]/2. - sqrt(determ[realroot]) 
!zeta_upper_save[realroot] = -b[realroot]/2. + sqrt(determ[realroot]) 
!zeta_lower_save[norealroot] = 1.
!zeta_upper_save[norealroot] = -1.

!save the minimum/maximum values of zeta
zeta_upper = max(zeta_upper,zeta_save)
!zeta_lower = min(zeta_lower,zeta_lower_save)
!zeta_upper = max(zeta_upper,zeta_upper_save)

enddo
enddo

m_r_rev = m_r
zeta_lower = -zeta_upper

!m_range = where((zeta_upper - zeta_lower)>1.e-10)
!m_rr = m_r[m_range]
!m_r_rev = m_rr
!zeta_upper = zeta_upper[m_range]
!zeta_lower = zeta_lower[m_range]

sort\down m_r_rev zeta_lower

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

m_low_nonmanif =  min(massvar_nonmanif)

display ` '
display `+++++++++++++++++++++++++++++++++++++++++++++++++++'
display `Minimum non-manifest LRS W_R mass is '//rchar(m_low_nonmanif)
display `for P_mu*xi*delta/rho limit = '//rchar(pmuxidelrholim)
!display `DO NOT TRUST THIS! It is certainly wrong.'
display `+++++++++++++++++++++++++++++++++++++++++++++++++++'
display ` '

destroy zeta_lower zeta_upper

skip1:



! for manifest LR symmetry, with V_ud^R = V_ud^L
m_low_manif = m_l*sqrt(sqrt(2/(1.-lim))) + 1.e-4
display ` '
display `+++++++++++++++++++++++++++++++++++++++++++++++++++'
display `Minimum pseudomanifest LRS W_R mass is '//rchar(m_low_manif)
display `for P_mu*xi*delta/rho limit = '//rchar(pmuxidelrholim)
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