function integral=ps5p2_stage1_integral(x1star); global ALPHAX ALPHAZ ALPHA BETA DELTA ALPHATILDEZ B C X1STAR SELECT OPTIONS... THETA2LOWER THETA2UPPER THETA1STAR Z1STAR ZUNDERBAR ZUPPERBAR; X1STAR=x1star; THETA1STAR=ps5p2_Theta(X1STAR); Z1STAR=ps5p2_z1star; lower=x1star-sqrt( ALPHAX*ALPHAZ/(1+ALPHAZ))*stdn_inv(1-1e-3/2); upper=x1star+sqrt( ALPHAX*ALPHAZ/(1+ALPHAZ))*stdn_inv(1-1e-3/2); ZBOUNDS=sort([ZUNDERBAR Z1STAR ZUPPERBAR]); integral1=quad8('ps5p2_stage1_integrand',lower,ZBOUNDS(1)-1e-5); integral2=quad8('ps5p2_stage1_integrand',ZBOUNDS(1)+1e-5,ZBOUNDS(2)-1e-5); integral3=quad8('ps5p2_stage1_integrand',ZBOUNDS(2)+1e-5,ZBOUNDS(3)-1e-5); integral4=quad8('ps5p2_stage1_integrand',ZBOUNDS(3)+1e-5,upper); integral=integral1+integral2+integral3+integral4;