# 
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

etas_m = ( 3 - sqrt(5) )/2
etas_p = ( 3 + sqrt(5) )/2

# A value of eta less than 1:
# 
eta = 0.75
xs_m = ( -eta - sqrt( -(1-3*eta+eta^2) ) )/( 1 - eta )
xs_p = ( -eta + sqrt( -(1-3*eta+eta^2) ) )/( 1 - eta ) 
print( xs_m )
print( xs_p )

# Test the signs in various cases : 
X = -7
( 1 - eta ) * X**2 + 2 * eta * X + ( 1 - 2 * eta ) # positive (H_1)

X = 0
( 1 - eta ) * X**2 + 2 * eta * X + ( 1 - 2 * eta ) # negative (H_0) 

X = 1.
( 1 - eta ) * X**2 + 2 * eta * X + ( 1 - 2 * eta ) # positive (H_1)


# A value of eta greater than 1:
# 
eta = 1.25
xs_m = ( -eta - sqrt( -(1-3*eta+eta^2) ) )/( 1 - eta )
xs_p = ( -eta + sqrt( -(1-3*eta+eta^2) ) )/( 1 - eta ) 
print( xs_m )
print( xs_p )

# Test the signs in various cases : 
X = -7
( 1 - eta ) * X**2 + 2 * eta * X + ( 1 - 2 * eta ) # negative (H_0)

X = 1.0 
( 1 - eta ) * X**2 + 2 * eta * X + ( 1 - 2 * eta ) # positive (H_1) 

X = 10.
( 1 - eta ) * X**2 + 2 * eta * X + ( 1 - 2 * eta ) # negative (H_0) 


# Compute P_F and P_D for 0.381966 < eta < 1.0 :
# 
etas = seq( from=etas_m*(1+0.01), to=1.0*(1-0.01), length.out=100 )

X_m = ( -etas - sqrt( -(1-3*etas+etas^2) ) )/( 1 - etas )
X_p = ( -etas + sqrt( -(1-3*etas+etas^2) ) )/( 1 - etas )

P_F_1 = 1 - atan(X_p)/pi + atan(X_m)/pi
P_D_1 = 1 - atan(X_p-1)/pi + atan(X_m-1)/pi

# Compute P_F and P_D for 1.0 < eta < 2.618034 :
# 
etas = seq( from=1.0*(1+0.01), to=etas_p*(1-0.01), length.out=100 )

X_m = ( -etas - sqrt( -(1-3*etas+etas^2) ) )/( 1 - etas )
X_p = ( -etas + sqrt( -(1-3*etas+etas^2) ) )/( 1 - etas )

P_F_2 = atan(X_m)/pi - atan(X_p)/pi
P_D_2 = atan(X_m - 1 )/pi - atan(X_p - 1)/pi


P_F = c(1.,P_F_1,P_F_2,0.)
P_D = c(1.,P_D_1,P_D_2,0.)

#postscript("../../WriteUp/Graphics/Chapter2/prob_2.2.8.eps", onefile=FALSE, horizontal=FALSE)
plot( P_F, P_D, xlim=c(0,1), ylim=c(0,1), type="l" )
#dev.off()